Insertion-deletion mutation data analysis method, device, system and readable storage medium

Abstract

The present invention provides an insertion and deletion mutation data analysis method, device, system and readable storage medium, which relate to the field of high-throughput second-generation sequencing. The method comprises: obtaining a sequence alignment file and extracting an intersection sequence; generating a target fragment according to the CIGAR value in the intersection sequence; performing a preliminary screening and elimination process on the target fragment; performing a re-alignment on the target fragment to obtain an alignment result; and obtaining mutation site information according to the alignment result. Based on the alignment result, the present invention can accurately find a representative sequence of InDel with high consistency in different sequences by processing the sequence near the InDel breakpoint, quickly obtain the target InDel by performing a local re-alignment, and then use the BAM re-alignment method to achieve accurate frequency statistics of the target InDel, solving the analysis problem of extremely low frequency or complex type mutations.

Description

Technical Field

[0001] The present invention relates to the field of high-throughput second-generation sequencing, and in particular to an insertion / deletion mutation data analysis method, device, system and readable storage medium. Background Art

[0002] Insertions and deletions (InDels) in local regions of the genome are common types of chromosomal mutations in many malignant tumors. They play an important role in the occurrence and development of tumors and often become targets for tumor drugs. For example, in lung cancer, partial sequence deletions in exon 19 of the EGFR gene or skipping mutations in exon 14 of the MET gene are key molecular markers. The accurate identification of these mutations is extremely valuable for clinical diagnosis, treatment selection, and efficacy monitoring. With the development of high-throughput sequencing technology, especially the application of next-generation sequencing (NGS), new possibilities have been provided for the detection of such mutations.

[0003] There are many existing InDel detection algorithms based on second-generation sequencing data, mainly including alignment-based methods, Split Read-based methods, Read Pair-based methods, and haplotype-based detection algorithms. These methods each have their own advantages and scope of application. Alignment-based methods rely on Bayesian models or heuristic algorithms to distinguish between true mutations and alignment errors, but this method is limited by read length and has limited detection capabilities for medium-sized or large InDels. Split Read-based methods can detect larger InDels, but require processing a large amount of reads data, which has a high computational cost. Read Pair-based methods are suitable for detecting large structural variations, but may not be sensitive enough for small InDels. Haplotype-based methods focus on the detection of small InDels, but also face challenges in accuracy and computational efficiency.

[0004] Although there are a variety of detection algorithms and software, their performance in specificity and sensitivity varies, and no algorithm or software can detect all real InDels. In particular, for some complex or low-frequency Small InDels, existing software may miss them during detection, or may not be accurate enough in determining the insertion and deletion breakpoints and estimating the frequency. In addition, the computing speed of the detection algorithm is also a major challenge, especially when dealing with large-scale clinical samples, where fast and accurate analysis results are crucial for clinical decision-making.

[0005] In addition to the limitations of the algorithm itself, the quality of the sample also has an important impact on the accurate detection of InDels. During the sampling, storage, transportation and experiment of tumor tissue or blood specimens, the accuracy of InDel detection may be affected by factors such as insufficient DNA content, oxidative damage or sequencing errors. The existence of these problems makes the rapid, sensitive and accurate detection of potential target InDels a long-standing technical problem.

[0006] In summary, although the existing InDel detection technology has met the needs of clinical and scientific research to a certain extent, there are still many challenges, including but not limited to the accuracy of the algorithm, the efficiency of operation, and the quality control of samples. These are the technical bottlenecks that need to be solved in the field of InDel detection by the current high-throughput sequencing technology. Summary of the invention

[0007] In view of this, in order to solve the above technical problems, the technical solution provided by the present invention includes:

[0008] In a first aspect, the present invention provides an insertion and deletion mutation data analysis method, comprising:

[0009] Obtaining a sequence alignment file corresponding to the target sample, and extracting the intersection sequence existing after the sequence alignment file is aligned with the first bed interval of the preset target capture interval, to obtain an intersection sequence set;

[0010] According to the CIGAR value of each sequence in the intersection sequence set, a target fragment based on the breakpoint extension of each sequence in the intersection sequence set is generated to obtain a first Kmer data set containing the target fragment; each target fragment in the first Kmer data set is preliminarily screened and eliminated to obtain a second Kmer data set;

[0011] Based on the preset targeted capture interval, re-aligning the target fragment to obtain an alignment result;

[0012] Mutation site information is obtained according to the comparison result, and a variation set is obtained based on the mutation site information.

[0013] In an optional embodiment, the method for acquiring the first bed interval of the preset targeted capture interval includes:

[0014] According to the preset targeted capture interval, the preset targeted capture interval is extended to the left and right by a first preset length respectively to obtain the first bed interval; wherein the first preset length represents a maximum InDel length threshold parameter that can be extended;

[0015] The calculation method of the genome start position and the genome end position of the first bed interval is:

[0016] The genomic starting position of the first bed interval is:

[0017] P 基因组起始 =start-L 第一 ;

[0018] The genomic end position of the first bed interval is:

[0019] P 基因组终止 =end+L 第一 ;

[0020] Among them, P 基因组起始 represents the genomic starting position of the first bed interval; P 基因组终止 represents the genomic end position of the first bed interval; start represents the genomic starting position of the preset targeted capture interval; end represents the genomic end position of the preset targeted capture interval; L 第一 Represents the first preset length.

[0021] In an optional embodiment, the obtaining of a sequence alignment file corresponding to the target sample, and extracting an intersection sequence existing after the sequence alignment file is aligned with the first bed interval of the preset target capture interval to obtain an intersection sequence set, comprises:

[0022] Taking the sequence alignment file as input, searching for a sequence matching the sequence alignment file in the first bed interval as a matching sequence;

[0023] Performing filtering and de-duplication processing on the matching sequence to obtain the intersection sequence;

[0024] Determining whether there is a sequence in the intersection sequence set that meets the correction condition;

[0025] If yes, the sequence satisfying the correction condition is used as the sequence to be corrected; and the sequences other than the sequence to be corrected in the intersection sequence set are used as the non-corrected sequence set;

[0026] The sequence to be corrected is subjected to CIGAR value correction processing to obtain a corrected sequence set; and the corrected sequence set and the non-corrected sequence set are taken together as the corrected intersection sequence set; and the correction condition is:

[0027] Chimeric alignment sequences that simultaneously meet the following conditions: A. with SA tag; B. FLAG value does not contain 0×100; C. aligned to the same chromosome; D. same chain type. In an optional embodiment, the CIGAR value correction process includes:

[0028] The read segments of the sequence to be corrected are classified based on the CIGAR value to obtain a classification result; based on the classification result, according to the linear alignment result of the chimeric alignment sequence, the CIGAR value and position of the representative alignment in the chimeric alignment sequence are corrected to obtain the corrected sequence set.

[0029] In an optional embodiment, the classification results include MS class and SM class;

[0030] In the MS class, the left side is the successfully mapped fragment, and the right side is the Softclip fragment;

[0031] In the SM class, the left side is the Softclip segment, and the right side is the successfully mapped segment.

[0032] In an optional embodiment, the linear alignment result includes a representative alignment and a substitute alignment; wherein the representative alignment is an alignment whose FLAG value of the linear alignment result does not include 0×100 and is the final alignment of the linear alignment result; the substitute alignment is an alignment whose FLAG value of the linear alignment result includes 0×100 and is the alignment corresponding to the SA tag;

[0033] The method of correcting the CIGAR value and position of the representative alignment in the chimeric alignment sequence based on the classification result according to the linear alignment result of the chimeric alignment sequence comprises:

[0034] Calculate the theoretical length of the target fragment;

[0035] According to the theoretical length of the target fragment and the actual length of the chimeric comparison sequence, confirm the gene variation type and variation length of the chimeric comparison sequence; wherein the variation length includes the deletion length and the insertion length;

[0036] Based on the gene variation type and variation length, the read length of the chimeric alignment sequence, the theoretical length of the target fragment, the length and position of the Match fragment corresponding to the representative alignment and the substitute alignment, the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence are modified.

[0037] In an optional embodiment, the calculation formula for the theoretical length of the target fragment is:

[0038] L 目标 =SM 终止 -MS 起始 ;

[0039] Among them, L 目标 represents the theoretical length of the target fragment; SM 终止 Represents the end position of the comparison of the SM class; MS 起始 represents the alignment start position of the MS class;

[0040] The step of confirming the gene variation type and variation length of the chimeric comparison sequence according to the theoretical length of the target fragment and the actual length of the chimeric comparison sequence comprises:

[0041] If the theoretical length of the target fragment is greater than the actual length of the chimeric comparison sequence, the gene variation type of the chimeric comparison sequence is determined to be deletion; the calculation formula for the deletion length of the chimeric comparison sequence whose gene variation type is deletion is: L 缺失 =L 目标 -L 序列 Among them, L 缺失 The deletion length of the chimeric alignment sequence representing the gene variation type as deletion; L 目标 represents the theoretical length of the target fragment; L 序列 is the actual length of the chimeric aligned sequence;

[0042] If the theoretical length of the target fragment is less than the actual length of the chimeric comparison sequence, the gene variation type of the chimeric comparison sequence is determined to be insertion. The calculation formula for the insertion length of the chimeric comparison sequence with the gene variation type being insertion is: L 插入 =L 序列 -L 目标 Among them, L 插入 The insertion length of the chimeric alignment sequence representing the genetic variation type as insertion.

[0043] In an optional embodiment, the modifying of the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence based on the gene variation type and variation length, the read length of the chimeric alignment sequence, the theoretical length of the target fragment, the length and position of the Match fragment corresponding to the representative alignment and the substitute alignment comprises:

[0044] If the classification result of the representative alignment sequence in the chimeric alignment sequence is the MS class, the starting position remains unchanged; if the variation type is deletion, the CIGAR value is modified from L1ML2S to L1ML 缺失 D(L 序列 -L1)M; if the mutation type is insertion, the CIGAR value is changed from L1ML2S to L1ML 插入 I(L 目标 -L1)M;

[0045] If the classification result of the representative alignment sequence in the chimeric alignment sequence is the SM class, the starting position is replaced with the position of the alternative alignment MS type; if the variation type is deletion, the CIGAR value is modified from L2SL1M to L3ML 缺失 D(L序列 -L3)M; if the mutation type is insertion, the CIGAR value is changed from L2SL1M to L3ML 插入 I(L 目标 -L3)M;

[0046] Wherein, L1 represents the length of the successfully mapped segment; L2 represents the length of the Softclip segment; L3 represents the length of the successfully mapped segment of the substitute comparison; L 缺失 is the missing length; L 插入 is the insertion length; L 序列 represents the actual length of the chimeric alignment sequence; L 目标 Represents the theoretical length of the target fragment.

[0047] In an optional embodiment, generating a target fragment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set includes:

[0048] In the intersection sequence set, the sequences whose screening gene variation type is deletion or insertion and whose deletion or insertion length is not less than 3 are used as screening sequences;

[0049] At the breakpoint position of the insertion or deletion of the screening sequence, a second preset length is extended to the left and right to serve as the target fragment.

[0050] In an optional embodiment, the performing preliminary screening and elimination processing on each of the target fragments in the first Kmer data set includes:

[0051] In the first Kmer data set, a preliminary screening and elimination process is performed according to at least one of the following elimination rules:

[0052] A. Eliminate the target fragments whose mapping quality of the sequence where the target fragment is located is less than a preset mapping quality threshold;

[0053] B. Based on the base composition, remove the target fragments with continuous repeated bases of not less than 12 bp;

[0054] C. Determine the target fragment to be removed according to the comparison information of the sequence where the target fragment is located.

[0055] In an optional embodiment, determining the target fragment to be removed according to the comparison information of the sequence where the target fragment is located includes:

[0056] Obtain the specific location of the sequence breakpoint where the target fragment is located in the reference genome;

[0057] According to the gene variation type and the specific position, the start and end positions of the breakpoint are calculated respectively; the calculation method includes:

[0058] If the gene variation type is insertion, the start and end positions of the breakpoint include:

[0059] Starting position: P 插入起始 =P 具体 -L 插入长度 ;

[0060] End position: P 插入终止 =P 具体 +L 插入长度 ;

[0061] If the gene variation type is deletion, the start and end positions of the breakpoints include:

[0062] Starting position: P 缺失起始 =P 具体 ;

[0063] End position: P 缺失终止 =P 具体 +L 缺失长度 ;

[0064] Among them, P 插入起始 Represents the starting position of the breakpoint where the gene variation type is insertion; P 插入终止 represents the end position of the breakpoint where the gene variation type is insertion; P 缺失起始 Represents the starting position of the breakpoint where the gene variation type is deletion; P 缺失终止 Represents the end position of the breakpoint where the gene variation type is deletion; P 具体 Represents the specific location; L 插入长度 represents the length of the inserted sequence; L 缺失长度 represents the length of the missing sequence;

[0065] If the interval where the start and end positions of the breakpoint are located does not intersect with the preset target capture interval, the target segment corresponding to the breakpoint is removed.

[0066] In an optional embodiment, after performing preliminary screening and elimination processing on each of the target fragments in the first Kmer data set to obtain the second Kmer data set, the method further includes:

[0067] Based on the second Kmer data set, according to the length of the inserted or deleted sequence in each of the target fragments, the sequence is extended to the breakpoints at both ends to generate Extend fragments, and a first Extend fragment data set is constructed;

[0068] In the first Extend segment data set, one or two representative Extend segments are obtained for each bed interval as target segments after screening of the interval, and the second Extend segment data set containing the target segments is generated.

[0069] In an optional embodiment, the method for calculating the extension length extending to the breakpoints at both ends of the sequence includes:

[0070] If the gene variation type corresponding to the target fragment is insertion, the length of the target fragment extending to the left and right is:

[0071] L 延伸 =C 可调 +L I ×2;

[0072] If the gene variation type corresponding to the target segment is deletion, the length of the target segment extending to the left and right is:

[0073] L 延伸 =C 可调 ;

[0074] Among them, L 延伸 Represents the length of the target segment extending to the left and right, and does not include the Softclip area; C 可调 Represents preset adjustable parameters; L I Represents the length of the inserted fragment in the CIGAR value.

[0075] In an optional implementation, in the first Extend segment data set, one or two representative Extend segments are obtained for each bed interval as target segments after screening of the interval, and the second Extend segment data set containing the target segments is generated, including:

[0076] The first Extend fragment data set is classified according to the length of the insertion or deletion sequence; wherein the sequences with the length of the insertion or deletion sequence not less than 6 are regarded as the E1a data set; and the sequences with the length of the insertion or deletion sequence less than 6 are regarded as the E1b data set.

[0077] The Extend fragments in the classified E1a dataset and the E1b dataset are sorted in order according to the number of Extend fragments, the number of target fragments, the quality of sequence alignment, and the length of the inserted or deleted sequence, and the first-ranked Extend fragment in the E1a dataset and the E1b dataset is selected as the best representative Extend fragment;

[0078] Comparing the similarity and the difference in the number of Extend segments between the longest common subsequences of the optimal representative Extend segments of the E1a data set and the E1b data set, and obtaining a comparison result; including:

[0079] A. If the similarity is greater than 80%, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set;

[0080] B. If the similarity is not greater than 80%, and the number of Extend segments is less than 3, and the difference in the number of Extend segments is more than 30 times, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set;

[0081] C. If both A and B do not match, the comparison result is: both the E1a data set and the E1b data set are retained to generate the second Extend segment data set.

[0082] In an optional embodiment, the re-aligning the target fragment based on the preset target capture interval to obtain an alignment result includes:

[0083] Extending the preset targeted capture interval to the left and right by a third preset length respectively to obtain a second bed interval;

[0084] Realigning the target fragment with the second bed interval to obtain the alignment result;

[0085] The scoring rule expression in the re-alignment algorithm is:

[0086] [match:1,mismatch:1,gap open:2,gap extend:0];

[0087] Among them, Match: 1 means that if the two aligned bases are the same, 1 point is obtained; Mismatch: -1 means that if the two aligned bases are different, 1 point is deducted; Gap Open: -2 means that if a gap is introduced, 2 points are deducted; Gap Extend: 0 means that if a gap already exists, no point is deducted.

[0088] In a second aspect, the present invention provides an insertion and deletion mutation data analysis device, comprising:

[0089] An alignment intersection module is used to obtain a sequence alignment file corresponding to a target sample, and extract an intersection sequence existing after the sequence alignment file is aligned with the first bed interval of a preset target capture interval, to obtain an intersection sequence set;

[0090] A preliminary screening module is used to generate a target fragment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set, so as to obtain a first Kmer data set containing the target fragment; and perform preliminary screening and elimination processing on each target fragment in the first Kmer data set to obtain a second Kmer data set;

[0091] A re-alignment module, used for re-aligning the target fragment based on the preset target capture interval to obtain an alignment result;

[0092] A result generation module is used to obtain mutation site information according to the comparison result, and obtain a variation set based on the mutation site information.

[0093] In a third aspect, the present invention provides an insertion-deletion mutation data analysis system, comprising a memory and a processor, wherein the memory stores an insertion-deletion mutation data analysis program, and the processor runs the insertion-deletion mutation data analysis program to enable the insertion-deletion mutation data analysis system to perform the insertion-deletion mutation data analysis method as described in any one of the aforementioned embodiments.

[0094] In a fourth aspect, the present invention provides a computer-readable storage medium, on which an insertion-deletion mutation data analysis program is stored. When the insertion-deletion mutation data analysis program is executed by a processor, the insertion-deletion mutation data analysis method as described in any one of the aforementioned embodiments is implemented.

[0095] The present invention provides an insertion and deletion mutation data analysis method, device, system and readable storage medium. Among them, the insertion and deletion mutation data analysis method significantly improves the efficiency, accuracy and speed of mutation detection. Specifically, the beneficial effects of the method include the following points:

[0096] (1) The present invention extracts the sequence alignment file of the target sample, generates the target fragment based on breakpoint extension and constructs the Kmer data set, which can quickly focus on the sequence region that may contain InDel mutations, thereby improving the detection efficiency of complex mutations.

[0097] (2) After obtaining the intersection sequence set, the present invention further achieves accurate frequency statistics of the target InDel by local realignment and acquisition of mutation site information. This method not only relies on local realignment of representative sequences, but also uses the original BAM file for realignment, ensuring that the statistics of mutation frequency are both fast and accurate. This is particularly important for extremely low-frequency or complex types of mutations, because these mutations are often difficult to accurately detect using traditional methods.

[0098] (3) The present invention avoids comprehensive local realignment of a large number of sequences, but instead reduces the amount of calculation and analysis time by screening and focusing on representative sequences. This method not only increases the speed of analysis, but also enables more samples to be processed with limited resources, which is particularly important for high-throughput sequencing data analysis.

[0099] (4) The present invention can accurately find representative InDel sequences with high consistency in different sequences, solving the problem of analyzing extremely low-frequency or complex mutations. This has important application value in the fields of studying tumor heterogeneity, tracking disease progression, and personalized medicine.

[0100] In summary, the method provided by the present invention achieves efficient, accurate and rapid analysis of insertion and deletion mutations through its unique steps and processing methods, and has significant beneficial effects. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] In order to more clearly illustrate the technical solution of the present invention, the following is a brief introduction to the drawings required for use in the embodiments. It should be understood that the following drawings only illustrate certain embodiments of the present invention and should not be regarded as limiting the scope of protection of the present invention. In each of the drawings, similar components are numbered similarly.

[0102] Figure 1 A schematic diagram of the structure of the hardware operating environment involved in the insertion / deletion mutation data analysis method embodiment of the present invention;

[0103] Figure 2 This is a schematic diagram of the process of Example 1 of the insertion / deletion mutation data analysis method of the present invention;

[0104] Figure 3 This is a schematic diagram of a detailed process of step S100 in Example 2 of the insertion / deletion mutation data analysis method of the present invention;

[0105] Figure 4 This is a schematic diagram of a detailed process of step S150 in Example 2 of the insertion / deletion mutation data analysis method of the present invention;

[0106] Figure 5 This is a schematic diagram of a detailed process of step S152 in Example 2 of the insertion / deletion mutation data analysis method of the present invention;

[0107] Figure 6 This is a schematic diagram of a detailed process of step S200 in Example 3 of the insertion / deletion mutation data analysis method of the present invention;

[0108] Figure 7 This is a schematic diagram of a detailed process of step S211 in Example 3 of the insertion / deletion mutation data analysis method of the present invention;

[0109] Figure 8 This is a schematic diagram of the process of Example 4 of the insertion and deletion mutation data analysis method of the present invention;

[0110] Fig. 9 This is a schematic diagram of a detailed process of step S600 in Example 4 of the insertion / deletion mutation data analysis method of the present invention;

[0111] Fig.10 This is a schematic diagram of a detailed process of step S300 in Example 5 of the insertion / deletion mutation data analysis method of the present invention;

[0112] Fig.11 This is an example diagram of the CIGAR value correction of the chimeric alignment sequence in Example 6 of the insertion and deletion mutation data analysis method of the present invention;

[0113] Fig.12 This is a specific sequence example diagram used for demonstration in Example 6 of the insertion / deletion mutation data analysis method of the present invention;

[0114] Fig.13 This is an example diagram of the generation of Kmer fragments (target fragments in the first Kmer data set) and Extend fragments (target fragments in the second Extend fragment data set) in Example 6 of the insertion and deletion mutation data analysis method of the present invention;

[0115] Fig.14 Schematic diagram of module connection of a device for processing curled images of intestinal spatiotemporal groups in an embodiment of the present invention. DETAILED DESCRIPTION

[0116] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0117] The components of the embodiments of the present invention generally described and shown in the drawings herein may be arranged and designed in a variety of different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative work are within the scope of protection of the present invention.

[0118] Hereinafter, the terms "including", "having" and their cognates, which may be used in various embodiments of the present invention, are intended only to indicate specific features, numbers, steps, operations, elements, components or combinations of the foregoing items, and should not be understood as first excluding the existence of one or more other features, numbers, steps, operations, elements, components or combinations of the foregoing items or adding the possibility of one or more features, numbers, steps, operations, elements, components or combinations of the foregoing items.

[0119] Furthermore, the terms “first”, “second”, “third”, etc. are merely used for distinguishing descriptions and are not to be understood as indicating or implying relative importance.

[0120] Unless otherwise defined, all terms (including technical terms and scientific terms) used herein have the same meanings as those generally understood by those skilled in the art to which the various embodiments of the present invention belong. The terms (such as those defined in generally used dictionaries) will be interpreted as having the same meanings as the contextual meanings in the relevant technical field and will not be interpreted as having idealized meanings or overly formal meanings unless clearly defined in the various embodiments of the present invention.

[0121] It should be understood that the specific embodiments described herein are only used to explain the present invention, and are not used to limit the present invention.

[0122] like Figure 1 , which is a schematic diagram of the structure of the hardware operating environment of the terminal involved in the embodiment of the present invention.

[0123] The insertion and deletion mutation data analysis system of the embodiment of the present invention can be a PC, or a mobile terminal device such as a smart phone, a tablet computer or a portable computer. The insertion and deletion mutation data analysis system may include: a processor 1001, such as a CPU, a network interface 1004, a user interface 1003, a memory 1005 and a communication bus 1002. Among them, the communication bus 1002 is used to realize the connection and communication between these components. The user interface 1003 may include a display screen, an input unit such as a keyboard, and a remote control. The optional user interface 1003 may also include a standard wired interface and a wireless interface. The network interface 1004 may optionally include a standard wired interface and a wireless interface (such as a WI-FI interface). The memory 1005 may be a high-speed RAM memory, or a stable memory, such as a disk memory. The memory 1005 may also be a storage device independent of the aforementioned processor 1001.

[0124] Those skilled in the art will understand that Figure 1 The insertion / deletion mutation data analysis system shown in the figure does not constitute a limitation thereof, and may include more or fewer components than shown in the figure, or combine certain components, or arrange the components differently. Figure 1 As shown, the memory 1005 as a computer-readable storage medium may include an operating system, a data interface control program, a network connection program, and an insertion and deletion mutation data analysis program.

[0125] In summary, the insertion and deletion mutation data analysis method provided by the present invention significantly improves the efficiency and accuracy of complex mutation detection by extracting and correcting chimeric alignment sequences that meet the correction conditions in the sequence alignment file of the target sample, and constructing a Kmer data set based on breakpoint extension. This method achieves rapid and accurate statistics of the target InDel mutation frequency through local realignment and realignment using the original BAM file, effectively solving the analysis problem of extremely low frequency or complex type mutations, and especially showing significant advantages in high-throughput sequencing data analysis.

[0126] The corresponding relationship between the names of multiple data sets and intervals and their English abbreviations in the embodiment is shown in the following table:

[0127] Table 1. Correspondence between data set and interval names and symbols in the embodiment

[0128] name symbol name symbol Preset target capture interval B0 The first Extend segment dataset E1 The first Kmer dataset K1 The second Extend segment data set E2 The second Kmer dataset K2 Dataset of sequences with insertion or deletion length no less than 6 E1a Intersection Sequence Set R1 Dataset of sequences with insertion or deletion length less than 6 E1b Corrected sequence set R2 First bed area B1 Uncorrected Sequence Set R3 Second bed area B2

[0129] Embodiment 1:

[0130] Reference Figure 2 Embodiment 1 of the present invention provides an insertion and deletion mutation data analysis method, comprising:

[0131] Step S100, obtaining a sequence alignment file corresponding to the target sample, and extracting the intersection sequence existing after the sequence alignment file is aligned with the first bed interval of the preset target capture interval to obtain an intersection sequence set;

[0132] As mentioned above, the target samples may include but are not limited to tumor samples that need to be tested, and may also be other biochemical samples, and their sources may also be blood, tissue, lavage fluid, urine, saliva, feces, pleural effusion, cerebrospinal fluid, etc.

[0133] The above-mentioned sequence alignment file can be a sequence alignment file in BAM format; it should be noted that the sequence alignment file should be a file that has been aligned and sorted by alignment software. Among them, BAM format is the standard format of data after data preprocessing and alignment of the second-generation sequencing.

[0134] The preset target capture interval can be a target capture bed interval of a bed format file of a genome target region. In this embodiment, the preset target capture interval can be defined as B0.

[0135] It should be noted that the targeted capture bed interval B0 refers to a specific genomic region selected by the researcher according to the research purpose in high-throughput sequencing analysis. These regions may include, but are not limited to, known genes, disease-related regions, or any other biologically interesting DNA fragments. In this embodiment, the B0 interval is defined in the BAM format file after alignment and sorting by the alignment software, which determines the specific position in the sequencing data that needs further analysis.

[0136] As mentioned above, the first bed interval can be a new interval obtained by extending the original input targeted capture bed interval B0. This extension operation is to ensure that even if the insertion and deletion (InDel) event occurs near the boundary of the original targeted capture interval during the analysis process, it can be included, thereby improving the comprehensiveness and accuracy of variation detection. In this embodiment, the first bed interval is defined as B1.

[0137] As mentioned above, step S100 involves obtaining sequence alignment files of the target sample, which files may be the result of aligning high throughput sequencing (HTS) data with a reference genome.

[0138] Step S200, generating a target segment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set, so as to obtain a first Kmer data set containing the target segment; performing a preliminary screening and elimination process on each target segment in the first Kmer data set, so as to obtain a second Kmer data set;

[0139] In this embodiment, the CIGAR value refers to a string used to describe the alignment relationship between each sequencing read and the reference genome in a sequence alignment file in the BAM format. The CIGAR string provides accurate alignment information for the read, including how the bases in the read match, insert, delete, or soft clip with the bases on the reference genome.

[0140] Based on the intersection sequence set obtained in step S100, this step generates target segments based on the extension of each sequence breakpoint to form a first Kmer data set. Subsequently, these Kmers are preliminarily screened to remove segments that do not meet the conditions to obtain a second Kmer data set.

[0141] In this step, a first Kmer data set is generated after preliminary screening, which contains Kmer fragments that may represent real InDel events. The screening process helps to eliminate possible false positives, so that subsequent analysis is more focused on real variations. In this embodiment, the Kmer fragment can be defined as the target fragment, and the first Kmer data set is defined as Kmer data set K1, and the second Kmer data set is defined as Kmer data set K2. It is used to characterize each breakpoint in the Kmer data set.

[0142] Step S300: re-align the target fragment based on the preset target capture interval to obtain an alignment result.

[0143] The preset targeted capture interval B0 or ​​its extended bed interval is re-aligned with the target fragment to obtain an alignment result.

[0144] Specifically, the preset target capture interval B0 may be extended to create a larger analysis interval.

[0145] As described above, based on the bed interval B0 of the preset target capture interval obtained in step S100, the target fragment is locally realigned to accurately determine the location and nature of the InDel event, thereby obtaining the result of local realignment and revealing detailed InDel information.

[0146] Local realignment in this step can improve the accuracy of variant detection, especially in complex or repetitive sequence regions.

[0147] It should be noted that the local realignment in this embodiment refers to the process of recalibrating those regions that may have insertion or deletion (InDel) mutations during the genome alignment process. The purpose of local realignment is to improve the accuracy of the alignment results through a more accurate alignment algorithm (such as the Smith-Waterman algorithm), thereby reducing erroneous variation detection.

[0148] Step S400, obtaining mutation site information according to the comparison result, and obtaining a variation set based on the mutation site information.

[0149] In the above, according to the comparison result in step S300, the mutation site information is identified, and a variation set containing all identified InDel events is generated accordingly, thereby obtaining a comprehensive variation set, which provides basic data for further biological analysis and research.

[0150] This step provides a comprehensive InDel variation catalog, which helps to understand the genetic variation characteristics of the sample.

[0151] For example, the indelpost standardized alignment results can be used to obtain accurate mutation site information, perform mutation depth statistics and mutation frequency calculations, and then collect standardized mutation information and mutation frequencies to obtain the final mutation set.

[0152] In summary, in this embodiment, based on the alignment results, by performing a series of processing on the sequences near the InDel breakpoints, representative InDel sequences with high consistency in different sequences can be accurately found, and then the target InDel can be quickly obtained by only performing local realignment on the representative sequences, and then the original BAM realignment method is used to achieve accurate frequency statistics of the target InDel. The above method solves the problem of analyzing extremely low frequency or complex types of mutations (such as point mutations and insertion and deletion monomer mutations, insertion and deletion combined monomer mutations).

[0153] Embodiment 2:

[0154] Embodiment 2 of the present invention provides an insertion and deletion mutation data analysis method. Based on the above embodiment 1, before step S100, the method may further include: obtaining a first bed interval of a preset target capture interval;

[0155] The method for obtaining the first bed interval of the preset target capture interval includes:

[0156] According to the preset targeted capture interval, the preset targeted capture interval is extended to the left and right by a first preset length respectively to obtain the first bed interval; wherein the first preset length represents a maximum InDel length threshold parameter that can be extended;

[0157] The calculation method of the genome start position and the genome end position of the first bed interval is:

[0158] The genomic starting position of the first bed interval is:

[0159] P 基因组起始 =start-L 第一 ;

[0160] The genomic end position of the first bed interval is:

[0161] P 基因组终止 =end+L 第一 ;

[0162] Among them, P 基因组起始 represents the genomic starting position of the first bed interval; P 基因组终止represents the genomic end position of the first bed interval; start represents the genomic starting position of the preset targeted capture interval; end represents the genomic end position of the preset targeted capture interval; L 第一 Represents the first preset length.

[0163] In the above, obtaining the first bed interval of the preset target capture interval involves determining the starting point of the analysis, that is, the preset target capture interval. This interval is predefined based on the research purpose and experimental design and is used to specify the genomic region to be analyzed.

[0164] In this step, the processing of extending the first preset length to the left and right is performed, that is, after determining the preset target capture interval, the start and end positions of the interval are extended. The extension length is preset and represents the possible maximum InDel length threshold parameter. In this embodiment, it is defined as max_del_len, which can be L in the formula 第一 This parameter is intended to ensure that even InDel events occurring near the boundaries of the targeted capture interval are included in the analysis.

[0165] As described above, the genomic starting position and ending position of the first bed interval are calculated: through mathematical calculation, the starting position of the preset target capture interval is extended to the left by max_del_len, and the ending position is extended to the right by max_del_len, thereby obtaining the extended first bed interval.

[0166] In this embodiment, the original input bed interval B0 can first be defined as [chrN:start-end], where chrN is the chromosome name, start is the genomic starting position of the preset target capture interval, and end is the genomic ending position of the preset target capture interval.

[0167] Then, define the length of the left and right extensions of max_del_len respectively. Where max_del_len is the adjustable maximum InDel length threshold parameter, and the expression for the extended bed interval B1 can be:

[0168] [chrN:start-max_del_len~end+max_del_len].

[0169] The advantage of this method is that it can expand the scope of analysis and increase the possibility of detecting InDel events near the boundary. By extending the operation, potential variable regions can be more comprehensively covered, thereby improving the accuracy and reliability of data analysis.

[0170] For further reference, Figure 3(The start and end in the figure are the start and end of the method for obtaining the intersection sequence set in this embodiment, respectively), the step S100, obtaining the sequence alignment file corresponding to the target sample, and extracting the intersection sequence existing after the sequence alignment file is aligned with the first bed interval of the preset target capture interval, to obtain the intersection sequence set, including:

[0171] Step S110, taking the sequence alignment file as input, searching in the first bed interval for a sequence that matches the sequence alignment file as a matching sequence;

[0172] The purpose of the above steps is to extract the sequences that match the preset target capture interval (first bed interval) from the sequence alignment file. This step can be implemented by the fetch function in the pysam.AlignmentFile class, so that the result obtained by matching is the matching sequence, that is, the reads that overlap with the preset target capture interval.

[0173] This step allows focusing on specific genomic regions, increasing the specificity and efficiency of subsequent analyses.

[0174] Step S120, filtering and deduplicating the matching sequence to obtain the intersection sequence;

[0175] It should be noted that in high-throughput sequencing, the same DNA molecule may be sequenced multiple times, resulting in duplicate reads in the comparison results. The purpose of this step is to perform deduplication on the matching sequence to ensure the accuracy of subsequent analysis. The result obtained is a qualified intersection sequence after deduplication, that is, a sequence set after excluding duplicate reads. In this embodiment, the intersection sequence can be defined as sequence set R1.

[0176] Deduplication helps reduce redundancy in data processing, improves computational efficiency, and avoids over-weighting duplicate reads during variant detection.

[0177] Step S130, determining whether there is a sequence satisfying a correction condition in the intersection sequence set;

[0178] Next, the CIGAR value correction process is performed on the sequences that meet the correction condition, thereby obtaining a corrected sequence set, which provides a more accurate data basis for subsequent analysis. In this embodiment, the intersection sequence can be defined as the sequence set R1.

[0179] It should be noted that the sequences that meet the correction conditions here are all Softclip sequences, which refer to those sequence fragments in the BAM format sequence alignment file of high-throughput sequencing data that are not aligned with the reference genome and are marked as soft clips. These Softclips usually occur at both ends of the read, especially when the read covers the insertion or deletion (InDel) variation region in the genome.

[0180] It should be noted that the matching sequences are filtered and deduplicated to obtain the intersection sequence R1. If step S200 is directly performed to generate a target fragment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set to obtain a first Kmer data set containing the target fragment, the true InDel variation in the Softclip sequence may be missed, resulting in incomplete detection results; at the same time, an erroneous CIGAR value may cause normal sequence changes to be misjudged as InDel variations, increasing the risk of false positives; in addition, an inaccurate CIGAR value may affect subsequent variation detection and analysis, such as mutation frequency calculation and variation biological interpretation, reduce the reliability and credibility of the detection, and have a negative impact on clinical decisions and research results.

[0181] Therefore, in this embodiment, CIGAR value correction processing is required, and the processing method of the correction processing is crucial to ensure the quality and effectiveness of high-throughput sequencing data analysis.

[0182] In the above step S130, it is necessary to check whether the sequences in the intersection sequence set meet specific correction conditions, which may involve specific features of the sequence, such as Softclip sequence, so as to determine which sequences need CIGAR value correction. Sequences that may affect the analysis results are processed in a targeted manner to improve data quality.

[0183] Step S140: If yes, the sequence satisfying the correction condition is used as a sequence to be corrected; and the sequences other than the sequence to be corrected in the intersection sequence set are used as a non-corrected sequence set;

[0184] Step S150, performing CIGAR value correction processing on the sequence to be corrected to obtain a corrected sequence set; and the corrected sequence set and the non-corrected sequence set are used together as the corrected intersection sequence set; the correction condition is: a chimeric alignment sequence that simultaneously meets the following conditions: A. with SA tag; B. FLAG value does not contain 0×100; C. Aligned to the same chromosome; D. The chain type is the same.

[0185] The correction condition refers to a special type of sequence alignment encountered in high-throughput sequencing data analysis, especially when processing sequence alignment files in BAM format, which can be called Softclip. Softclip means that a part of the sequence in the sequence alignment does not fully match the reference genome, but is "softly clipped" to the end of the reference sequence, usually indicating a possible insertion or deletion mutation.

[0186] This correction condition contains an SA tag: SA, or Secondary Alignment, refers to the fact that this sequence has multiple matching regions on the genome. It can be that the same part of the sequence has different matching regions, or that different regions on a sequence have different matching regions. The present invention is defined as Softclip sequences with only two different matching regions on the genome in different regions. In the sequence alignment with an SA tag, the SA tag provides additional alignment information. The FLAG value does not contain 0×100: In the sequence alignment file in the BAM format, the FLAG value is used to indicate the characteristics and status of the alignment. Each alignment has a corresponding FLAG value, which contains multiple flags to describe the different properties of the alignment. The hexadecimal flag 0×100 here represents not primary alignment, that is, non-primary alignment. If the FLAG value does not contain 0×100, it means that the alignment is a primary alignment. Aligned to the same chromosome: It means that these different Softclip sequences must be aligned to the same chromosome, which means that they are all located on the same chromosome of the reference genome. Same strand type: refers to the strand type (strand) refers to the positive strand (Watson strand) or the negative strand (Crick strand) of DNA. This condition requires that these Softclip sequences must be located on the same strand type, that is, they are all located on the positive strand or all located on the negative strand. Softclip sequences that meet these conditions at the same time are described as chimeric alignment sequences.

[0187] When these conditions are met, it means that these Softclip sequences may be caused by real biological variations (such as insertions or deletions) rather than sequencing or alignment errors. Therefore, these sequences need special treatment to ensure that these variants can be accurately identified and quantified in subsequent variant detection analysis.

[0188] In actual data processing, this condition is used to identify sequences that may need further analysis and correction. By correcting the CIGAR values ​​of these sequences, the location and type of insertion or deletion mutations can be more accurately determined, thereby improving the accuracy of variant detection.

[0189] For sequences that meet the correction conditions, their CIGAR values ​​need to be adjusted to correct possible alignment errors. For the remaining sequences that do not meet the correction conditions, they are directly retained, thereby obtaining a corrected sequence set and a non-corrected sequence, forming a new corrected intersection sequence set. This improves the accuracy of InDel detection, especially in the Softclip area, and reduces false positive or false negative results caused by alignment errors.

[0190] In addition, if there is no sequence satisfying the correction condition in the intersection sequence set, no CIGAR value correction processing is performed.

[0191] In this embodiment, by analyzing and adjusting the CIGAR values ​​of the sequences that meet the correction conditions, the alignment results can be corrected, thereby providing a more accurate data basis for subsequent mutation detection and analysis. This step is a key step in the insertion and deletion mutation data analysis method and helps to improve the accuracy and reliability of mutation detection.

[0192] For further reference, Figure 4 In step S150, the CIGAR value correction process includes:

[0193] Step S151 , classifying the read segments of the sequence to be corrected based on the CIGAR values ​​to obtain classification results.

[0194] The purpose of CIGAR value correction is to correct possible alignment errors, reduce the identification of false positive variants, and thus improve the reliability of variant detection. Through this correction, we can obtain a corrected sequence set, which helps to more accurately reflect the actual variation of the sequence, thereby improving the accuracy of subsequent variant detection.

[0195] In the alignment results, some reads may match multiple positions of the reference genome, forming chimeric alignments. The CIGAR values ​​of these sequences can reveal different types of variants, such as insertions or deletions. By analyzing the CIGAR values, we can classify these chimeric alignment sequences, which helps to identify and process those sequences that may contain complex variants, thereby improving the accuracy and sensitivity of variant detection.

[0196] In one embodiment, the classification results include SM and MS categories;

[0197] In the MS class, the left side is a successfully mapped segment, and the right side is a Softclip segment; for example, 40S110M or 60S40M50S.

[0198] In the SM class, the left side is the Softclip segment, and the right side is the successfully mapped segment, such as 40M110S or 50S40M60S.

[0199] Step S152, based on the classification result and the linear alignment result of the chimeric alignment sequence, correct the CIGAR value and position of the representative alignment in the chimeric alignment sequence to obtain the corrected sequence set.

[0200] As described above, based on the classification, this step involves correcting the CIGAR values ​​and alignment positions of the representative alignments in the chimeric alignment sequence, thereby obtaining a corrected sequence set, which includes the corrected CIGAR values ​​and alignment positions. In this embodiment, the corrected sequence set R2 and the non-corrected sequence set R3 can be defined, and the intersection sequence set after the CIGAR value correction process can include the corrected sequence set R2 and the non-corrected sequence set R3.

[0201] For example, for the following sequence reads that meet the correction conditions:

[0202] "TPNB500261AR:910:H5ML7BGXV:1:22202:22681:4089 163 1

[0203] 11856256 60 116M34S=11856291 164

[0204] CCCTCACCTGGATGGGAAAGATCCCGGGGACGATGGGGCAAGTGATGCCCATGTCGGTGCATGCCTTCACAAAGCGGAAGAATGTGTCAGCCTCAAAGAAAAGCTGCGTGATGATGGGCAAGTGATGCCCATGTCGGTGCATGCCTTCAC

[0205] ? >? CCAB@BCA@? CAAABBC@@B? @9AAAAC8A? DBABBACD? D@? DB@@B? D? C9B>DBA@D@@CBCBCBCCDB9BACDAC@D? D? BBB@>CBACCD@CBCDBCDA7? DA@D@@DBA=@CD? 4A14A? @BAB2C)A? 4BB1B=05CC30 NM:i:0MD:Z:116MC:Z:22S129M AS:i:116XS:i:0

[0206] RG:Z:sample SA:Z:1,11856291,+,115S35M,60,0;".

[0207] In the above sequence read segment, the FLAG value of the second column is 163, which does not contain 0x100, indicating that the alignment information corresponding to the 1st to 17th columns of the main row of the sequence is the representative alignment information; and the last column "SA:Z:1,11856291,+,115S35M,60,0" is the substitute alignment information.

[0208] In this embodiment, by classification and correction based on the classification results, the ability to identify real variations in the Softclip region is improved, and false positives caused by alignment errors are reduced. The corrected CIGAR value and position can more accurately reflect the actual variation of the sequence.

[0209] Furthermore, the linear alignment result includes a representative alignment and a substitute alignment; wherein, the linear alignment result includes a representative alignment and a substitute alignment; wherein, the representative alignment is an alignment that does not contain 0×100 in the FLAG value of the linear alignment result, and is the final alignment of the linear alignment result; the substitute alignment is an alignment that contains 0×100 in the FLAG value of the linear alignment result, and is the alignment corresponding to the SA tag; it should be noted that in high-throughput sequencing data analysis, especially when processing sequence alignment files in BAM format, FLAG values ​​are used to indicate the characteristics and status of the alignment. Each alignment has a corresponding FLAG value, which is a decimal number containing multiple flag bits used to describe the different attributes of the alignment. In the representative alignment, the main and optimal alignment result does not contain the label "not primary alignment", i.e. 0×100, in its FLAG value; while the substitute alignment refers to the alignment results other than the main alignment in the alignment results, and the FLAG values ​​of these alignments contain the label "not primary alignment", i.e. 0×100.

[0210] As mentioned above, in high-throughput sequencing data analysis, a read may be aligned to multiple positions of the reference genome, forming a mosaic alignment. These alignment results include representative alignments (RL) and alternative alignments (AL). The representative alignment is the optimal alignment selected by the alignment software, while the alternative alignment is an additional alignment position provided by the SA tag.

[0211] Specifically, they may include:

[0212] (1) RL is the representative alignment, i.e. the final alignment in the results;

[0213] (2) AL is the substitute alignment, i.e. the alignment corresponding to the SA tag.

[0214] Thus, representative alignment and substitute alignment results are obtained in the linear alignment results, which are used for subsequent CIGAR value correction.

[0215] By distinguishing between representative and alternative alignments, chimeric alignment sequences can be more accurately identified and processed, improving the accuracy of variant detection.

[0216] Reference Figure 5 The step S152, based on the classification result and according to the linear alignment result of the chimeric alignment sequence, corrects the CIGAR value and position of the representative alignment in the chimeric alignment sequence, including:

[0217] In this step, the InDel type is determined and reconstructed based on the chromosome number, position, positive and negative strand information of the linear alignment results (RL and AL) in the chimeric alignment, and the CIGAR value corresponding to RL in the chimeric alignment is modified.

[0218] Step S1521, calculating the theoretical length of the target segment;

[0219] As mentioned above, in the chimeric alignment sequence, the target fragment length is calculated based on the position information of the representative alignment and the substitute alignment. This step involves comparing the base differences between the two alignment positions. The target fragment length is thus obtained, which is an important parameter for subsequent confirmation of the type of gene variation. Accurate calculation of the target fragment length helps to identify true insertion or deletion events and reduce false positives.

[0220] In this step, the theoretical length of the target fragment can be calculated based on the premise that RL and AL are on the same chromosome and the chain type is the same positive / negative. Specifically, the theoretical length of the target fragment can be calculated by comparing the CIGAR strings of the representative alignment and the substitute alignment.

[0221] Step S1522, confirming the gene variation type and variation length of the chimeric alignment sequence according to the theoretical length of the target fragment and the actual length of the chimeric alignment sequence; wherein the variation length includes the deletion length and the insertion length;

[0222] As described above, the gene variation type (insertion or deletion) of the chimeric alignment sequence is determined based on the target fragment length and the total length of the read segment, thereby clarifying the type of gene variation and providing a basis for subsequent CIGAR value correction. In this step, accurate confirmation of the gene variation type helps to improve the specificity and sensitivity of variation detection. Specifically, the variation type can be determined by analyzing the relationship between the target fragment length and the read segment length.

[0223] Step S1523, based on the gene variation type and variation length, as well as the read length of the chimeric alignment sequence, the theoretical length of the target fragment, and the length and position of the Match fragment corresponding to the representative alignment and the alternative alignment, modify the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence.

[0224] As described above, after confirming the gene variation type and variation length of the chimeric alignment sequence, the CIGAR value and position of the representative alignment can be corrected to generate a corrected sequence set R2.

[0225] As described above, based on the gene variation type and the alignment result, the CIGAR value and position of the representative alignment in the chimeric alignment sequence are corrected. This step group involves adjusting the insertion, deletion or matching part in the CIGAR string to obtain the corrected CIGAR values, which more accurately reflect the real gene variation.

[0226] In this step, correcting the CIGAR value helps improve the accuracy of subsequent variant detection and analysis.

[0227] Furthermore, the calculation formula of the target fragment length is:

[0228] L 目标 =SM 终止 -MS 起始 ;

[0229] Among them, L 目标 represents the target fragment length; SM 终止 Represents the end position of the SM class; MS 起始 represents the starting position of the MS class;

[0230] The above formula involves calculating the target fragment length between the SM class (Softclip-Match) and the MS class (Match-Softclip) in the chimeric alignment sequence. Specifically, the target fragment length is obtained by subtracting the starting position of the MS class from the ending position of the SM class. This length is used to subsequently determine the type of genetic variation. By calculating the target fragment length, the type of genetic variation represented by the chimeric alignment sequence can be more accurately determined, whether it is an insertion or a deletion.

[0231] The step S1522, confirming the gene variation type and variation length of the chimeric comparison sequence according to the theoretical length of the target fragment and the actual length of the chimeric comparison sequence, includes:

[0232] Step S15221, if the theoretical length of the target fragment is greater than the actual length of the chimeric comparison sequence, the gene variation type of the chimeric comparison sequence is determined to be deletion; the calculation formula for the deletion length of the chimeric comparison sequence with the gene variation type being deletion is: L 缺失 =L 目标 -L 序列 Among them, L 缺失 The deletion length of the chimeric alignment sequence representing the gene variation type as deletion; L 目标represents the theoretical length of the target fragment; L 序列 is the actual length of the chimeric aligned sequence;

[0233] Step S15222: If the theoretical length of the target fragment is less than the actual length of the chimeric alignment sequence, the gene variation type of the chimeric alignment sequence is determined to be insertion. The calculation formula for the insertion length of the chimeric alignment sequence with the gene variation type being insertion is: L 插入 =L 序列 -L 目标 Among them, L 插入 The insertion length of the chimeric alignment sequence representing the genetic variation type as insertion.

[0234] In this step, it is necessary to compare the theoretical length of the target fragment with the actual length of the chimeric comparison sequence to determine the type of variation. If the theoretical length of the target fragment is greater than the actual length of the chimeric comparison sequence, it is determined to be a deletion; if the theoretical length of the target fragment is less than the actual length of the chimeric comparison sequence, it is determined to be an insertion. The type of genetic variation represented by the chimeric comparison sequence is clarified, providing a basis for the subsequent correction of the CIGAR value. Then, according to the type of genetic variation of the deletion or insertion, the corresponding deletion length and insertion length are determined. This determination method improves the accuracy of variation type identification and helps to more accurately correct the CIGAR value in subsequent steps.

[0235] Further, the step S1523 modifies the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence based on the gene variation type and variation length, the read length of the chimeric alignment sequence, the theoretical length of the target fragment, and the length and position of the Match fragment corresponding to the representative alignment and the substitute alignment, including:

[0236] Step S15231, if the classification result of the representative alignment in the chimeric alignment sequence is the MS class, the starting position is kept unchanged; if the variation type is deletion, the CIGAR value is modified from L1ML2S to L1ML 缺失 D(L 序列 -L1)M; if the mutation type is insertion, the CIGAR value is changed from L1ML2S to L1ML 插入 I(L 目标 -L1)M;

[0237] Step S15232, if the classification result of the representative alignment in the chimeric alignment sequence is the SM class, the starting position is replaced with the position of the alternative alignment MS type; if the variation type is deletion, the CIGAR value is modified from L2SL1M to L3ML 缺失 D(L 序列-L3)M; if the mutation type is insertion, the CIGAR value is changed from L2SL1M to L3ML 插入 I(L 目标 -L3)M;

[0238] Wherein, L1 represents the length of the successfully mapped segment; L2 represents the length of the Softclip segment; L3 represents the length of the successfully mapped segment of the substitute comparison; L 缺失 is the missing length; L 插入 is the insertion length; L 序列 represents the actual length of the chimeric alignment sequence; L 目标 Represents the theoretical length of the target fragment.

[0239] As mentioned above, step S1523 involves correcting the CIGAR value of the chimeric alignment sequence, which is performed after determining the gene variation type (deletion or insertion) and the variation length. Specifically, this step includes two sub-steps, which process the chimeric alignment sequences of the MS class and the SM class respectively.

[0240] In the above step S15231, the representative alignment sequences in the MS-type chimeric alignment sequences are processed. The MS-type chimeric alignment sequences refer to those sequences with successful mapping fragments on the left and Softclip fragments on the right. In this step, if the variation type is deletion, the original CIGAR value (expressed as L1ML2S, where L1 is the length of the successfully mapped fragment and L2 is the length of the Softclip fragment) needs to be modified to reflect the actual deletion length. If the variation type is insertion, the CIGAR value also needs to be modified accordingly to reflect the actual insertion length.

[0241] For missing, CIGAR value is modified to L1ML 缺失 D(L 序列 -L1)M, where L 缺失 is the length of the deletion, calculated as the length of the entire read segment L 序列 Subtract L1.

[0242] For insert, CIGAR value is changed to L1ML 插入 I(L 目标 -L1)M, where L 插入 is the length of insertion, calculated as the theoretical length L 目标 Subtract L1.

[0243] The above step S15232 processes the representative alignment sequences in the SM chimeric alignment sequences. The SM chimeric alignment sequences refer to those sequences with Softclip fragments on the left and successfully mapped fragments on the right. In this step, the starting position needs to be replaced with the position of the MS type in the substitute alignment. Then, the CIGAR value is modified accordingly according to the variation type (deletion or insertion).

[0244] For missing, CIGAR value is modified to L3ML 缺失 D(L 序列 -L3)M, where L3 is the length of the successfully mapped fragment of the alternative alignment.

[0245] For insert, CIGAR value is changed to L3ML 插入 I(L 目标 -L3)M.

[0246] Through the above method in this embodiment, it can be ensured that the CIGAR value and the alignment position accurately reflect the actual variation of the sequence, thereby improving the accuracy and reliability of subsequent variation detection.

[0247] Embodiment 3:

[0248] Reference Figure 6 Embodiment 3 of the present invention provides an insertion and deletion mutation data analysis method. In step S200, according to the CIGAR value of each sequence in the intersection sequence set, a target fragment based on the breakpoint extension of each sequence in the intersection sequence set is generated, including:

[0249] Step S210, in the intersection sequence set, selecting sequences whose gene variation type is deletion or insertion and whose deletion or insertion length is not less than 3 as screening sequences;

[0250] In the intersection sequence set, a screening sequence is obtained by screening. The specific screening conditions include:

[0251] Condition (1) The type of gene variation is deletion (D) or insertion (I).

[0252] Condition (2) The length of the missing or inserted sequence is not less than 3.

[0253] It should be noted that conditions (1) and (2) need to be met at the same time, so that the screening sequences obtained are sequences that may contain real insertion or deletion mutations, which will be used for subsequent target fragment generation. Through screening, those mutations that are too short can be excluded. These short mutations may be errors in the sequencing or alignment process, thereby improving the accuracy of subsequent analysis.

[0254] The intersection sequence set mentioned above may be the intersection sequence set obtained after step S100, or may be the intersection sequence set obtained after performing CIGAR value correction processing based on the intersection sequence set in embodiment 2.

[0255] Step S220, extending the insertion or deletion breakpoint position of the screening sequence to the left and right by a second preset length to serve as the target fragment.

[0256] As described above, at the insertion or deletion breakpoint position of the screening sequence, the second preset length is extended to the left and right to generate the target fragment. This step is to capture the complete variant sequence and its surrounding context information so as to perform subsequent comparison and analysis more accurately.

[0257] It should be noted that the second preset length extending to the left and right should be as short as possible to ensure that the Kmer fragment that most accurately represents the target deletion or insertion sequence is obtained to the greatest extent.

[0258] Extending the target fragment helps capture the complete variant sequence, including the sequence near the variant, which is important for accurately identifying and characterizing the variant.

[0259] Furthermore, in step S200, each target fragment in the first Kmer data set is subjected to preliminary screening and elimination, including:

[0260] Step S211, in the first Kmer data set, perform preliminary screening and elimination according to at least one of the following elimination rules:

[0261] A. Eliminate the target fragments whose mapping quality of the sequence where the target fragment is located is less than a preset mapping quality threshold;

[0262] B. Based on the base composition, remove the target fragments with continuous repeated bases of not less than 12 bp;

[0263] C. Determine the target fragment to be removed according to the comparison information of the sequence where the target fragment is located.

[0264] In the above, step S211 is performed to perform a preliminary screening and elimination process on the target fragments in the first Kmer data set K1 in the insertion / deletion mutation data analysis method.

[0265] In the elimination rule A, the target fragments with low mapping quality are eliminated. In the first Kmer dataset K1, the target fragments whose mapping quality of the sequence is lower than the preset threshold are removed. The mapping quality refers to the reliability of the sequence alignment with the reference genome. Thus, a higher quality Kmer dataset is obtained, in which the sequence alignment contained is more reliable.

[0266] The preset mapping quality threshold is a preset mapping quality evaluation threshold, for example, it is set to 50 by default in this embodiment.

[0267] Eliminating rule A improves the accuracy of subsequent variant detection and avoids false positive results caused by low-quality alignments.

[0268] In the elimination rule B, in the first Kmer data set K1, those target fragments with the number of consecutive repeated bases greater than or equal to 12 bp are removed. Such duplication may be caused by sequencing or alignment errors. Thus, a Kmer data set with possible erroneous duplications removed is obtained, which improves the reliability of the data. It reduces misjudgments caused by repeated sequences and improves the accuracy of variant detection.

[0269] In the elimination rule C, in the first Kmer dataset K1, based on the alignment information of the sequence where the target fragment is located, the fragments that do not meet the specific conditions, such as abnormal alignment positions or fragments with low alignment scores, are determined and eliminated. Thus, a more accurate Kmer dataset is obtained, which contains more accurate sequence alignment information. The accuracy of variant detection is further improved, and possible false positive results are eliminated.

[0270] As mentioned above, the elimination rules ABC can be selected any one of them as the rule in this embodiment, or three elimination rules can be selected. In terms of execution logic, they can be executed simultaneously or in any order.

[0271] For further reference, Figure 7 In step S211, determining the target fragment to be removed according to the comparison information of the sequence where the target fragment is located includes:

[0272] The following steps of this embodiment illustrate how to determine and eliminate target fragments that do not meet the conditions based on the alignment information of the sequence where the target fragment is located in the insertion and deletion mutation data analysis method.

[0273] Step S2111, obtaining the specific position of the sequence breakpoint where the target fragment is located in the reference genome;

[0274] First, it is necessary to determine the exact position of the insertion or deletion breakpoint in the target fragment in the reference genome. This is obtained by analyzing the alignment results, where the breakpoint position is where the variation occurs. The exact position of the breakpoint is obtained, providing a basis for the subsequent calculation of the start and end positions.

[0275] Step S2112, according to the gene variation type and the specific position, respectively calculate the start and end positions of the breakpoint; the calculation method includes:

[0276] If the gene variation type is insertion, the start and end positions of the breakpoint include:

[0277] Starting position: P 插入起始 =P 具体 -L 插入长度 ;

[0278] End position: P 插入终止 =P 具体 +L 插入长度 ;

[0279] The expression of the start and end positions can be: [indel_pos-len(ins):indel_pos+len(ins)].

[0280] If the gene variation type is deletion, the start and end positions of the breakpoints include:

[0281] Starting position: P 缺失起始 =P 具体 ;

[0282] End position: P 缺失终止 =P 具体 +L 缺失长度 ;

[0283] The expression of the start and end positions is [indel_pos:indel_pos+len(del)].

[0284] Among them, P 插入起始 Represents the starting position of the breakpoint where the gene variation type is insertion; P 插入终止 represents the end position of the breakpoint where the gene variation type is insertion; P 缺失起始 Represents the starting position of the breakpoint where the gene variation type is deletion; P 缺失终止 Represents the end position of the breakpoint where the gene variation type is deletion; P 具体 Represents the specific location; L 插入长度 represents the length of the inserted sequence; L 缺失长度 represents the length of the missing sequence;

[0285] In the above expression, indel_pos represents the specific position; len(ins) and len(del) represent the length of the inserted sequence and the length of the deleted sequence, respectively.

[0286] As described above, the start and end positions of the breakpoints are calculated according to the gene variation type (insertion or deletion) and the specific position of the breakpoints, thereby obtaining the start and end positions of the breakpoints, which are used to define the range of the target fragment.

[0287] In this step, by precisely defining the range of the target fragment, potential insertion or deletion events can be more accurately analyzed and evaluated. The specific calculation method can be:

[0288] (1) If the mutation type is insertion, the start position is the breakpoint position minus the length of the inserted sequence, and the end position is the breakpoint position plus the length of the inserted sequence.

[0289] (2) If the variant type is deletion, the starting position is the breakpoint position, and the ending position is the breakpoint position plus the length of the deletion sequence.

[0290] Step S2113: if the interval where the start and end positions of the breakpoint are located does not intersect with the preset target capture interval, the target segment corresponding to the breakpoint is removed.

[0291] As mentioned above, if the interval where the start and end positions of the breakpoints are located does not intersect with the preset target capture interval, it means that the target fragment is not located in the genomic region of interest, so it is removed. The target fragments that do not meet the conditions are removed, and the target fragments that intersect with the preset target capture interval are retained.

[0292] The above steps help to exclude irrelevant or unrepresentative target fragments, thereby improving the accuracy and efficiency of data analysis.

[0293] Embodiment 4:

[0294] Reference Figure 8 Embodiment 4 of the present invention provides an insertion and deletion mutation data analysis method. In step S200, after performing preliminary screening and elimination processing on each of the target fragments in the first Kmer data set to obtain a second Kmer data set, the method further includes:

[0295] Step S500, based on the second Kmer data set, according to the length of the inserted or deleted sequence in each of the target fragments, extending to the breakpoints at both ends of the sequence to generate Extend fragments, and constructing a first Extend fragment data set;

[0296] In this embodiment, after step S200, for the second Kmer dataset K2, a method for further screening target fragments after the initial screening of the Kmer dataset is performed to improve the accuracy and efficiency of insertion / deletion mutation data analysis.

[0297] It should be noted that the unscreened Kmer dataset may contain too short sequences, which may lead to inaccurate alignment results and fail to ensure accurate variation information. Therefore, sequence extension is performed based on the second Kmer dataset K2 to generate Extend fragments, the purpose of which is to ensure the accuracy of downstream local realignment.

[0298] To solve the above problem, in this embodiment, based on the second Kmer dataset K2, according to the length of the inserted or deleted sequence in each target fragment, the fragments are extended to the breakpoints at both ends of the sequence to generate Extend fragments. This step is to capture a more complete variant sequence, including the contextual information around the variant. Thus, a set of Extend fragments is obtained, which covers the original target fragment and the sequence around it.

[0299] The Extend fragment provides more comprehensive information about the variant region, which helps to improve the sensitivity and accuracy of variant detection.

[0300] As described above, the generated Extend fragments and their related information (such as position, length, type, etc.) are stored in a data set to form a first Extend fragment data set E1. Thus, a data set containing detailed variant region information is obtained, which provides a basis for subsequent variant analysis. Organizing and storing detailed variant region information facilitates subsequent processing and analysis.

[0301] Step S600: In the first Extend segment data set, one or two representative Extend segments are obtained for each bed interval as target segments after screening of the interval, and the second Extend segment data set including the target segments is generated.

[0302] As described above, in the first Extend segment dataset E1, the Extend segments in each bed interval are counted, and the most representative 1 or 2 are selected as the representative segments of the interval to form the second Extend segment dataset E2. Thus, a streamlined representative segment dataset is obtained, in which each bed interval contains the representative segment that best characterizes the target InDel.

[0303] By selecting representative fragments, the computational effort of subsequent analysis is reduced and the analysis efficiency is improved.

[0304] In this embodiment, based on the initial screening in step S200, the method in steps S600-S700 is used for re-screening. By screening representative fragments, unnecessary calculations are reduced and the analysis speed is improved; representative fragments are more likely to reflect real mutations, thereby improving the accuracy of mutation detection. The streamlined data set reduces the demand for computing resources and makes the analysis process more efficient.

[0305] In summary, in this embodiment, an effective strategy is provided for further screening and optimizing target fragments after the initial screening of the Kmer dataset to improve the overall performance of insertion / deletion mutation data analysis.

[0306] Furthermore, in step S500, the method for calculating the extension length extending to the breakpoints at both ends of the sequence includes:

[0307] (1) If the gene mutation type corresponding to the target fragment is insertion, the length of the target fragment extending to the left and right is:

[0308] L 延伸 =C 可调 +L I ×2;

[0309] (2) If the gene variation type corresponding to the target segment is deletion, the length of the target segment extending to the left and right is:

[0310] L 延伸 =C 可调 ;

[0311] Among them, L 延伸 Represents the length of the target segment extending to the left and right, and does not include the Softclip area; C 可调 Represents preset adjustable parameters; L I Represents the sum of the lengths of the inserted fragments in the CIGAR value.

[0312] Above, C 可调 Represents a preset adjustable parameter, which can be represented by extend_seq_size. The default value used in this embodiment is 15.

[0313] As described above, for each target fragment, the length extending to the breakpoints at both ends of the sequence is calculated according to the corresponding gene variation type (insertion or deletion), thereby obtaining the extension length of each target fragment, which is used to generate the Extend fragment.

[0314] By accurately calculating the extension length, we can ensure that the Extend fragment can fully cover the variant region, improve the accuracy of variant detection, and thus ensure the accuracy of downstream local realignment.

[0315] As described above, the calculated extension length is used to extend to the breakpoints at both ends of the target fragment to generate Extend fragments, thereby obtaining a set of Extend fragments, which will be used to construct the first representative fragment data set. By generating Extend fragments, the variant region can be analyzed more comprehensively, improving the sensitivity and accuracy of variant detection.

[0316] Then, the generated Extend fragments and their related information (such as position, length, type, etc.) are stored in a data set to form a first representative fragment data set E1, thereby obtaining a data set containing detailed variant region information, providing a basis for subsequent variant analysis.

[0317] For further reference, Fig. 9, the step S600, in which one or two representative Extend segments are obtained from each bed interval in the first Extend segment data set as target segments after screening of the interval, and the second Extend segment data set containing the target segments is generated, includes:

[0318] Step S610, classify the first Extend fragment data set according to the length of the insertion or deletion sequence; wherein the sequences with the length of the insertion or deletion sequence not less than 6 are used as the E1a data set; and the sequences with the length of the insertion or deletion sequence less than 6 are used as the E1b data set.

[0319] Based on the first Extend segment data set E1, two data subsets E1a and E1b are formed according to the length of the insertion or deletion sequence in the target segment. E1a contains data with an insertion or deletion sequence length of not less than 6, while E1b contains data with a length less than 6.

[0320] Through classification, two classified data subsets are obtained, which facilitates the specialized analysis of variants of different lengths. Based on the above classification steps, different types of variants can be processed more specifically, improving the accuracy and efficiency of variant detection.

[0321] Step S620, sorting the Extend fragments in the classified E1a data set and the E1b data set according to the number of Extend fragments, the number of target fragments, the quality of sequence alignment, and the length of the inserted or deleted sequence, and selecting the first-ranked Extend fragment in the E1a data set and the E1b data set as the best representative Extend fragment;

[0322] As described above, the Extend fragments in the classified E1a and E1b data sets were sorted according to the number of Extend fragments, the number of target fragments, the quality of sequence alignment, and the length of the inserted or deleted sequences. According to the sorting results, the first-ranked Extend fragment was selected as the best representative Extend fragment in the E1a and E1b data sets, respectively.

[0323] In this step, one or two optimal representative Extend fragments are selected in each data set as representative fragments of the interval. Selecting representative fragments helps to focus resources on analyzing sequences that are most likely to contain real variants, thereby improving the efficiency and accuracy of the analysis.

[0324] Step S630, comparing the similarity and the difference in the number of Extend segments between the longest common subsequences of the optimal representative Extend segments of the E1a data set and the E1b data set, and obtaining a comparison result; including:

[0325] A. If the similarity is greater than 80%, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set;

[0326] B. If the similarity is not greater than 80%, and the number of Extend segments is less than 3, and the difference in the number of Extend segments is more than 30 times, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set;

[0327] C. If both A and B do not match, the comparison result is: both the E1a data set and the E1b data set are retained to generate the second Extend segment data set.

[0328] As described above, the similarities between the longest common subsequences of the best representative Extend segments selected from the E1a and E1b data sets were compared, as well as the differences in the number of Extend segments.

[0329] In this step, it is determined whether to retain one or two representative fragments based on the similarity and quantity difference. By comparing the similarity and quantity difference, the most representative fragments can be selected more accurately, reducing false positive and false negative results.

[0330] Then, based on the above comparison results, the second Extend fragment dataset E2 is generated, which contains 1 or 2 representative fragments in each bed interval. A streamlined and representative variation dataset is obtained, providing high-quality input for subsequent analysis. Through screening and comparison, it is ensured that the dataset contains the most representative and reliable variation information.

[0331] Embodiment 5:

[0332] Reference Fig.10 Embodiment 5 of the present invention provides an insertion and deletion mutation data analysis method, wherein step S300, based on the preset target capture interval, re-aligns the target fragment to obtain an alignment result, including:

[0333] Step S310, extending the preset target capture interval to the left and right by a third preset length to obtain a second bed interval;

[0334] As mentioned above, the bed interval of the preset targeted capture interval is the targeted capture bed interval B0 in this embodiment. In this embodiment, the second bed interval can be defined as B2.

[0335] As mentioned above, the preset targeted capture interval (bed interval B0) is extended to create a larger analysis interval (second bed interval). This step takes into account the possible existence of larger InDel events and ensures that these events can be fully included in the analysis interval, thereby obtaining an extended bed interval, which provides a wider analysis range for local realignment.

[0336] In this step, by extending the bed interval, the possibility of detecting InDel events near the boundary is increased.

[0337] In the above, the bed interval of the preset target capture interval is extended to obtain a second bed interval;

[0338] As described above, based on the defined preset targeted capture interval B0, in order to more comprehensively capture possible insertion or deletion mutations, these regions need to be further extended. By extending the third preset length to the left and right, a second bed interval B2 is obtained. This interval is larger than the original targeted capture bed interval and can cover a wider range of genomic regions.

[0339] The third preset length mentioned above can be represented by extend_bed_size, which is an adjustable length parameter and is 75 by default in this embodiment.

[0340] By extending the process, the possibility of capturing mutations near the boundary can be increased, thereby improving the sensitivity and accuracy of mutation detection.

[0341] Step S320, re-aligning the target fragment with the second bed interval to obtain the alignment result.

[0342] It should be noted that the target fragment in this embodiment can be the target fragment after the initial screening process contained in the second Kmer data set in step S200 in Example 1, or it can be the target fragment (Extend) in the second Extend fragment data set in Example 4.

[0343] As described above, a local alignment algorithm (such as the Smith-Waterman algorithm) is used to locally realign the target fragment with the second bed interval. This step is to accurately identify and correct insertion or deletion events in the sequence.

[0344] First, you can select the alignment algorithm and choose a suitable local alignment algorithm, such as the Smith-Waterman algorithm, which is a commonly used local sequence alignment algorithm that is particularly suitable for accurate alignment of short sequences.

[0345] During the local realignment process, specific scoring rules are used.

[0346] Finally, the alignment results are analyzed to analyze the results of the local realignment and identify the best alignment scheme, which is usually the scheme with the highest cumulative score based on the above scoring rules.

[0347] Local realignment allows for more accurate identification of insertion or deletion events in a sequence, especially those that may have been missed or misidentified in the original alignment. Optimized alignment algorithms and scoring rules help reduce false positive results, which are incorrect identifications of non-real variants. Local realignment can improve the detection of complex variants, such as mosaic insertions and deletions, which is critical to understanding the diversity and complexity of the genome. By focusing on local realignment of target fragments, computing resources can be used more efficiently, improving the efficiency of analysis.

[0348] In summary, local realignment plays a vital role in the insertion and deletion mutation data analysis method. It improves the accuracy and efficiency of variation detection through precise alignment algorithms and optimized scoring rules.

[0349] Furthermore, the scoring rule expression in the re-alignment algorithm is:

[0350] [match:1,mismatch:1,gap open:2,gap extend:0];

[0351] Among them, Match: 1 means that if the two aligned bases are the same, 1 point is obtained; Mismatch: -1 means that if the two aligned bases are different, 1 point is deducted; Gap Open: -2 means that if a gap is introduced, 2 points are deducted; Gap Extend: 0 means that if a gap already exists, no point is deducted.

[0352] The above steps are used to define the scoring rules of the local realignment algorithm in the insertion and deletion mutation data analysis method, and these rules are used to evaluate and optimize the quality of the alignment results.

[0353] As mentioned above, Match score means that during the local realignment process, if the two aligned bases are the same, that is, they match in the sequence alignment, then this situation is considered a Match and gets a score of 1. The Match score is used to encourage the matching of bases in the alignment, thereby improving the accuracy of the alignment.

[0354] As mentioned above, the mismatch score means that if two bases are different in the alignment, i.e. they do not match, this situation is considered a mismatch and 1 point is deducted. The mismatch deduction is used to penalize mismatched bases and reduce alignment errors.

[0355] As mentioned above, Gap Open scoring means that if a new gap is introduced in the alignment (whether in the reference sequence or the read sequence), this situation is considered a Gap Open and a score of 2 is deducted. Gap Open deduction is used to avoid unnecessary gaps, because gaps usually indicate insertions or deletions in the sequence.

[0356] The above, Gap Extend scoring, means that if there is already a gap and this gap is extended further, this situation is considered a Gap Extend and no points are deducted. The rule of no points deduction for Gap Extend encourages extending existing gaps rather than introducing new gaps.

[0357] The way these scoring rules are implemented usually involves the use of dynamic programming algorithms, such as the Smith-Waterman algorithm, which is an algorithm for local sequence alignment that finds the optimal alignment by computing and comparing all possible alignments.

[0358] Through these scoring rules, the algorithm can evaluate different alignment schemes and select the scheme with the highest score as the final alignment result. This method helps to improve the accuracy and reliability of indel mutation detection.

[0359] Embodiment 6:

[0360] In order to make the purpose, technical scheme and advantages of the present invention more clear and to better illustrate the insertion and deletion mutation data analysis method provided by the present invention, in this example, the specific scheme is described in detail in combination with the specific implementation method and the accompanying drawings. It should be noted that the specific genomic positions and sequence contents listed in the description are only used to illustrate the logical processing. The data content that actually conforms to the rules can be applied to this analysis method without special indication.

[0361] This embodiment provides an insertion and deletion mutation data analysis method, comprising:

[0362] Step S1, initial screening and preprocessing of BAM sequence set.

[0363] (1) First, according to the chromosome number and start and end positions of the input bed interval (preset target capture interval) B0, a bed interval [7:116411861-116411910] on chromosome 7 of the human genome database is taken as an example (the human genome database has multiple versions, and the b37 version is used in this embodiment: [ftp: / / ftp.broadinstitute.org / pub / seq / references / Homo_sapiens_assembly19.fasta, MD5sum:886ba1559393f75872c1cf459eb57f2d]), the reads (intersection sequence R1) that have an intersection with the bed interval (first bed interval B1) [7:116411711-116412060] after extending 150 bp to the left and right are searched in the BAM file (sequence alignment file).

[0364] It should be noted that the choice of extending the input bed interval by 150 bp to the left and right is to take into account all InDels with a length of 3-150 bp so that the data set obtained from the initial screening is not missed.

[0365] (2) Then, the sequences marked as duplicates are filtered out, and the chimeric alignment sequences that meet the correction conditions are subjected to CIGAR value correction processing to obtain an intersection sequence set including a corrected sequence set and a non-corrected sequence set. For example, in the sample, there is a chimeric alignment sequence with the original chromosome number, alignment start position, and CIGAR value of {chr:7,pos:117411706,cigar:42S108M}, which is corrected to {chr:7,pos:117411727,cigar:45M63I42M} after correction processing. The processing process is as follows: Fig.11 shown.

[0366] It should be noted that the reason why the CIGAR value needs to be corrected here is that the present invention identifies InDels by CIGAR values, and the omission of InDels in chimeric alignment sequences can be avoided by correcting the CIGAR value.

[0367] Step S2, generation and initial screening of Kmer data set.

[0368] (1) Based on the intersection sequence set processed in step S1, the sequences with CIGAR values ​​containing D or I and InDel length greater than or equal to 3 are selected, and Kmer fragments are generated by extending 17 bp to the left and right along each breakpoint of each sequence as the Kmer data set for characterizing each breakpoint.

[0369] by Fig.12As an example, according to the CIGAR value 127M4D24M of the sequence, a Kmer fragment "CTTTCTTTCTCTCTGGTGATCTGGGCAGTGAATT" with a length of 34 bp is obtained from the position between the 127th and 128th bases of the sequence to the left and right, respectively, to characterize the 4 bp deletion at the breakpoint. Fig.13 For specific demonstration.

[0370] It should be noted that the Kmer extension length should be as short as possible to obtain the specific Kmer fragment that most accurately represents the target InDel and reduce the possibility of introducing other mutations. The default extension of 17 bp is the best parameter obtained through testing. It can be modified through testing in actual use. The recommended parameter value range is [15-18].

[0371] (2) Then eliminate Kmers that do not meet the requirements according to the following conditions:

[0372] According to the mapping quality, Kmers whose mapping quality of the sequence in which the Kmer is located is less than 50 are eliminated to obtain Kmers with higher credibility;

[0373] According to the base composition of Kmer, Kmers with the number of consecutive repeated bases greater than or equal to 12 bp were eliminated, for example, “CTTTCTTTCTAAAAAAAAAAAAGGGCAGTGAATT” would be eliminated.

[0374] It should be noted that this is also to ensure that the Kmer obtained is credible.

[0375] According to the alignment information of the Kmer sequence, obtain the specific position of the breakpoint on the reference genome, indel_pos;

[0376] A. If the mutation type is insertion, the start and end positions are [indel_pos-len(ins):indel_pos+len(ins)];

[0377] B. If the mutation type is deletion, the start and end positions are [indel_pos:indel_pos+len(del)];

[0378] C. If the interval where the start and end positions are located has no intersection with the input bed interval, the Kmer will be eliminated.

[0379] by Fig.12Taking the sequence as an example, the start and end positions of the deletion in the sequence are [7:116411898-116411902], which is within the input bed interval [7:116411861-116411910], so the Kmer where the deletion is located will be retained.

[0380] It should be noted that the processing here is to accurately obtain the target InDel within the specified input interval, eliminate the influence of other irrelevant InDels, and improve the downstream analysis speed.

[0381] Step S3, generation of Extend data set and screening of interval representative fragments.

[0382] (1) Based on the Kmer dataset generated in step S2, the Kmer is dynamically extended to both ends of the sequence again according to the length of the InDel to generate an Extend dataset (the first Extend fragment dataset E1).

[0383] Specifically, if it is an insertion type mutation, the extension length of the left and right ends of the Kmer is extend_seq_size+length*2;

[0384] If it is a deletion type mutation, the extension length of the left and right ends of the Kmer is extend_seq_size;

[0385] During this period, it is necessary to ensure that the extended sequence does not contain the Softclip area. The extend_seq_size is an adjustable parameter, the default is 15, and the length is the corresponding length of InDel.

[0386] by Fig.12 Take the sequence as an example. Because it is a deletion mutation, it will be further extended 15bp along the left and right ends of the sequence based on the Kmer fragment "CTTTCTTTCTCTCTGGTGATCTGGGCAGTGAATT" to generate the Extend fragment "CGTCTTTAACAAGCTCTTTCTTTCTCTCTGGTGATCTGGGCAGTGAAT TAGTTCGC". If the distance to the beginning or end of the sequence is less than 15, it will be extended to the beginning or end of the sequence. Fig.13 For specific demonstration.

[0387] It should be noted that the downstream local realignment result is closely related to the length of the input sequence, especially the insertion type of repeated mutations, which require at least a sequence greater than 2 times the length of the repeated fragment to be identified by local realignment. Therefore, the Kmer sequence needs to be extended here to ensure the accuracy of the downstream local realignment. After testing, the effect is best when the extend_seq_size parameter length is 15bp. It can also be adjusted according to actual conditions. The recommended value range is [15-18].

[0388] (2) Based on the statistical information in the Extend data set, obtain one or two most representative Extends as representative segments of the interval.

[0389] Specifically, the Extend dataset E1 is first divided into two parts according to whether the length of InDel is greater than or equal to 6, which are respectively represented as dataset E1a and dataset E1b; then E1a and E1b are sorted according to the number of Extend fragments, the number of Kmers, the quality of sequence alignment, and the length of InDel, and the first one in E1a and E1b is selected as the best representative Extend fragment; then the similarity between the longest common subsequences of the best representative Extend fragments of E1a and E1b and the difference in the number of Extend fragments are compared:

[0390] A. If the similarity is greater than 80%, only the ones with a larger number will be retained;

[0391] B. Otherwise, if there are less than 3 fragments and the difference between the two is more than 30 times, the one with the larger number will also be retained;

[0392] C. In other cases, both are retained.

[0393] In this way, the interval representative segment data set (the second Extend segment data set E1) is obtained.

[0394] Here is an example of the result: The interval [7:116411861-116411910] represents the fragment dataset: [('CGTCTTTAACAAGCTCTTTCTTTCTCTCTGGTGATCTGGGCAGTGAA TTAGTTCGCTACGATGC',[118,60.0,172,4,'CTTTCTTTCTCTCTGGTGATCTGGGCAGTGAATT',0])].

[0395] It should be noted that the screening model here is constructed to obtain the 1 to 2 most likely real mutations within the target interval.

[0396] Step S4: local re-alignment.

[0397] (1) First, extend the input bed interval [7:116411861-116411910] to the left and right by extend_bed_size: 75 to obtain the extended bed interval [7:116411786-116411985], and obtain the reference sequence "TTGTGTGCTGTCTTATATGTAGTCCATAAAACCCATGAGTTCTGGGCA CTGGGTCAAAGTCTCCTGGGGCCCATGATAGCCGTCTTTAACAAGCTCTTTCTTTCTCTCTGTTTTAAGATCTGGGCAGTGAATTAGTTCGCTACGATGCAAGAGTACACACTCCTCATTTGGATAGGCTTGTAAGTGCCCGAAGTGTAAGC" within the extended bed interval.

[0398] (2) Then, the interval representative Extend fragment "CGTCTTTAACAAGCTCTTTCTTTCTCTCTGGTGATCTGGGCAGTGAAT TAGTTCGCTACGATGC" obtained in step S3 is locally aligned with the reference sequence (SW algorithm), wherein the SW algorithm uses the [match: 1, mismatch: 1, gap open: 2, gap extend: 0] scoring rule.

[0399] After alignment, the mutation site {chr:7,pos:116411896,ref:TTTTAA,alt:GT} was obtained.

[0400] It should be noted that the SW local realignment is performed here to obtain more accurate mutation information, and the alignment parameters are the best parameter combination after multiple tests. In this embodiment, it is not recommended to modify it.

[0401] Step S5, mutation frequency calculation.

[0402] (1) Using indelpost, the alignment results are first standardized to obtain accurate mutation site information. For example, in other embodiments, there is such a mutation site {chr:7,pos:55242474,ref:AGAAGC,alt:C}, which will become a more accurate representation method {chr:7,pos:55242470,ref:TAAGAG,alt:T} after standardized left alignment processing.

[0403] (2) Then the mutation depth statistics and mutation frequency calculations are performed on the mutation sites.

[0404] The specific implementation method is as follows:

[0405] Step S5-1, determination and processing of complex insertions and deletions: If a mutation event consists of only a single insertion or deletion, it is determined to be a simple InDels (simple_indels), and the remaining InDel types are determined to be complex InDels (complex_indels).

[0406] For complex InDels, the decompose_complex_variant method is used to decompose them into smaller units, and the largest unit is analyzed; for simple InDels, they are analyzed directly.

[0407] Step S5-2, target sequence identification: check whether the current variation is equivalent to the target variation, and set the target state of the sequence (i.e., whether it is the target sequence, whether it is worth re-aligning, etc.) according to the CIGAR string, equivalence check results, sequence quality and other information.

[0408] Step S5-3, consensus template sequence (Contig) construction: randomly select a portion of the target sequence to construct a Contig sequence.

[0409] Step S5-4, grid search: Use the Smith-Waterman algorithm for local alignment and perform a grid search to automatically select the best alignment parameters from multiple combinations corresponding to the alignment parameters gap_open_penalty and gap_extension_penalty [(3,0), (3,1), (5,1), (5,0), (4,1), (4,0)], and update the target state of the sequence according to the alignment results corresponding to the best alignment parameters.

[0410] Step S5-5, local realignment: using the Smith-Waterman algorithm and the best alignment parameters obtained in the grid search step, local realignment is performed on the sequences whose target states are marked as worthy of realignment, and the target states of the sequences are updated according to the alignment results.

[0411] Step S5-6, depth statistics and mutation frequency calculation: The target mutation depth and reference sequence depth before and after deduplication are counted according to the target state of the sequence, and the mutation frequency is calculated.

[0412] Step S6, the final variation set is obtained.

[0413] The standardized mutation information and mutation frequency were collected to obtain the final mutation set {Gene:MET,Exon:14,cHGVS:c.2888-7_2888-2delinsGT,AF:0.121,chr:7,start:116411896,end:116411901,ref:TTTTAA,alt:GT}. The result here is the reference result based on the above example.

[0414] The method of the present invention is further described below through specific investigation experiments, but it should be understood that these investigation experiments are only used for more detailed description and should not be understood as limiting the present invention in any form.

[0415] Method performance investigation:

[0416] In order to evaluate the detection efficiency of the present invention for Small InDels, in this experiment, 120 variant sites including 71 simple InDels (simple_indels) and 49 complex InDels (complex_indels) were randomly selected from the GIAB standard indels variation set as the standard variation set, and the sensitivity, precision and F1score of the analysis method of the present invention and other commonly used indel detection tools pindel and mutect2 were calculated for a comprehensive investigation.

[0417] Experimental results:

[0418] According to the comparison results in Table 2, we can find that:

[0419] (1) For simple_indels, the sensitivity of the analysis method of the present invention (0.634) is lower than that of pindel (0.986), but the precision (0.957) is higher than that of pindel (0.673). The analysis method of the present invention and pindel are higher than mutect2 (0.423, 0.231) in detection sensitivity and precision.

[0420] (2) For complex_indels, the sensitivity and accuracy of the analysis method of the present invention (0.429, 0.382) are much higher than those of pindel (0, 0) and mutect2 (0, 0).

[0421] (3) Based on the F1score corresponding to all variants (all_records) including simple_indels and complex_indels, the overall performance of the analysis method of the present invention (0.592) is better than pindel (0.522) and mutect2 (0.2).

[0422] (4) In terms of running speed, the analysis method of the present invention (27.5S) is also faster than pindel (439.3S) and mutect2 (86S).

[0423] Table 2. Comparison results of small_indel, pindel and mutect2

[0424]

[0425] It should be noted that in Table 2, the GIAB standard indels variation set selected by the three methods is "HG001_GRCh37_benchmark_indel"; "this application" represents the insertion and deletion mutation data analysis method provided by this application.

[0426] The unit of total running speed is S.

[0427] “total.truth” represents the known true results of the standard data set; “total.query” represents the results corresponding to the detection tool to be evaluated; “tp” represents the true positive result; “fp” represents the false positive result; and “fn” represents the false negative result.

[0428] “recall” means the recall rate (ie, sensitivity); the calculation formula is: recall = tp / (tp+fn).

[0429] “Precision” refers to the accuracy, and the calculation formula is: precision = tp / (tp+fp).

[0430] "F1score" is the harmonic mean of precision and recall, which strikes a balance between the two. When both precision and recall are high, the F1 score will also be high. The calculation formula is: F1score = 2*(precision*recall) / (precision+recall).

[0431] “all_records” indicates all variant sites including simple_indels and complex_indels.

[0432] The results show that the data analysis method of local variations containing long insertion and deletion fragments based on targeted high-throughput sequencing of the present invention can quickly, accurately and sensitively detect Small InDels in tissue or blood tumor specimens (including but not limited to tumor samples), thereby improving the detection efficiency of Small InDels, especially for complex variations.

[0433] In addition, reference Fig.14In an embodiment of the present application, a device for analyzing insertion and deletion mutation data is also provided, comprising:

[0434] The alignment intersection module 10 is used to obtain a sequence alignment file corresponding to the target sample, and extract the intersection sequence existing after the alignment of the sequence alignment file with the first bed interval of the preset target capture interval to obtain an intersection sequence set;

[0435] A preliminary screening module 20 is used to generate a target fragment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set to obtain a first Kmer data set containing the target fragment; perform preliminary screening and elimination processing on each of the target fragments in the first Kmer data set to obtain a second Kmer data set;

[0436] A re-alignment module 30, used to re-align the target fragment based on the preset target capture interval to obtain an alignment result;

[0437] The result generating module 40 is used to obtain the mutation site information according to the comparison result, and obtain the variation set based on the mutation site information.

[0438] In addition, an insertion-deletion mutation data analysis system is also provided in an embodiment of the present application, comprising a memory and a processor, wherein an insertion-deletion mutation data analysis program is stored in the memory, and the processor runs the insertion-deletion mutation data analysis program so that the insertion-deletion mutation data analysis system performs the insertion-deletion mutation data analysis method as described above.

[0439] In addition, a computer-readable storage medium is provided in an embodiment of the present application, on which an insertion-deletion mutation data analysis program is stored. When the insertion-deletion mutation data analysis program is executed by a processor, the insertion-deletion mutation data analysis method as described above is implemented.

[0440] In several embodiments provided in the present application, it should be understood that the disclosed devices and methods can also be implemented in other ways. The device embodiments described above are merely schematic. For example, the flowcharts and structure diagrams in the accompanying drawings show the possible architecture, functions and operations of the devices, methods and computer program products according to multiple embodiments of the present invention. In this regard, each box in the flowchart or block diagram can represent a module, a program segment or a part of a code, and the module, a program segment or a part of a code contains one or more executable instructions for implementing the specified logical function. It should also be noted that in an alternative implementation, the functions marked in the box can also occur in a different order from the order marked in the accompanying drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the structure diagram and / or the flow chart, and the combination of boxes in the structure diagram and / or the flow chart, can be implemented with a dedicated hardware-based system that performs a specified function or action, or can be implemented with a combination of dedicated hardware and computer instructions.

[0441] In addition, the functional modules or units in the various embodiments of the present invention may be integrated together to form an independent part, or each module may exist independently, or two or more modules may be integrated to form an independent part.

[0442] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art or the part of the technical solution, can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for a computer device (which can be a smart phone, a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the methods described in each embodiment of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk.

[0443] The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of changes or substitutions within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A method for analyzing insertion and deletion mutation data, characterized in that: include: Obtaining a sequence alignment file corresponding to the target sample, and extracting the intersection sequence existing after the sequence alignment file is aligned with the first bed interval of the preset target capture interval, to obtain an intersection sequence set; According to the CIGAR value of each sequence in the intersection sequence set, a target fragment based on the breakpoint extension of each sequence in the intersection sequence set is generated to obtain a first Kmer data set containing the target fragment; each target fragment in the first Kmer data set is preliminarily screened and eliminated to obtain a second Kmer data set; Based on the preset targeted capture interval, re-aligning the target fragment to obtain an alignment result; Mutation site information is obtained according to the comparison result, and a variation set is obtained based on the mutation site information.

2. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: The method for obtaining the first bed interval of the preset targeted capture interval includes: According to the preset targeted capture interval, the preset targeted capture interval is extended to the left and right by a first preset length respectively to obtain the first bed interval; wherein the first preset length represents a maximum InDel length threshold parameter that can be extended; The calculation method of the genome start position and the genome end position of the first bed interval is: The genomic starting position of the first bed interval is: P 基因组起始 =start-L 第一 ; The genomic end position of the first bed interval is: P 基因组终止 =end+L 第一 ; Among them, P 基因组起始 represents the genomic starting position of the first bed interval; P 基因组终止 represents the genomic end position of the first bed interval; start represents the genomic starting position of the preset targeted capture interval; end represents the genomic end position of the preset targeted capture interval; L 第一 Represents the first preset length.

3. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: The step of obtaining a sequence alignment file corresponding to the target sample and extracting an intersection sequence existing after aligning the sequence alignment file with the first bed interval of the preset target capture interval to obtain an intersection sequence set includes: Taking the sequence alignment file as input, searching for a sequence matching the sequence alignment file in the first bed interval as a matching sequence; Performing filtering and de-duplication processing on the matching sequence to obtain the intersection sequence; Determining whether there is a sequence in the intersection sequence set that meets the correction condition; If yes, the sequence satisfying the correction condition is used as the sequence to be corrected; and the sequences other than the sequence to be corrected in the intersection sequence set are used as the non-corrected sequence set; Performing CIGAR value correction processing on the sequence to be corrected to obtain a corrected sequence set; and taking the corrected sequence set and the non-corrected sequence set together as the corrected intersection sequence set; The correction conditions are: chimeric alignment sequences that simultaneously meet the following conditions: A. with SA tag; B. FLAG value does not contain 0×100; C. aligned to the same chromosome; D. same chain type.

4. The insertion / deletion mutation data analysis method according to claim 3, characterized in that: The CIGAR value correction process includes: Classifying the read segments of the sequence to be corrected based on the CIGAR values ​​to obtain a classification result; Based on the classification results, according to the linear alignment results of the chimeric alignment sequences, the CIGAR values ​​and positions of the representative alignments in the chimeric alignment sequences are corrected to obtain the corrected sequence set.

5. The method for analyzing insertion-deletion mutation data according to claim 4, characterized in that: The classification results include MS and SM categories; In the MS class, the left side is the successfully mapped fragment, and the right side is the Softclip fragment; In the SM class, the left side is the Softclip segment, and the right side is the successfully mapped segment.

6. The method for analyzing insertion-deletion mutation data according to claim 5, characterized in that: The linear alignment result includes a representative alignment and a substitute alignment; wherein the representative alignment is an alignment whose FLAG value of the linear alignment result does not contain 0×100 and is the final alignment of the linear alignment result; the substitute alignment is an alignment whose FLAG value of the linear alignment result contains 0×100 and is the alignment corresponding to the SA tag; The method of correcting the CIGAR value and position of the representative alignment in the chimeric alignment sequence based on the classification result according to the linear alignment result of the chimeric alignment sequence comprises: Calculate the theoretical length of the target fragment; According to the theoretical length of the target fragment and the actual length of the chimeric comparison sequence, confirm the gene variation type and variation length of the chimeric comparison sequence; wherein the variation length includes the deletion length and the insertion length; Based on the gene variation type and variation length, the read length of the chimeric alignment sequence, the theoretical length of the target fragment, the length and position of the Match fragment corresponding to the representative alignment and the substitute alignment, the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence are modified.

7. The method for analyzing insertion-deletion mutation data according to claim 6, characterized in that: The calculation formula of the theoretical length of the target fragment is: L 目标 =SM 终止 -MS 起始 ; Among them, L 目标 represents the theoretical length of the target fragment; SM 终止 Represents the end position of the comparison of the SM class; MS 起始 represents the alignment start position of the MS class; The step of confirming the gene variation type and variation length of the chimeric comparison sequence according to the theoretical length of the target fragment and the actual length of the chimeric comparison sequence comprises: If the theoretical length of the target fragment is greater than the actual length of the chimeric comparison sequence, the gene variation type of the chimeric comparison sequence is determined to be deletion; the calculation formula for the deletion length of the chimeric comparison sequence whose gene variation type is deletion is: L 缺失 =L 目标 -L 序列 Among them, L 缺失 The deletion length of the chimeric alignment sequence representing the gene variation type as deletion; L 目标 represents the theoretical length of the target fragment; L 序列 is the actual length of the chimeric alignment sequence; If the theoretical length of the target fragment is less than the actual length of the chimeric comparison sequence, the gene variation type of the chimeric comparison sequence is determined to be insertion. The calculation formula for the insertion length of the chimeric comparison sequence with the gene variation type being insertion is: L 插入 =L 序列 -L 目标 Among them, L 插入 The insertion length of the chimeric alignment sequence representing the genetic variation type as insertion.

8. The method for analyzing insertion-deletion mutation data according to claim 6, characterized in that: The modifying of the CIGAR value and position corresponding to the representative alignment in the chimeric alignment sequence based on the gene variation type and variation length, the read length of the chimeric alignment sequence, the theoretical length of the target fragment, the length and position of the Match fragment corresponding to the representative alignment and the substitute alignment comprises: If the classification result of the representative alignment sequence in the chimeric alignment sequence is the MS class, the starting position remains unchanged; if the variation type is deletion, the CIGAR value is modified from L1ML2S to L1ML 缺失 D(L 序列 -L1)M; if the mutation type is insertion, the CIGAR value is changed from L1ML2S to L1ML 插入 I(L 目标 -L1)M; If the classification result of the representative alignment sequence in the chimeric alignment sequence is the SM class, the starting position is replaced with the position of the alternative alignment MS type; if the variation type is deletion, the CIGAR value is modified from L2SL1M to L3ML 缺失 D(L 序列 -L3)M; if the mutation type is insertion, the CIGAR value is changed from L2SL1M to L3ML 插入 I(L 目标 -L3)M; Wherein, L1 represents the length of the successfully mapped segment; L2 represents the length of the Softclip segment; L3 represents the length of the successfully mapped segment of the substitute comparison; L 缺失 is the missing length; L 插入 is the insertion length; L 序列 represents the actual length of the chimeric alignment sequence; L 目标 Represents the theoretical length of the target fragment.

9. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: The step of generating, according to the CIGAR value of each sequence in the intersection sequence set, a target segment based on the breakpoint extension of each sequence in the intersection sequence set comprises: In the intersection sequence set, the sequences whose screening gene variation type is deletion or insertion and whose deletion or insertion length is not less than 3 are used as screening sequences; At the breakpoint position of the insertion or deletion of the screening sequence, a second preset length is extended to the left and right to serve as the target fragment.

10. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: The performing preliminary screening and elimination processing on each of the target fragments in the first Kmer data set comprises: In the first Kmer data set, a preliminary screening and elimination process is performed according to at least one of the following elimination rules: A. Eliminate the target fragments whose mapping quality of the sequence where the target fragment is located is less than a preset mapping quality threshold; B. Based on the base composition, remove the target fragments with continuous repeated bases of not less than 12 bp; C. Determine the target fragment to be removed according to the comparison information of the sequence where the target fragment is located.

11. The method for analyzing insertion-deletion mutation data according to claim 10, characterized in that: The step of determining the target fragment to be removed according to the comparison information of the sequence where the target fragment is located comprises: Obtain the specific location of the sequence breakpoint where the target fragment is located in the reference genome; According to the gene variation type and the specific position, the start and end positions of the breakpoint are calculated respectively; the calculation method includes: If the gene variation type is insertion, the start and end positions of the breakpoint include: Starting position: P 插入起始 =P 具体 -L 插入长度 ; End position: P 插入终止 =P 具体 +L 插入长度 ; If the gene variation type is deletion, the start and end positions of the breakpoints include: Starting position: P 缺失起始 =P 具体 ; End position: P 缺失终止 =P 具体 +L 缺失长度 ; Among them, P 插入起始 Represents the starting position of the breakpoint where the gene variation type is insertion; P 插入终止 Represents the end position of the breakpoint where the gene variation type is insertion; P 缺失起始 Represents the starting position of the breakpoint where the gene variation type is deletion; P 缺失终止 Represents the end position of the breakpoint where the gene variation type is deletion; P 具体 Represents the specific location; L 插入长度 represents the length of the inserted sequence; L 缺失长度 represents the length of the missing sequence; If the interval where the start and end positions of the breakpoint are located does not intersect with the preset target capture interval, the target segment corresponding to the breakpoint is removed.

12. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: After performing preliminary screening and elimination processing on each of the target fragments in the first Kmer data set to obtain a second Kmer data set, the method further includes: Based on the second Kmer data set, according to the length of the inserted or deleted sequence in each of the target fragments, the sequence is extended to the breakpoints at both ends to generate Extend fragments, and a first Extend fragment data set is constructed; In the first Extend segment data set, one or two representative Extend segments are obtained for each bed interval as target segments after screening of the interval, and the second Extend segment data set containing the target segments is generated.

13. The method for analyzing insertion-deletion mutation data according to claim 12, characterized in that: The method for calculating the extension length extending to the breakpoints at both ends of the sequence includes: If the gene variation type corresponding to the target fragment is insertion, the length of the target fragment extending to the left and right is: L 延伸 =C 可调 +L I ×2; If the gene variation type corresponding to the target segment is deletion, the length of the target segment extending to the left and right is: L 延伸 =C 可调 ; Among them, L 延伸 Represents the length of the target segment extending to the left and right, and does not include the Softclip area; C 可调 Represents preset adjustable parameters; L I Represents the length of the inserted fragment in the CIGAR value.

14. The method for analyzing insertion-deletion mutation data according to claim 12, characterized in that: In the first Extend segment data set, one or two representative Extend segments are obtained for each bed interval as target segments after screening of the interval, and the second Extend segment data set containing the target segments is generated, including: The first Extend fragment data set is classified according to the length of the insertion or deletion sequence; wherein the sequence with the length of the insertion or deletion sequence not less than 6 is regarded as the E1a data set; and the sequence with the length of the insertion or deletion sequence less than 6 is regarded as the E1b data set; The Extend fragments in the classified E1a dataset and the E1b dataset are sorted in order according to the number of Extend fragments, the number of target fragments, the quality of sequence alignment, and the length of the inserted or deleted sequence, and the first-ranked Extend fragment in the E1a dataset and the E1b dataset is selected as the best representative Extend fragment; Comparing the similarity and the difference in the number of Extend segments between the longest common subsequences of the optimal representative Extend segments of the E1a data set and the E1b data set, and obtaining a comparison result; including: A. If the similarity is greater than 80%, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set; B. If the similarity is not greater than 80%, and the number of Extend segments is less than 3, and the difference in the number of Extend segments is more than 30 times, the comparison result is: retaining the data set with the largest number of Extend segments to generate the second Extend segment data set; C. If both A and B do not match, the comparison result is: both the E1a data set and the E1b data set are retained to generate the second Extend segment data set.

15. The insertion / deletion mutation data analysis method according to claim 1, characterized in that: The re-aligning of the target fragments based on the preset targeted capture interval to obtain an alignment result includes: Extending the preset targeted capture interval to the left and right by a third preset length respectively to obtain a second bed interval; Realigning the target fragment with the second bed interval to obtain the alignment result; The scoring rule expression in the re-alignment algorithm is: [match:1,mismatch:1,gap open:2,gap extend:0]; Among them, Match: 1 means that if the two aligned bases are the same, 1 point is obtained; Mismatch: -1 means that if the two aligned bases are different, 1 point is deducted; Gap Open: -2 means that if a gap is introduced, 2 points are deducted; Gap Extend: 0 means that if a gap already exists, no point is deducted.

16. An insertion-deletion mutation data analysis device, characterized in that: include: An alignment intersection module is used to obtain a sequence alignment file corresponding to a target sample, and extract an intersection sequence existing after the sequence alignment file is aligned with the first bed interval of a preset target capture interval, to obtain an intersection sequence set; A preliminary screening module is used to generate a target fragment based on the breakpoint extension of each sequence in the intersection sequence set according to the CIGAR value of each sequence in the intersection sequence set, so as to obtain a first Kmer data set containing the target fragment; and perform preliminary screening and elimination processing on each target fragment in the first Kmer data set to obtain a second Kmer data set; A re-alignment module, used for re-aligning the target fragment based on the preset target capture interval to obtain an alignment result; A result generation module is used to obtain mutation site information according to the comparison result, and obtain a variation set based on the mutation site information.

17. An insertion-deletion mutation data analysis system, characterized in that: It comprises a memory and a processor, wherein the memory stores an insertion-deletion mutation data analysis program, and the processor runs the insertion-deletion mutation data analysis program so that the insertion-deletion mutation data analysis system performs the insertion-deletion mutation data analysis method according to any one of claims 1 to 15.

18. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores an insertion-deletion mutation data analysis program, and when the insertion-deletion mutation data analysis program is executed by the processor, the insertion-deletion mutation data analysis method according to any one of claims 1 to 15 is implemented.

Images