A method and system for detecting splicing interval variation based on allele perception

Abstract

The present invention belongs to the technical field of genome structural variation detection, and specifically relates to a splice interval variation detection method and system based on allele perception. The method first calculates an effective high splice signal interval, and then marks the splice interval, so that the variation in the splice interval can be more comprehensively detected. Before the variation feature comparison, for the same variation in the splice interval, which is manifested as the simultaneous existence of an insertion or duplication variation signal and a splice signal, or the simultaneous existence of a deletion variation signal and a splice signal, an insertion signal and a deletion signal are constructed, and the signals are unified to make the clustering result more accurate. According to the clustering result in the splice interval, a consensus sequence is constructed, and the consensus sequence is aligned to a reference genome. The minimap2 alignment result, the new alignment result obtained by re-alignment, or the read segment alignment information and other multi-level information are comprehensively applied to comprehensively detect various types of variations and complex heterozygous variations, so as to improve the detection capability of structural variations.

Description

Technical Field

[0001] The present invention belongs to the technical field of genome structural variation detection, and specifically relates to a splicing interval variation detection method and system based on allele perception. Background Art

[0002] In genomics, the main technical problem for variation detection based on clipping intervals is to accurately identify and analyze the location and characteristics of complex structural variations (such as large insertions or deletions, inversions, and duplications). These variants are located in areas of high clipping signals, and traditional alignment algorithms tend to ignore or misinterpret such fragments. In addition, the diversity of the characteristics of variation events (such as length uncertainty, directionality, position overlap, etc.) increases the complexity of identification, resulting in certain limitations in existing methods in variation detection, such as the omission of complex variations and inaccurate identification results. Summary of the invention

[0003] The present invention provides a method and system for detecting splicing interval variation based on allele perception.

[0004] The technical solution of the present invention is as follows:

[0005] The present invention provides a method for detecting splicing interval variation based on allele perception, comprising the following steps:

[0006] S1: Align the sequencing data with the reference genome, analyze the alignment information, extract the variation signal of each base position, and obtain the effective high shearing signal interval;

[0007] Based on the effective high shearing signal interval, all alignment fragments of a read segment are obtained, and according to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked;

[0008] S2: Based on the splicing interval position and splicing interval variation type, the variation features of the splicing intervals of the two reads are extracted respectively. For the same variation, if the insertion or duplication variation signal and the splicing signal exist at the same time, the insertion signal is constructed; if the deletion variation signal and the splicing signal exist at the same time, the deletion signal is constructed;

[0009] Sorting by genomic coordinates, constructing variant feature sequences, clustering them, and obtaining cluster sets;

[0010] S3: All reads in the cluster set are locally spliced ​​to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splicing interval to obtain the alignment result;

[0011] If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval;

[0012] If the alignment result does not contain the shear interval, the shear interval is expanded to both sides, and the new alignment result is processed to obtain the mutated base sequence;

[0013] S4: According to the merging conditions, the shear interval mutations are merged to determine the effective mutations;

[0014] Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.

[0015] The step S1, obtaining an effective high shear signal interval, is specifically:

[0016] After obtaining the shearing signal according to the variation signal, the proportion of reads with shearing signals in the to-be-tested interval is calculated, the to-be-tested interval is divided according to the sliding step, and the proportion of reads with shearing signals in the sliding step interval is calculated. If the proportion of reads with shearing signals in the sliding step interval is greater than the preset proportion, the sliding step interval is recorded as the target shearing interval. According to the sliding sequence, the obtained shearing intervals are sequentially recorded as the first target shearing interval, the second target shearing interval, ..., the Nth target shearing interval, and the interval between the start position of the first target shearing interval and the end position of the Nth target shearing interval is taken as the high shearing signal interval;

[0017] Obtaining the alignment fragments and supplementary alignment fragments of the high shearing signal interval, calculating the number of non-supplementary alignment fragments whose alignment fragment length is greater than a preset length, if the number of non-supplementary alignment fragments is not less than the preset number of supporting reads, and the ratio of the number of non-supplementary alignment fragments to the number of all non-supplementary alignment fragments is not less than the preset ratio, then recording the shearing signal interval as a valid high shearing signal interval;

[0018] If the total number of aligned fragments is not less than the first number and not greater than the second number, and the ratio of the number of non-complementary aligned fragments to the total number of non-complementary aligned fragments is not less than a preset ratio, the shearing signal interval is recorded as a valid high shearing signal interval.

[0019] The step S1, according to the shearing signal at one end of the comparison segment, obtains the comparison segment adjacent to the other end of the comparison segment, and marks the shearing interval position and shearing interval variation type, specifically:

[0020] Compare the directions, breakpoint positions and overlaps of the two aligned fragments, obtain the variation types of the aligned fragments and record the number of reads; until all aligned fragments are processed, obtain the number of reads of insertion or deletion variations greater than the preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal interval;

[0021] Based on the variation types of the aligned fragments, the number of reads of insertion or deletion variations greater than a preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal intervals are compared, and the shearing interval position and shearing interval variation type are marked.

[0022] Preferably, the directions, breakpoint positions and overlaps of the two aligned fragments are compared to obtain the variation types of the aligned fragments, specifically:

[0023] If the two aligned fragments have the same direction, no overlap, and there is an insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints belong to the same variation interval, and the variation type of the aligned fragment is an insertion or deletion variation;

[0024] If the two aligned fragments have the same direction and overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragment is a duplication variation;

[0025] If the directions of the two aligned fragments are opposite, there is no overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragments is an inversion variation.

[0026] Preferably, based on the variation type of the aligned fragments, the number of reads containing insertion or deletion variation greater than a preset size, the number of reads of duplication variation, and the number of reads of inversion variation in the effective high shearing signal interval are compared, and the shearing interval position and shearing interval variation type are marked, specifically:

[0027] If the number of read segments containing insertion and deletion variations larger than a preset size in the effective high shearing signal interval is not less than a preset number, the minimum breakpoint position of the interval is taken as the starting position of the interval, and the maximum breakpoint position of the interval is taken as the ending position of the interval, which are used as the shearing interval position;

[0028] Compare the number of reads for insertion variants to the number of reads for deletion variants:

[0029] If the number of reads of the insertion variant is not less than the number of reads of the deletion variant, the variant type of the splicing interval is marked as an insertion variant;

[0030] If the number of reads for the insertion variant is less than the number of reads for the deletion variant, the variant type of the splicing interval is marked as a deletion variant;

[0031] If the effective high shear signal interval does not contain insertion or deletion mutations larger than the preset size, the average values ​​of the breakpoint positions of the left and right intervals are used as the starting and ending positions of the shear interval respectively;

[0032] Compare the number of reads for duplication variants to the number of reads for inversion variants:

[0033] If the number of reads for duplication variation is not less than the number of reads for inversion variation, the variation type of the splicing interval is marked as duplication variation;

[0034] If the number of reads for duplication variation is less than the number of reads for inversion variation, the variation type of the splicing interval is marked as inversion variation.

[0035] In step S2, if the insertion or duplication mutation signal and the shearing signal exist simultaneously, an insertion signal is constructed, specifically:

[0036] If the insertion or duplication variation signal and the shearing signal exist at the same time, the two closest fragments on the same chromosome are selected, and the reference distance and the read distance of the breakpoint are calculated. If the reference distance is less than the first threshold, and the two fragments overlap with each other, and there is a comparison fragment with a spacing greater than the preset distance between the two fragments, or, if the read distance is less than the second threshold, and the absolute value of the distance difference between the reference distance and the read distance is greater than the third threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read position at the breakpoint at the starting position, with a length of the absolute value of the distance difference between the reference distance and the read distance plus 1, the base sequence in the read is taken as the sequence base information of the variation, and the insertion signal is constructed;

[0037] If the deletion variation signal and the shearing signal exist at the same time, the deletion signal is constructed as follows:

[0038] If the deletion variation signal and the splicing signal exist at the same time, the two closest fragments on the same chromosome are selected. If the two fragments are in the same direction, and the read distance is greater than the fourth threshold, and the reference distance is less than the fifth threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read position at the breakpoint of the starting position, the base sequence in the reference genome is taken as the sequence base information of the variation with a length of 1 plus the absolute value of the difference between the reference distance and the read distance to construct a deletion signal.

[0039] The step S2 constructs the variant feature sequences respectively according to the genome coordinates, clusters them, and obtains a cluster set, specifically:

[0040] The variation characteristics of the splicing interval are divided into two categories, the first type and the second type, wherein the variation signal of the first type is an insertion or deletion, and the variation signal of the second type is a splicing signal;

[0041] Sorting by genomic coordinates, constructing variant feature sequences respectively, performing feature matching according to matching conditions, calculating feature matching degrees, and clustering allele-aware variant features of read segments in the splicing interval according to the feature matching degrees to obtain a cluster set;

[0042] The matching conditions are specifically:

[0043] For the first type, if two variants are of the same type, the variant length ratio is not less than the sixth threshold, the reference distance is not greater than the seventh threshold, and the sequence identity is not less than the eighth threshold, then the two variants are considered to match;

[0044] For the second type, if the two variants are of the same type, and the absolute value of the distance difference between the breakpoint positions is not higher than the ninth threshold, and the two variants The ratio of the second largest value to the maximum value is not less than the tenth threshold, where Indicates the reference distance between the breakpoint of the same read and its next adjacent breakpoint, and the two variants are considered to match.

[0045] In step S3, if the comparison result does not contain the shear interval, the shear interval is expanded to both sides, and the new comparison result is processed to obtain the mutated base sequence, specifically:

[0046] If the alignment result does not include the shear interval, the shear interval is extended to both sides, and the read segments of the extended interval are obtained according to the clustering information. The consensus sequence is constructed by applying abPOA, and the consensus sequence is aligned to the reference sequence of the shear interval to obtain a new alignment result.

[0047] If the new alignment result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval;

[0048] If the new alignment result does not include the shear interval, the first read that can completely include the shear interval read information is obtained based on all the read information in the cluster set, and the sequence of the read located in the shear interval is extracted as the mutated base sequence.

[0049] In step S3, if the comparison result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval, specifically:

[0050] If the comparison result includes a clipping interval, then determine the comparison block at the start position and the comparison block at the end position of the clipping interval;

[0051] Traverse the comparison information of the comparison block at the starting position of the clipping interval to determine the matching block before the starting position and closest to the starting position;

[0052] Traverse the comparison information of the comparison block at the end position of the clipping interval, and determine the matching block after the end position and closest to the end position;

[0053] According to the comparison information of the two matching blocks, the starting and ending positions of the consensus sequence located in the shear interval are calculated, and the number of reads supporting the variation is recorded. If the number of reads supporting the variation is greater than the preset support number, the variation type marked by the shear interval is determined as the detection type; according to the variation position information, the variation base sequence is extracted.

[0054] The present invention also provides a splicing interval variation detection system based on allele perception, comprising:

[0055] Shear interval position and shear interval variation type marking module: used to compare sequencing data with the reference genome, analyze the comparison information, extract the variation signal of each base position, and obtain the effective high shear signal interval;

[0056] Based on the effective high shearing signal interval, all alignment fragments of a read segment are obtained, and according to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked;

[0057] Variation feature clustering module: Based on the shear interval position and shear interval variation type, the variation features of the shear intervals of the two reads are extracted respectively. For the same variation, if the insertion or duplication variation signal and the shear signal exist at the same time, the insertion signal is constructed; if the deletion variation signal and the shear signal exist at the same time, the deletion signal is constructed;

[0058] Sorting by genomic coordinates, constructing variant feature sequences, clustering them, and obtaining cluster sets;

[0059] Splice interval variation detection module: used to perform local splicing of all reads in the clustering set to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splice interval to obtain the alignment result;

[0060] If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval;

[0061] If the alignment result does not contain the shear interval, the shear interval is expanded to both sides, and the new alignment result is processed to obtain the mutated base sequence;

[0062] Genotype calculation module: according to the merging conditions, the shear interval mutations are merged to determine the effective mutations;

[0063] Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.

[0064] Beneficial effects: The present invention first calculates the effective high shearing signal interval, and then marks the shearing interval, so as to more comprehensively detect the variation in the shearing interval; before the variation feature comparison, for the same variation in the shearing interval, which may be manifested as the simultaneous existence of an insertion or duplication variation signal and a shearing signal, or the simultaneous existence of a deletion variation signal and a shearing signal, an insertion signal and a deletion signal are constructed, and the signals are unified to make the clustering result more accurate; according to the clustering result in the shearing interval, a consensus sequence is constructed, and the consensus sequence is aligned to the reference genome, and the minimap2 alignment result, the new alignment result obtained by re-alignment, or the read segment alignment information and other multi-level information are comprehensively applied to comprehensively detect various types of variations and complex heterozygous variations, thereby improving the detection capability of structural variations. DETAILED DESCRIPTION

[0065] The following examples are intended to illustrate the present invention rather than to further limit the present invention.

[0066] The present invention takes the human genome as the research object and studies the variation detection in the genome clipping interval of the third-generation sequencing data, aiming to solve the detection of duplication, inversion and large-size variation in genome variation detection.

[0067] The present invention provides a method for detecting splice interval variation based on allele perception, comprising the following steps:

[0068] S1: Align the sequencing data with the reference genome, analyze the alignment information, extract the variation signal of each base position, and obtain the effective high shearing signal interval;

[0069] Based on the effective high shear signal interval, all the alignment fragments of a read segment are obtained, and according to the shear signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shear interval position and shear interval variation type are marked.

[0070] For large-scale structural variations, the read alignment information is usually broken into several segments, which appear as clipping signals. Therefore, the genomic interval containing a large number of clipping signals is called a clipping interval.

[0071] There are usually a large number of clipping signals or insertion or deletion variation signals of a large length in the clipping interval, which means that there is a large-scale structural variation in the interval. In step S1, the present invention can more comprehensively detect the variation in the clipping interval by calculating the effective clipping signal interval and then calculating the clipping interval, and will not miss the large-scale insertion or deletion variation in the interval.

[0072] Step S1, the operation is as follows:

[0073] 1. Align the sequencing data with the reference genome, analyze the alignment information, extract the variation signal of each base position, and obtain the effective high shear signal interval.

[0074] First, each chromosome is detected one by one according to the order of chromosomes in the genome, and the genomic regions containing mutations are calculated. In order to finely describe the alignment information of each base position in the reference genome, the base coverage is defined as the following five-tuple ( Occ A , Occ C , Occ G , Occ T , Occ N ),in Occ X Indicates the base in the alignment data aligned to this base position The number of

[0075] Then, based on the alignment information to the reference sequence, the following variation signals are extracted for each base position of the reference sequence: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive clipping data, and excessive short insertions or deletions. The above variation signals are described in detail in the "Insertion or deletion variation detection method and system based on allele perception" with publication number CN118969073A, and will not be repeated here.

[0076] Afterwards, to reduce memory overhead, each chromosome is divided into many 1Mb intervals, which are called blocks or detection blocks. For the detection of each block, the alignment data in the block is first loaded, and then a sliding window of 500 bp is used to detect the block in order from 5' to 3' of the reference genome, and the genomic region with high clipping signal is calculated. The calculation method is as follows:

[0077] After obtaining the shearing signal according to the variation signal, the proportion of reads with shearing signals in the inspection interval is calculated, the inspection interval is divided according to the sliding step, and the proportion of reads with shearing signals in the sliding step interval is calculated. If the proportion of reads with shearing signals in the sliding step interval is greater than the preset proportion, the sliding step interval is recorded as the target shearing interval. According to the sliding order, the obtained shearing intervals are recorded as the first target shearing interval, the second target shearing interval, ..., the Nth target shearing interval, and the interval between the start position of the first target shearing interval and the end position of the Nth target shearing interval is used as the high shearing signal interval.

[0078] The determined high clipping signal interval indicates that the interval contains too much clipping data (for example, the proportion of read segments with clipping signals is 20% or more). For example, it can be:

[0079] Each base position p in the 500 bp detection window is detected in turn. If the proportion of reads with clipping signals in the interval [p, p+100] is 20% or more, the position i where the high clipping interval first appears is recorded as the starting position of the clipping interval, and the position j where the high clipping interval last appears is recorded as the ending position of the clipping interval. The interval [i, j] is the high clipping signal interval.

[0080] After the high clipping signal interval is determined, the insertion or deletion variation interval overlapping with it is deleted.

[0081] Next, the alignment fragments and the supplementary alignment fragments of the high shear signal interval are obtained, and the number of non-supplementary alignment fragments whose alignment fragment length is greater than the preset length is calculated. If the number of non-supplementary alignment fragments is not less than the preset number of supporting reads, and the ratio of the number of non-supplementary alignment fragments to the number of all non-supplementary alignment fragments is not less than the preset ratio, then the shear signal interval is recorded as a valid high shear signal interval;

[0082] If the total number of aligned fragments is not less than the first number and not greater than the second number, and the ratio of the number of non-complementary aligned fragments to the total number of non-complementary aligned fragments is not less than a preset ratio, the shearing signal interval is recorded as a valid high shearing signal interval.

[0083] Regarding the relationship between alignment segments and reads: a read segment that is aligned to the reference segment is called an alignment segment. A read segment may be broken into multiple segments due to alignment problems, and each segment is an alignment segment.

[0084] The process of obtaining the effective high shear signal interval according to the high shear signal interval is as follows:

[0085] In the first step, the alignment fragments within the high clipping signal interval are obtained, and the supplementary alignment fragments (SA) are obtained according to the alignment fragment labels, and the read segment information that does not contain clipping signals and the read segment alignment information whose clipping sizes at both ends are less than 200bp are removed.

[0086] In the second step, the number of non-complementary alignment fragments (longer_seg_num) whose alignment fragment length is greater than a preset length (e.g., half of the read length) is calculated.

[0087] The third step is to obtain the effective high shearing signal interval according to the size of the number of non-complementary alignment fragments and the preset number of supporting reads, or the ratio of the number of non-complementary alignment fragments to the total number of non-complementary alignment fragments, specifically:

[0088] In one case, if the number of non-complementary alignment fragments is not less than the preset number of supporting reads, and the ratio of the number of non-complementary alignment fragments to the total number of non-complementary alignment fragments is not less than the preset ratio, then the shearing signal interval is recorded as a valid high shearing signal interval; for example, it can be:

[0089] The number of non-supplementary alignment fragments (longer_seg_num) ≥ the preset number of supported reads, and the ratio of longer_seg_num to the total number of non-supplementary alignment fragments ≥ 0.1;

[0090] In one case, if the total number of aligned fragments is not less than the first number and not greater than the second number, and the ratio of the number of non-complementary aligned fragments to the total number of non-complementary aligned fragments is not less than a preset ratio, then the shearing signal interval is recorded as a valid high shearing signal interval; for example, it can be:

[0091] If the coverage in this interval is too small, such as 1≤total number of aligned fragments≤10, and the ratio of longer_seg_num to the total number of non-supplementary aligned fragments is ≥0.1.

[0092] If neither of the above two conditions is met, the interval is considered an invalid clipping signal interval, and the data in this interval is too messy to produce high-quality mutation detection results. Therefore, this interval is excluded from the mutation detection and will not be discussed later.

[0093] 2. Based on the effective high shear signal interval, obtain all the alignment fragments of a read segment, obtain the alignment fragment adjacent to the other end of the alignment fragment according to the shear signal at one end of the alignment fragment, and mark the shear interval position and shear interval variation type.

[0094] According to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked, specifically:

[0095] Compare the directions, breakpoint positions and overlaps of the two aligned fragments, obtain the variation types of the aligned fragments and record the number of reads; until all aligned fragments are processed, obtain the number of reads of insertion or deletion variations greater than the preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal interval;

[0096] Based on the variation types of the aligned fragments, the number of reads of insertion or deletion variations greater than a preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal intervals are compared, and the shearing interval position and shearing interval variation type are marked.

[0097] Among them, the directions, breakpoint positions and overlaps of the two aligned fragments are compared to obtain the variation types of the aligned fragments, specifically:

[0098] If the two aligned fragments have the same direction, no overlap, and there is an insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints belong to the same variation interval, and the variation type of the aligned fragment is an insertion or deletion variation;

[0099] If the two aligned fragments have the same direction and overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragment is a duplication variation;

[0100] If the directions of the two aligned fragments are opposite, there is no overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragments is an inversion variation.

[0101] If two breakpoints do not belong to the same variation interval, the two adjacent breakpoints are recorded in different intervals according to the size of the breakpoint positions.

[0102] In addition, based on the variation type of the aligned fragments, the number of reads containing insertion or deletion variations greater than the preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shear signal interval are compared, and the shear interval position and shear interval variation type are marked. Specifically:

[0103] If the number of read segments containing insertion and deletion variations larger than a preset size in the effective high shearing signal interval is not less than a preset number, the minimum breakpoint position of the interval is taken as the starting position of the interval, and the maximum breakpoint position of the interval is taken as the ending position of the interval, which are used as the shearing interval position;

[0104] Compare the number of reads for insertion variants to the number of reads for deletion variants:

[0105] If the number of reads of the insertion variant is not less than the number of reads of the deletion variant, the variant type of the splicing interval is marked as an insertion variant;

[0106] If the number of reads for the insertion variant is less than the number of reads for the deletion variant, the variant type of the splicing interval is marked as a deletion variant;

[0107] If the effective high shear signal interval does not contain insertion or deletion mutations larger than the preset size, the average values ​​of the breakpoint positions of the left and right intervals are used as the starting and ending positions of the shear interval respectively;

[0108] Compare the number of reads for duplication variants to the number of reads for inversion variants:

[0109] If the number of reads for duplication variation is not less than the number of reads for inversion variation, the variation type of the splicing interval is marked as duplication variation;

[0110] If the number of reads for duplication variation is less than the number of reads for inversion variation, the variation type of the splicing interval is marked as inversion variation.

[0111] S2: Based on the splicing interval position and splicing interval variation type, the variation features of the splicing intervals of the two reads are extracted respectively. For the same variation, if the insertion or duplication variation signal and the splicing signal exist at the same time, the insertion signal is constructed; if the deletion variation signal and the splicing signal exist at the same time, the deletion signal is constructed;

[0112] The variant feature sequences are constructed and clustered according to the genomic coordinate sorting to obtain a cluster set.

[0113] In step S2, the insertion signal and deletion signal construction method provided by the present invention can better handle more complex situations in the clipping interval, such as the situation where large-sized insertions or deletions and clipping signals coexist in an interval. Unifying the signals can make the clustering results more accurate and correctly detect large-sized insertions or deletions in the clipping interval.

[0114] Clustering of splice interval allele-aware variant features is performed as follows:

[0115] 1. Based on the splicing interval position and splicing interval variation type, the variation features of the splicing intervals of the two reads are extracted respectively. For the same variation, if the insertion or duplication variation signal and the splicing signal exist at the same time, the two closest fragments on the same chromosome are selected, and the reference distance of the breakpoint and the read distance are calculated. If the first preset condition is met, the insertion signal is constructed according to the breakpoint position; if the deletion variation signal and the splicing signal exist at the same time, the two closest fragments on the same chromosome are selected. If the second preset condition is met, the deletion signal is constructed according to the breakpoint position.

[0116] The same variant in the clipping interval sometimes manifests as the coexistence of an insertion or duplication variant signal and a shearing signal, or the coexistence of a deletion variant signal and a shearing signal. For example, the same variant manifests both 0 and 2, or both 1 and 2 (0 is the insertion type, 1 is the deletion type, and 2 is the clipping type). Therefore, the different manifestations of the same variant need to be unified for easy comparison.

[0117] Preferably, if the insertion or duplication variation signal and the shearing signal exist simultaneously, the two closest fragments on the same chromosome are selected, the reference distance and the read distance of the breakpoint are calculated, and if the first preset condition is met, the insertion signal is constructed according to the breakpoint position, specifically:

[0118] If the insertion or duplication variation signal and the shearing signal exist at the same time, the two closest fragments on the same chromosome are selected, and the reference distance and the read distance of the breakpoint are calculated. If the reference distance is less than the first threshold, and the two fragments overlap with each other, and there is a comparison fragment with a spacing greater than the preset distance between the two fragments, or, if the read distance is less than the second threshold, and the absolute value of the difference between the reference distance and the read distance is greater than the third threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read position at the breakpoint at the starting position, with a length of the absolute value of the difference between the reference distance and the read distance plus 1, the base sequence in the read is taken as the sequence base information of the variation, and an insertion signal is constructed.

[0119] Among them, the first preset condition is:

[0120] The reference distance is less than a first threshold value (e.g., 500 bp), the two fragments overlap each other, and there is a comparison fragment with a spacing greater than a preset distance between the two fragments;

[0121] Or, the read distance is less than a second threshold (such as 500 bp), and the absolute value of the distance difference between the reference distance and the read distance is greater than a third threshold (500 bp), and there is no alignment fragment with a spacing greater than a preset distance between the two fragments.

[0122] The insertion signal is expressed as, , indicating that the end position of the inserted signal is the same as the start position; Indicates that the mutation type is insertion; The length that represents the absolute value of the difference between the reference distance and the read distance plus 1. To take the absolute value operation; It means starting from the read position at the start position breakpoint, the base sequence in the read is taken as the sequence base information of the variant with the absolute value of the distance difference between the reference distance and the read distance plus 1.

[0123] In addition, if the deletion variation signal and the splicing signal exist at the same time, the two closest fragments on the same chromosome are selected. If the second preset condition is met, the deletion signal is constructed according to the breakpoint position, specifically:

[0124] If the deletion variation signal and the splicing signal exist at the same time, the two closest fragments on the same chromosome are selected. If the two fragments are in the same direction, and the read distance is greater than the fourth threshold, and the reference distance is less than the fifth threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read position at the breakpoint of the starting position, the base sequence in the reference genome is taken as the sequence base information of the variation with a length of 1 plus the absolute value of the difference between the reference distance and the read distance to construct a deletion signal.

[0125] Among them, the second preset condition is:

[0126] The two fragments have the same direction, and the read distance is greater than a fourth threshold (such as 500 bp), and the reference distance is less than a fifth threshold (such as 500 bp), and there is no alignment fragment with a spacing greater than a preset distance between the two fragments.

[0127] The missing signal is expressed as, , which means the end position is the absolute value of the difference between the start position plus the reference distance and the read distance; Indicates that the variant type is deletion; The length that represents the absolute value of the difference between the reference distance and the read distance plus 1. To take the absolute value operation; It means starting from the read position at the start position breakpoint, the length of the absolute value of the distance difference between the reference distance and the read distance plus 1 is used to take the base sequence in the reference genome as the sequence base information of the variant.

[0128] 2. Sort by genomic coordinates to construct variant feature sequences, cluster them, and obtain cluster sets.

[0129] Specifically:

[0130] The variation characteristics of the splicing interval are divided into two categories, the first type and the second type, wherein the variation signal of the first type is an insertion or deletion, and the variation signal of the second type is a splicing signal;

[0131] The variant feature sequences are constructed respectively by sorting the genomic coordinates, feature matching is performed according to the matching conditions, the feature matching degree is calculated, and the allele-aware variant feature clustering is performed on the read segments of the splicing interval according to the feature matching degree to obtain a cluster set.

[0132] The variation characteristics of the shear interval are divided into two categories, the first type and the second type, wherein the variation signal of the first type is an insertion or deletion, and the variation signal of the second type is a shear signal. The specific steps are as follows:

[0133] In order to distinguish different allele variations in the clipping interval, the variation characteristics in the clipping interval are divided into two categories: the first type (denoted as Indel type) and the second type (denoted as clipping type). The variation characteristics of the Indel type are defined as a sextuple , the mutation characteristics of the clipping type are defined as a four-tuple ,in, is the name of the chromosome where the mutation is located; and Chromosome The starting and ending positions on the is the variable length; is the sequence base information of the variation; is the reference distance between the breakpoint of the same read segment and its next adjacent breakpoint, and the last breakpoint The value is 0; is the type of mutation signal: 0 is insertion type, 1 is deletion type, and 2 is clipping type, representing duplication or inversion.

[0134] In addition, variant feature sequences are constructed by sorting the genomic coordinates, feature matching is performed according to the matching conditions, the feature matching degree is calculated, and the read segments of the splicing interval are clustered according to the allele-aware variant features to obtain a cluster set. The specific steps are as follows:

[0135] The two variant feature sequences can be expressed as: , .

[0136] Preferably, the matching condition is specifically:

[0137] For the first type (Indel type), if the two variants are of the same type, and the variant length ratio is not less than the sixth threshold (such as 0.7), and the reference distance is not greater than the seventh threshold (such as (500 bp), and the sequence identity is not less than the eighth threshold (such as 0.8), then the two variants are considered to match;

[0138] For the second type (clipping type), if the two variants are of the same type, and the absolute value of the distance difference between the breakpoint positions is not higher than the ninth threshold (such as 500 bp), and the two variants The ratio of the second largest value to the maximum value is not less than the tenth threshold (such as 0.9), where Indicates the reference distance between the breakpoint of the same read and its next adjacent breakpoint, and the two variants are considered to match.

[0139] The match score is calculated as: match +2, mismatch -5, gap extension -2, gap opening -4.

[0140] Then, a variation feature comparison matrix is ​​constructed based on the matching scores. According to the variation feature comparison matrix, the variation feature sequences of the two reads are compared to obtain a matching vector; the feature matching degree MR is calculated based on the matching vector. The value range of MR is 0-1. In particular, when the length of the variation feature sequence is 0, it means that neither of the two reads contains variation features and is filtered. The process of calculating the feature matching degree is described in detail in the "Insertion or deletion variation detection method and system based on allele perception" with publication number CN118969073A, which will not be repeated here.

[0141] Finally, the reads in the shearing interval are clustered according to the allele-aware variation characteristics based on the feature matching degree to obtain a cluster set. The first is to initialize the cluster read set, and its goal is to find two reads with different variation characteristics as the initial cluster sets A and B. Secondly, after initializing A and B, the read set R is clustered according to the allele-aware variation characteristics. The clustering results are saved in a file, and the information includes the number of reads and the name of the reads. This process is described in detail in the "Allele-aware Insertion or Deletion Variation Detection Method and System" with publication number CN118969073A, so I will not go into details here.

[0142] S3: All reads in the cluster set are locally spliced ​​to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splicing interval to obtain the alignment result;

[0143] If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval;

[0144] If the alignment result does not include the shear interval, the shear interval is extended to both sides, and the new alignment result is processed to obtain the mutated base sequence.

[0145] The variation in the clipping interval is more complex, and the detection of variation is more challenging. In step S3, the present invention constructs a consensus sequence using different methods for different situations based on the clustering results in the clipping interval, and aligns the consensus sequence to the reference genome, and comprehensively applies the minimap2 alignment results, the new alignment results obtained by re-alignment, or the multi-level information such as the read segment alignment information, which can comprehensively detect various types of variations and complex heterozygous variations, and improve the detection ability of structural variations.

[0146] The mutation detection of each clipping interval is calculated independently, including the mutation position, mutation sequence, and subsequent genotype calculation, and different intervals are not affected. The operation is as follows:

[0147] 1. All reads in the cluster set are locally spliced ​​to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the shear interval to obtain the alignment result.

[0148] The specific method is as follows: the read sequence information after clustering is saved in the FastA file, including the description information and base sequence information of each read. Then, the wtdbg2 tool is used to locally splice the sequences in the read set to obtain the splicing result (contig), i.e. the consensus sequence, which is stored in the FastA file, including the description information and base sequence of each consensus sequence.

[0149] In addition, the read name information of all reads in the clustered read set is further stored in the description information of each contig, so as to facilitate the subsequent calculation of the genotype of the detected variation.

[0150] Since wtdbg2 has limited sequence splicing capabilities in complex intervals, some complex intervals may not be able to successfully construct the splicing results of variant intervals. Therefore, if the splicing result does not contain the variant interval, it is necessary to use abPOA to reconstruct the consensus sequence of the variant interval and perform minimap2 alignment, and detect the variation based on the regenerated alignment information. If the above strategies are unsuccessful, the variation is detected directly based on the read segment information.

[0151] 2. If the alignment result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval;

[0152] If the alignment result does not include the shear interval, the shear interval is extended to both sides, and the new alignment result is processed to obtain the mutated base sequence.

[0153] If the alignment result does not contain the shear interval, the shear interval is extended to both sides, and the new alignment result is processed to obtain the mutated base sequence, specifically:

[0154] If the alignment result does not include the shear interval, the shear interval is extended to both sides, and the read segments of the extended interval are obtained according to the clustering information. The consensus sequence is constructed by applying abPOA, and the consensus sequence is aligned to the reference sequence of the shear interval to obtain a new alignment result.

[0155] If the new alignment result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval;

[0156] If the new alignment result does not include the shear interval, the first read that can completely include the shear interval read information is obtained based on all the read information in the cluster set, and the sequence of the read located in the shear interval is extracted as the mutated base sequence.

[0157] Among them, if the alignment result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval, specifically:

[0158] If the comparison result includes a clipping interval, then determine the comparison block at the start position and the comparison block at the end position of the clipping interval;

[0159] Traverse the comparison information of the comparison block at the starting position of the clipping interval to determine the matching block before the starting position and closest to the starting position;

[0160] Traverse the comparison information of the comparison block at the end position of the clipping interval, and determine the matching block that is closest to the end position;

[0161] According to the comparison information of the two matching blocks, the starting and ending positions of the consensus sequence located in the shear interval are calculated, and the number of reads supporting the variation is recorded. If the number of reads supporting the variation is greater than the preset support number, the variation type marked by the shear interval is determined as the detection type; according to the variation position information, the variation base sequence is extracted.

[0162] From the above, we can see that for the detection method of duplication and inversion mutation, the specific operation can be expressed as:

[0163] Step 1: Parse the minimap2 alignment information PAF file, traverse the CIGAR field of each alignment block (a line of the PAF file), and extract the alignment information of each alignment block;

[0164] Step 2: Traverse the alignment information of each alignment block and determine the two alignment blocks located at the start and end positions of the clipping interval (since the minimap2 alignment information of the same consensus sequence is sometimes split into multiple alignment blocks, the start and end positions of the clipping interval may correspond to different alignment blocks);

[0165] Step 3: For the matching block at the start position of the clipping interval, traverse its matching information and determine the matching block before the start position and closest to the start position, requiring the matching type of the matching block to be matching (CIGAR operation type is BAM_CMATCH or BAM_CEQUAL); similarly, determine the matching block after the end position of the clipping interval and closest to the end position;

[0166] Step 4: Based on the alignment information of the two matching blocks, calculate the start and end positions of the consensus sequence located in the clipping interval; record the number of reads supporting the variation. If the number of reads supporting the variation is greater than the preset support number, it means that the variation is successfully detected, and extract the variant base sequence based on the variation position information.

[0167] If the minimap2 alignment information does not include the clipping interval, you need to regenerate the alignment information of the mutation interval and detect the mutation based on the alignment information. The method is as follows:

[0168] Step 5: The shearing interval was extended by 1 kb on both sides, and the read segments of the extended interval were obtained according to the clustering information. The consensus sequence was constructed by abPOA, and the consensus sequence was aligned to the reference sequence by minimap2 to obtain a new alignment result.

[0169] Step 6: If the new alignment result completely contains the clipping interval, execute steps 1-4. If the variant sequence is successfully obtained, the mutation detection of the clipping interval is successful.

[0170] Step 7: If the new alignment result does not include the clipping interval, then based on all the read information in the cluster set, obtain the first read that can completely include the clipping interval read information, extract the sequence of the read in the clipping interval as the variant base sequence. If there is no read that can completely include the variant interval, the specific information of the variant cannot be determined, and the variant in the clipping interval cannot be successfully detected.

[0171] Step 8: For duplication mutations, if the mutation sequence is successfully detected, the mutation detection of the clipping interval is successful. For inversion mutations, the mutation sequence needs to be further reverse-complemented and compared with the reference sequence to calculate the sequence consistency of the two. If the sequence consistency is ≥0.8, the inversion mutation detection is successful, otherwise the inversion mutation cannot be successfully detected.

[0172] For the detection method of insertion and deletion variation, the specific operation can be expressed as:

[0173] Step 1: Parse the minimap2 alignment information PAF file, traverse the CIGAR field of each alignment block (a line of the PAF file), and extract the alignment information of each alignment block;

[0174] Step 2: Traverse the comparison information of each comparison block and determine the two comparison blocks located at the start and end positions of the clipping interval;

[0175] Step 3: For the matching block at the start position of the clipping interval, traverse its matching information and determine the matching block before the start position and closest to the start position, requiring the matching type of the matching block to be matching (CIGAR operation type is BAM_CMATCH or BAM_CEQUAL); similarly, determine the matching block after the end position of the clipping interval and closest to the end position;

[0176] Step 4: Based on the alignment information of the two matching blocks, calculate the starting and ending positions of the consensus sequence located in the clipping interval; record the number of reads supporting the variation. If the number of reads supporting the variation is greater than the preset support number, it means that the variation is successfully detected, and extract the variant base sequence based on the variation position information.

[0177] If the mutation is not successfully detected, you need to regenerate the comparison information of the mutation interval and detect the mutation based on the comparison information. The method is as follows:

[0178] Step 5: The shearing interval was extended by 2 kb on both sides, and the read segments of the extended interval were obtained based on the clustering information. A consensus sequence was constructed using abPOA, and the consensus sequence was aligned to the reference sequence using minimap2 to obtain a new alignment result.

[0179] Step 6: If the new alignment result completely contains the clipping interval, execute steps 1-4. If the variant sequence is successfully obtained, the mutation detection of the clipping interval is successful.

[0180] Step 7: If the new alignment result does not include the clipping interval, then based on all the read information in the cluster set, obtain the first read that can completely include the clipping interval read information, extract the sequence of the read in the clipping interval as the variant base sequence. If there is no read that can completely include the variant interval, the specific information of the variant cannot be determined, and the variant in the clipping interval cannot be successfully detected.

[0181] S4: Calculate genotype:

[0182] According to the merging conditions, the shear interval mutations are merged to determine the effective mutations;

[0183] Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.

[0184] The specific operations are:

[0185] Temporarily save the detected mutations in the clipping interval into the allele array;

[0186] Step 1: Sort the variants in the allele array according to the number of reads that support the variant;

[0187] Step 2: Compare any two variants in the allele array. If the merging conditions are met, that is, the following two conditions: ① The distance deviation between the starting position and the ending position of the variant is not greater than the eleventh threshold (100 bp); ② The consistency of the contig base sequence of the variant sequence is not less than the twelfth threshold (such as 0.9), then the two variants are merged. When merging, the variant with a large number of reads supporting the variant is retained, and its number of reads supporting the variant is updated to the sum of the two numbers, and the variant with a small number of reads supporting the variant is deleted.

[0188] Step 3: The first two mutations in the allele array are retained as valid mutations, and the remaining mutations are considered invalid mutations and removed.

[0189] Step 4: Calculate the genotype, which is divided into two cases according to the number of alleles detected in the splicing interval:

[0190] ① When one allele is detected in the shear interval, first use the formula: , calculate allele frequencies ;

[0191] In the formula, To support this variation The number of read segments, For variation The sequencing coverage depth at the variant position.

[0192] According to the formula: , calculate the variation Genotype ;

[0193] Among them, the genotype "1 / 1" indicates a homozygous variation, "0 / 1" indicates a heterozygous variation, and "0 / 0" indicates that the variation is not obvious and the genotype of the variation needs to be further confirmed.

[0194] ② When two alleles are detected in the splicing interval, the genotype is a heterozygous variation. According to the formula , calculate allele frequencies , and marked the genotype as 0 / 1.

[0195] Compared with the prior art, the present invention provides a variation detection method for splicing intervals based on allele perception, which first calculates the effective high splicing signal interval and then marks the splicing interval, so that the variation in the splicing interval can be detected more comprehensively; before the variation feature comparison, for the same variation in the splicing interval, which will be manifested as the simultaneous existence of an insertion or duplication variation signal and a splicing signal, or the simultaneous existence of a deletion variation signal and a splicing signal, an insertion signal and a deletion signal are constructed, and the signals are unified to make the clustering result more accurate; according to the clustering result in the splicing interval, a consensus sequence is constructed, and the consensus sequence is aligned to the reference genome, and the minimap2 alignment result, the new alignment result obtained by re-alignment, or the read segment alignment information and other multi-level information are comprehensively applied to comprehensively detect various types of variations and complex heterozygous variations, thereby improving the detection capability of structural variations.

[0196] The present invention also provides a splicing interval variation detection system based on allele perception, comprising:

[0197] Shear interval position and shear interval variation type marking module: used to compare sequencing data with the reference genome, analyze the comparison information, extract the variation signal of each base position, and obtain the effective high shear signal interval;

[0198] Based on the effective high shearing signal interval, all alignment fragments of a read segment are obtained, and according to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked;

[0199] Variation feature clustering module: Based on the shear interval position and shear interval variation type, the variation features of the shear intervals of the two reads are extracted respectively. For the same variation, if the insertion or duplication variation signal and the shear signal exist at the same time, the insertion signal is constructed; if the deletion variation signal and the shear signal exist at the same time, the deletion signal is constructed;

[0200] Sorting by genomic coordinates, constructing variant feature sequences, clustering them, and obtaining cluster sets;

[0201] Splice interval variation detection module: used to perform local splicing of all reads in the clustering set to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splice interval to obtain the alignment result;

[0202] If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval;

[0203] If the alignment result does not contain the shear interval, the shear interval is expanded to both sides, and the new alignment result is processed to obtain the mutated base sequence;

[0204] Genotype calculation module: used to calculate genotype:

[0205] According to the merging conditions, the shear interval mutations are merged to determine the effective mutations;

[0206] Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.

Claims

1. A method for detecting splice interval variation based on allele perception, characterized in that: The following steps are involved: S1: The sequencing data is compared with the reference genome, the comparison information is analyzed, the variation signal of each base position is extracted, and after the shearing signal is obtained according to the variation signal, the proportion of reads with shearing signals in the interval to be tested is calculated, the interval to be tested is divided according to the sliding step, and the proportion of reads with shearing signals in the sliding step interval is calculated. If the proportion of reads with shearing signals in the sliding step interval is greater than the preset proportion, the sliding step interval is recorded as the target shearing interval. According to the sliding order, the obtained shearing intervals are recorded as the first target shearing interval, the second target shearing interval, ..., the Nth target shearing interval, and the interval between the start position of the first target shearing interval and the end position of the Nth target shearing interval is taken as the high shearing signal interval; Obtaining the alignment fragments and supplementary alignment fragments of the high shearing signal interval, calculating the number of non-supplementary alignment fragments whose alignment fragment length is greater than a preset length, if the number of non-supplementary alignment fragments is not less than the preset number of supporting reads, and the ratio of the number of non-supplementary alignment fragments to the number of all non-supplementary alignment fragments is not less than the preset ratio, then recording the shearing signal interval as a valid high shearing signal interval; If the total number of aligned fragments is not less than the first number and not greater than the second number, and the ratio of the number of non-complementary aligned fragments to the total number of non-complementary aligned fragments is not less than a preset ratio, the shearing signal interval is recorded as a valid high shearing signal interval; Based on the effective high shearing signal interval, all alignment fragments of a read segment are obtained, and according to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked; S2: Based on the splicing interval position and the splicing interval variation type, respectively extract the variation features of the splicing intervals of the two reads. For the same variation, if the insertion or duplication variation signal and the splicing signal exist at the same time, select the two closest fragments on the same chromosome, calculate the reference distance of the breakpoint and the read distance, if the reference distance is less than the first threshold, and the two fragments overlap with each other, and there is a comparison fragment with a spacing greater than the preset distance between the two fragments, or, if the read distance is less than the second threshold, and the absolute value of the distance difference between the reference distance and the read distance is greater than the third threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then the read position at the breakpoint at the starting position is selected. Starting from the position of the read segment, the base sequence in the read segment is taken as the sequence base information of the variation with the length of the absolute value of the distance difference between the reference distance and the read segment distance plus 1, and the insertion signal is constructed; if the deletion variation signal and the shearing signal exist at the same time, the two closest fragments on the same chromosome are selected, if the two fragments are in the same direction, and the read segment distance is greater than the fourth threshold, and the reference distance is less than the fifth threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read segment position at the breakpoint of the starting position, the base sequence in the reference genome is taken as the sequence base information of the variation with the length of the absolute value of the distance difference between the reference distance and the read segment distance plus 1, and the deletion signal is constructed; Sorting by genomic coordinates, constructing variant feature sequences, clustering them, and obtaining cluster sets; S3: All reads in the cluster set are locally spliced ​​to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splicing interval to obtain the alignment result; If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval; If the alignment result does not contain the shear interval, the shear interval is expanded to both sides, and the new alignment result is processed to obtain the mutated base sequence; S4: According to the merging conditions, the shear interval mutations are merged to determine the effective mutations; Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.

2. The method for detecting splicing interval variation based on allele perception according to claim 1, characterized in that: The step S1, according to the shearing signal at one end of the comparison segment, obtains the comparison segment adjacent to the other end of the comparison segment, and marks the shearing interval position and shearing interval variation type, specifically: Compare the directions, breakpoint positions and overlaps of the two aligned fragments, obtain the variation types of the aligned fragments and record the number of reads; until all aligned fragments are processed, obtain the number of reads of insertion or deletion variations greater than the preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal interval; Based on the variation types of the aligned fragments, the number of reads of insertion or deletion variations greater than a preset size, the number of reads of duplication variations, and the number of reads of inversion variations in the effective high shearing signal intervals are compared, and the shearing interval position and shearing interval variation type are marked.

3. The method for detecting splice interval variation based on allele perception according to claim 2, characterized in that: Compare the directions, breakpoint positions and overlaps of the two aligned fragments to obtain the variation types of the aligned fragments, specifically: If the two aligned fragments have the same direction, no overlap, and there is an insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints belong to the same variation interval, and the variation type of the aligned fragment is an insertion or deletion variation; If the two aligned fragments have the same direction and overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragment is a duplication variation; If the directions of the two aligned fragments are opposite, there is no overlap, and there is no insertion or deletion variation greater than the preset size between the two breakpoints, then the two breakpoints do not belong to the same variation interval, and the variation type of the aligned fragments is an inversion variation.

4. The method for detecting splicing interval variation based on allele perception according to claim 2, characterized in that: Based on the variation type of the aligned fragments, the number of reads of insertion or deletion variation, the number of reads of duplication variation, and the number of reads of inversion variation containing effective high shear signal intervals greater than the preset size are compared, and the shear interval position and shear interval variation type are marked, specifically: If the number of read segments containing insertion and deletion variations larger than a preset size in the effective high shearing signal interval is not less than a preset number, the minimum breakpoint position of the interval is taken as the starting position of the interval, and the maximum breakpoint position of the interval is taken as the ending position of the interval, which are used as the shearing interval position; Compare the number of reads for insertion variants to the number of reads for deletion variants: If the number of reads of the insertion variant is not less than the number of reads of the deletion variant, the variant type of the splicing interval is marked as an insertion variant; If the number of reads for the insertion variant is less than the number of reads for the deletion variant, the variant type of the splicing interval is marked as a deletion variant; If the effective high shear signal interval does not contain insertion or deletion mutations larger than the preset size, the average values ​​of the breakpoint positions of the left and right intervals are used as the starting and ending positions of the shear interval respectively; Compare the number of reads for duplication variants to the number of reads for inversion variants: If the number of reads for duplication variation is not less than the number of reads for inversion variation, the variation type of the splicing interval is marked as duplication variation; If the number of reads for duplication variation is less than the number of reads for inversion variation, the variation type of the splicing interval is marked as inversion variation.

5. The method for detecting splice interval variation based on allele perception according to claim 1, characterized in that: The step S2 constructs the variant feature sequences respectively according to the genome coordinates, clusters them, and obtains a cluster set, specifically: The variation characteristics of the splicing interval are divided into two categories, the first type and the second type, wherein the variation signal of the first type is an insertion or deletion, and the variation signal of the second type is a splicing signal; Sorting by genomic coordinates, constructing variant feature sequences respectively, performing feature matching according to matching conditions, calculating feature matching degrees, and performing allele-aware variant feature clustering on the read segments of the splicing interval according to the feature matching degrees to obtain a cluster set; The matching conditions are specifically: For the first type, if two variants are of the same type, the variant length ratio is not less than the sixth threshold, the reference distance is not greater than the seventh threshold, and the sequence identity is not less than the eighth threshold, then the two variants are considered to match; For the second type, if the two variants are of the same type, and the absolute value of the distance difference between the breakpoint positions is not higher than the ninth threshold, and the two variants The ratio of the second largest value to the maximum value is not less than the tenth threshold, where Indicates the reference distance between the breakpoint of the same read and its next adjacent breakpoint, and the two variants are considered to match.

6. The method for detecting splice interval variation based on allele perception according to claim 1, characterized in that: In step S3, if the comparison result does not contain the shear interval, the shear interval is expanded to both sides, and the new comparison result is processed to obtain the mutated base sequence, specifically: If the alignment result does not include the shear interval, the shear interval is extended to both sides, and the read segments of the extended interval are obtained according to the clustering information. The consensus sequence is constructed by applying abPOA, and the consensus sequence is aligned to the reference sequence of the shear interval to obtain a new alignment result. If the new alignment result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval; If the new alignment result does not include the shear interval, the first read that can completely include the shear interval read information is obtained based on all the read information in the cluster set, and the sequence of the read located in the shear interval is extracted as the mutated base sequence.

7. The method for detecting splice interval variation based on allele perception according to claim 1, characterized in that: In step S3, if the comparison result includes a shear interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the shear interval, specifically: If the comparison result includes a clipping interval, then determine the comparison block at the start position and the comparison block at the end position of the clipping interval; Traverse the comparison information of the comparison block at the starting position of the clipping interval to determine the matching block before the starting position and closest to the starting position; Traverse the comparison information of the comparison block at the end position of the clipping interval, and determine the matching block that is closest to the end position; According to the comparison information of the two matching blocks, the starting and ending positions of the consensus sequence located in the shear interval are calculated, and the number of reads supporting the variation is recorded. If the number of reads supporting the variation is greater than the preset support number, the variation type marked by the shear interval is determined as the detection type; according to the variation position information, the variation base sequence is extracted.

8. A splicing interval variation detection system based on allele perception, characterized in that: include: Cleavage interval position and cleavage interval variation type marking module: used to compare sequencing data with the reference genome, analyze the comparison information, extract the variation signal of each base position, calculate the proportion of reads with cleavage signals in the interval to be tested after obtaining the cleavage signal according to the variation signal, divide the interval to be tested according to the sliding step, calculate the proportion of reads with cleavage signals in the sliding step interval, if the proportion of reads with cleavage signals in the sliding step interval is greater than the preset proportion, the sliding step interval is recorded as the target cleavage interval, and according to the sliding order, the obtained cleavage intervals are recorded as the first target cleavage interval, the second target cleavage interval, ..., the Nth target cleavage interval in sequence, and the interval between the start position of the first target cleavage interval and the end position of the Nth target cleavage interval is taken as the high cleavage signal interval; Obtaining the alignment fragments and supplementary alignment fragments of the high shearing signal interval, calculating the number of non-supplementary alignment fragments whose alignment fragment length is greater than a preset length, if the number of non-supplementary alignment fragments is not less than the preset number of supporting reads, and the ratio of the number of non-supplementary alignment fragments to the number of all non-supplementary alignment fragments is not less than the preset ratio, then recording the shearing signal interval as a valid high shearing signal interval; If the total number of aligned fragments is not less than the first number and not greater than the second number, and the ratio of the number of non-complementary aligned fragments to the total number of non-complementary aligned fragments is not less than a preset ratio, the shearing signal interval is recorded as a valid high shearing signal interval; Based on the effective high shearing signal interval, all alignment fragments of a read segment are obtained, and according to the shearing signal at one end of the alignment fragment, the alignment fragment adjacent to the other end of the alignment fragment is obtained, and the shearing interval position and shearing interval variation type are marked; Mutation feature clustering module: Based on the shear interval position and shear interval variation type, the variation features of the two read shear intervals are extracted respectively. For the same variation, if the insertion or duplication variation signal and the shear signal exist at the same time, the two closest fragments on the same chromosome are selected to calculate the reference distance of the breakpoint and the read distance. If the reference distance is less than the first threshold, and the two fragments overlap with each other, and there is a comparison fragment with a spacing greater than the preset distance between the two fragments, or, if the read distance is less than the second threshold, and the absolute value of the difference between the reference distance and the read distance is greater than the third threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then the breakpoint at the starting position is selected. Starting from the read position, the base sequence in the read is taken as the sequence base information of the variation with the length of the absolute value of the distance difference between the reference distance and the read distance plus 1 to construct an insertion signal; if the deletion variation signal and the shearing signal exist at the same time, the two closest fragments on the same chromosome are selected, if the two fragments are in the same direction, and the read distance is greater than the fourth threshold, and the reference distance is less than the fifth threshold, and there is no comparison fragment with a spacing greater than the preset distance between the two fragments, then starting from the read position at the breakpoint of the starting position, the base sequence in the reference genome is taken as the sequence base information of the variation with the length of the absolute value of the distance difference between the reference distance and the read distance plus 1 to construct a deletion signal; Sorting by genomic coordinates, constructing variant feature sequences, clustering them, and obtaining cluster sets; Splice interval variation detection module: used to perform local splicing of all reads in the clustering set to obtain a consensus sequence; the consensus sequence is aligned to the reference sequence of the splice interval to obtain the alignment result; If the alignment result includes a splicing interval, the base sequence of the insertion, deletion, duplication or inversion variation is obtained according to the positions of the two breakpoints in the splicing interval; If the alignment result does not contain the shear interval, the shear interval is expanded to both sides, and the new alignment result is processed to obtain the mutated base sequence; Genotype calculation module: according to the merging conditions, the shear interval mutations are merged to determine the effective mutations; Based on the number of alleles in the splicing interval, the genotype of the effective variant is determined.