Method and system for detecting insertion or deletion variation based on allele perception

Abstract

The present invention belongs to the technical field of genome structural variation detection, and specifically relates to an insertion or deletion variation detection method and system based on allele perception. By extracting the insertion or deletion variation characteristics of sequencing data at the variation position, a variety of matching index calculation methods are used to fully utilize the diploid variation feature information in the genome, a variation feature clustering method of different alleles is designed, and an accurate insertion or deletion variation detection method is designed based on the clustering results, thereby improving variation detection in complex regions, especially in multi-alleles, thereby improving variation detection accuracy. The detection method can be applied to complex genome regions, including those with multiple alleles and medically related genes, and improves the detection capabilities in these regions.

Description

Technical Field

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

[0002] Structural variation (SV) accounts for a large part of the entire genome sequence variation and is of great significance in human genomics and clinical medicine. There are currently two mainstream variation detection methods. One is a variation detection method based on traditional clustering, which first extracts the original variation signal by analyzing the internal alignment and split alignment of sequencing reads, and then clusters and filters the variation features to detect the variation; the other is a variation detection method based on allele-aware clustering, which first obtains structural variation feature information from the alignment data, clusters and refines the filtering based on the information, and then obtains allele information, optimizes the structural variation information based on the allele information, and thus identifies the variation.

[0003] However, these two methods also have some limitations. First, they do not fully utilize the diploid information of the human genome, and are prone to miss some complex heterozygous insertion or deletion variants; second, the existing variant detection method based on clustering variant features clusters the similarity of the size and position of the variant features, and the recognition results are inaccurate, and the accuracy of the variant breakpoints still needs to be further improved. Summary of the invention

[0004] In order to solve the problems raised by the background technology, the present invention takes the human genome as the research object and studies the detection of genome insertion or deletion variation in third-generation sequencing data, aiming to solve the problems of missing some complex heterozygous insertion or deletion variation and inaccurate variation breakpoints in genome variation detection, especially insertion or deletion variation at multi-allelic loci and some highly repetitive complex regions. To this end, the present invention provides an insertion or deletion variation detection method and system based on allele perception.

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

[0006] The present invention provides an insertion or deletion variation detection method based on allele perception, comprising the following steps:

[0007] S1: Sequencing data is aligned with the reference genome, and variation signals are extracted from the aligned data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive shearing data, and excessive short insertions or deletions;

[0008] S2: Based on the variation signal, the interval between two adjacent normal regions is calculated to obtain the variation region;

[0009] S3: performing pairwise comparisons on all reads of the variant region. If there is an intersection, extracting the variant feature sequences of the two reads, calculating the variant feature matching degree, and storing the reads in the first set if the variant feature matching degree is less than a preset variant feature matching degree;

[0010] Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected, and the two are stored in the first cluster set and the second cluster set respectively;

[0011] The reads in the variant region that are not stored in the first cluster set and the second cluster set are recorded as the to-be-assigned reads, and are compared with the reads of the first cluster set and the second cluster set, respectively. The reads of the first cluster set are recorded as the first representative reads, and the reads of the second cluster set are recorded as the second representative reads.

[0012] According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and according to the variation feature matching degree, the read segment to be assigned is stored in the first cluster set and the second cluster set;

[0013] S4: The variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to a partial order alignment algorithm or genome assembly to construct consensus sequences, which are recorded as the first consensus sequence and the second consensus sequence;

[0014] S5: After the first consensus sequence and the second consensus sequence are aligned with the local reference genome sequence of the variable region through the whole genome, a first alignment result and a second alignment result are obtained;

[0015] Extracting variation features based on the first comparison results and the second comparison results, performing pairwise comparisons of the extracted variation features with all variations in the first cluster set and the second cluster set to determine effective variations;

[0016] S6: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

[0017] In step S3, according to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and according to the variation feature matching degree, the read segment to be assigned is stored in the first cluster set and the second cluster set, specifically:

[0018] If the read to be assigned has an intersection with both the first representative read and the second representative read, then the variation feature matching degrees of the read to be assigned, the first representative read and the second representative read are calculated respectively to obtain the first variation feature matching degree and the second variation feature matching degree, and the first variation feature matching degree is compared with the second variation feature matching degree;

[0019] If the first variation feature matching degree is not equal to the second variation feature matching degree, the read segment to be allocated is stored in a cluster set of read segments with a large variation feature matching degree;

[0020] If the read to be assigned has only an intersection with the first representative read, the variation feature matching degree between the read to be assigned and the first representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the second cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the first cluster set;

[0021] If the read to be assigned has only an intersection with the second representative read, the variation feature matching degree between the read to be assigned and the second representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the first cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the second cluster set.

[0022] In step S5, based on the first comparison result and the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the first cluster set and the second cluster set in pairs to determine the effective variation, specifically:

[0023] Based on the first comparison result, extract the variation feature, compare the extracted variation feature with all the variations in the first cluster set in pairs, and record it as a match if the first matching condition is met; if the number of recorded matches is not less than the preset number of read segments, record the variation feature as a valid variation;

[0024] Based on the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the second clustering set. If the first matching condition is met, it is recorded as a match; if the number of recorded matches is not less than the preset number of read segments, the variation feature is recorded as a valid variation.

[0025] In step S6, based on the first consensus sequence and the second consensus sequence, allele determination is performed to determine the genotype of the effective variation, specifically:

[0026] If two valid variants, respectively from the first consensus sequence and the second consensus sequence, can overlap each other after being extended to both sides by a preset length at the variant position, the two valid variants are alleles, and the variant type is determined;

[0027] If both valid variants are insertion variants and the sequence consistency is not less than the sequence consistency threshold, the two valid variants are merged and the genotype is recorded as a homozygous variant;

[0028] If both effective variants are insertion variants and the sequence identity is less than the sequence identity threshold, the genotypes of the two effective variants are recorded as heterozygous variants;

[0029] If both valid variants are deletion variants and the variant length ratio is not less than the variant length ratio threshold, the two valid variants are merged and the genotype is recorded as a homozygous variant;

[0030] If both effective variants are deletion variants and the variant size ratio is less than the variant size ratio threshold, the genotypes of the two effective variants are recorded as heterozygous variants;

[0031] If the two effective variants are of different types, the genotypes of the two effective variants are recorded as heterozygous variants.

[0032] In step S3, the variation feature sequences of the two read segments are extracted and the variation feature matching degree is calculated, specifically:

[0033] Extracting the variation feature sequences of the two read segments, calculating the matching scores according to the second matching condition, and constructing the variation feature comparison matrix;

[0034] According to the variation feature alignment matrix, the variation feature sequences of the two read segments are aligned to obtain a matching vector, and the variation feature matching degree is calculated according to the matching vector.

[0035] The second matching condition is that the variation types are the same, the variation length ratio is not less than a preset variation length ratio, the reference distance is not greater than a reference distance threshold, and the sequence consistency of the inserted variation is not less than a preset sequence consistency.

[0036] The variation feature matching degree is calculated according to the matching vector, specifically,

[0037] According to the formula: , calculate the matching degree of variant features;

[0038] In the formula, is the matching degree of variant features; The value of is 1 or 0, indicating a match or a mismatch. ; is the length of the matching vector.

[0039] In step S5, the first matching condition is that the variation type is the same, the variation size is not less than the preset variation length, the variation length ratio is not less than the preset variation length ratio, the reference distance is not greater than the preset reference distance, and the sequence consistency of the inserted variation is not less than the preset sequence consistency.

[0040] In step S4, the variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to the partial alignment algorithm or the genome assembly, including smoothing the reads in the first cluster set and the second cluster set, specifically:

[0041] Based on the CIGAR data of the reads in the first cluster set and the second cluster set, if the CIGAR operand is "M", the variant feature sequence of the read is replaced with the corresponding sequence of the reference genome;

[0042] If the CIGAR operand is "I" or "D", the position is indicated as an insertion or deletion; for an insertion, the base sequence of the insertion variation is deleted from the read segment; for a deletion, the deleted base sequence is copied from the reference genome and inserted into the read segment.

[0043] The present invention also provides an insertion or deletion variation detection system based on allele perception, comprising:

[0044] Variation signal extraction module: used to compare sequencing data with the reference genome, and extract variation signals from the obtained comparison data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive shearing data, and excessive short insertions or deletions;

[0045] Mutation region calculation module: based on the mutation signal, calculate the interval between two adjacent normal regions to obtain the mutation region;

[0046] Clustering module: used to compare all reads in the variant region in pairs. If there is an intersection, the variant feature sequences of the two reads are extracted, and the variant feature matching degree is calculated. If the variant feature matching degree is less than the preset variant feature matching degree, the reads are stored in the first set;

[0047] Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected, and the two are stored in the first cluster set and the second cluster set respectively;

[0048] The reads in the variant region that are not stored in the first cluster set and the second cluster set are recorded as the to-be-assigned reads, and are compared with the reads of the first cluster set and the second cluster set, respectively. The reads of the first cluster set are recorded as the first representative reads, and the reads of the second cluster set are recorded as the second representative reads.

[0049] According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and according to the variation feature matching degree, the read segment to be assigned is stored in the first cluster set and the second cluster set;

[0050] Consensus sequence construction module: used for constructing consensus sequences of variant feature sequences of reads in the first cluster set and the second cluster set respectively through a partial order alignment algorithm or genome assembly, and recorded as the first consensus sequence and the second consensus sequence;

[0051] Effective variation determination module: used to compare the first consensus sequence and the second consensus sequence with the local reference genome sequence of the variation region through the whole genome, and obtain the first comparison result and the second comparison result;

[0052] Extracting variation features based on the first comparison results and the second comparison results, performing pairwise comparisons of the extracted variation features with all variations in the first cluster set and the second cluster set to determine effective variations;

[0053] Genotype determination module: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

[0054] Beneficial effects: The present invention extracts the insertion or deletion variation characteristics of the sequencing data at the variation position, uses a variety of matching index calculation methods, makes full use of the diploid variation characteristic information in the genome, designs a variation characteristic clustering method for different alleles, and designs an accurate insertion or deletion variation detection method based on the clustering results, thereby improving the variation detection in complex regions, especially in multi-alleles, thereby improving the variation detection accuracy. This detection method can be applied to complex genomic regions, including those with multiple alleles and medically related genes, and improves the detection capabilities in these regions. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] Figure 1 Schematic diagram of a flow chart of an insertion or deletion variation detection method based on allele perception in an embodiment. DETAILED DESCRIPTION

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

[0057] The present invention provides a method for detecting insertion or deletion variation based on allele perception, such as Figure 1 As shown, the following steps are included:

[0058] 1. Calculate genomic regions with insertion or deletion variants

[0059] S1: Sequencing data is aligned with the reference genome, and variation signals are extracted from the aligned data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, too many insertions or deletions, too much clipping data, and too many short insertions or deletions.

[0060] Specifically, each chromosome is detected one by one according to the order of chromosomes in the reference genome, and the genomic regions containing the variation 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 , OccC , Occ G , Occ T , Occ N ),in Occ X Indicates the base in the alignment data aligned to this base position X The number of X for A , C , G , T , N To reduce memory usage, each chromosome is divided into many 1Mb intervals, called blocks or detection blocks. For each block detection, the alignment data in the block is first loaded, and then a sliding window of 500bp is used to detect the block in the order of 5' to 3' of the reference genome, and the genomic regions containing insertion or deletion variations are calculated.

[0061] Extract the variation signal from the obtained comparison data, where:

[0062] ① Non-consistent bases: All reads covering a certain base in the reference sequence are counted. If the maximum number of bases of a type at the base position is different from the bases in the reference sequence, or the ratio of the maximum number of bases of a type at the base position to the total number of bases at the position is lower than the threshold (such as <0.8), the base at the base position is called a non-consistent base. Non-consistent bases will appear in most variant positions, so the region containing this variant signal usually corresponds to the region of insertion or deletion variants.

[0063] ② Abnormal base coverage depth: represents high (or low) base coverage depth, which means that the ratio of the coverage depth within the 500bp detection window to the average coverage depth of the entire block (1Mb) is >2 (or <0.5). This signal is mainly used to detect structural variations of insertion, deletion and duplication types.

[0064] ③Zero coverage depth: refers to the lack of reads data aligned to a certain base position in the genome. Zero coverage depth is a special type of abnormal base coverage depth. At certain insertion and deletion variation positions, zero coverage depth coverage loss will occur.

[0065] ④ Excessive insertions or deletions: refers to excessive insertions (deletions) at a certain base position of the reference genome, that is, the number of insertions (deletions) at this base position is greater than the preset insertion (deletion) number threshold. This variation signal will appear at most insertion and deletion variation positions, and other types of variation will also usually have excessive insertions and deletions.

[0066] ⑤ Excessive clipping data: refers to the number of clipping data (terminal length greater than 100bp) contained at a certain base position of the reference genome is greater than 3 times the average number of clippings per base in the entire detection block (1Mb). This variation signal is mainly used to detect long sequence structural variation. Long sequence structural variation usually leads to a large amount of clipping data at the variation boundary.

[0067] ⑥ Too many short insertions or deletions: refers to too many short insertions or deletions at a certain base position of the reference genome, that is, the sum of the number of short insertions or deletions at the base position is greater than twice the average number of short insertions or deletions per base in the detection block (1Mb). This variation signal is mainly used to detect variations whose length is less than the insertion or deletion length threshold. Such variations will be mistakenly filtered out as sequencing noise in the previous noise filtering step due to their small length.

[0068] The present invention describes the variation characteristics at the base level, which can more accurately extract the variation signal and facilitate more accurate subsequent allele perception clustering results.

[0069] S2: Based on the variation signal, the interval between two adjacent normal regions is calculated to obtain the variation region.

[0070] On the basis of extracting the variation signal of each base position in the genome, the minimum distance d between two adjacent variations is set (for example, 10bp), and the genomic position containing the variation is indirectly calculated by detecting two adjacent non-variant genomic regions (normal regions) to obtain the variation region. Therefore, the normal region is defined as: the region size is at least d (i.e., the minimum distance between two adjacent variations), and each base position in the region does not contain variation signals: non-consistent bases, abnormal base coverage depth, zero coverage depth, too many insertions or deletions, too much clipping data, and too many short insertions or deletions.

[0071] Set the minimum variant length for insertion or deletion variant detection, denoted as minSVsize. Only variants with a length greater than this threshold will be detected.

[0072] Afterwards, the specific steps to obtain the variant region are as follows: within each 500bp detection window, check two adjacent normal regions, calculate the interval size between the two adjacent normal regions, and the region where the interval is located is the genomic region containing insertion or deletion mutations, referred to as the variant region.

[0073] While calculating the variation region, the size of the variation in the region can also be recorded. Specifically:

[0074] Record the average length of the largest variation in the region and store it; if the region has both insertion and deletion variation (that is, allele variation), take the maximum value of the average length of the two variations and store it. Because the stored variation size is an average value, when clustering is performed later, the expansion distance when extracting variation feature sequences can be dynamically adjusted according to the stored average length, saving resources and time consumed by calculation.

[0075] 2. After obtaining the variant region, the next step is to cluster the variant features of allele perception in the insertion or deletion region:

[0076] S3: All reads in the variant region are compared pairwise. If there is an intersection, the variant feature sequences of the two reads are extracted and the variant feature matching degree is calculated. If the variant feature matching degree is less than the preset variant feature matching degree, it can be expressed as , the read segment is stored in the first set.

[0077] Regarding the variant feature sequence, there are the following definitions:

[0078] The mutation feature is regarded as an element in the sequence, and the mutation feature is defined as a six-tuple ,in, is the name of the chromosome where the mutation occurs, and Chromosome The starting and ending positions on is the structural variation type, is the length of the structural variation, The variant base information is the sequence base information of the variant. The variant features in the 1kb region before and after the variant of the two read segments are extracted respectively, and the variant feature sequence is obtained by sorting the variant features according to the genome coordinates. .

[0079] Preferably, in step S3, the variation feature sequences of the two read segments are extracted and the variation feature matching degree is calculated, specifically:

[0080] Extracting the variation feature sequences of the two read segments, calculating the matching scores according to the second matching condition, and constructing the variation feature comparison matrix;

[0081] According to the variation feature comparison matrix, the variation feature sequences of the two read segments are compared to obtain a matching vector, where the matching vector , represents the result of each pair of variant feature comparison, there are 1 to n comparison results in total, The value of is 1 or 0, indicating a match or no match. Represents the matching vector According to the matching vector, the matching degree of the mutation feature is calculated.

[0082] Among them, the second matching condition is that the variation type is the same, the variation length ratio is not less than a preset variation length ratio (such as 0.7), the reference distance is not greater than a reference distance threshold (such as 100bp), and the sequence consistency of the inserted variation is not less than a preset sequence consistency (such as 0.7).

[0083] For the detection of insertion and deletion mutations, the present invention adopts more specific matching conditions, especially for insertion mutations, the sequence consistency is calculated, and the corresponding reference sequence is extracted according to the mutation position between the upstream starting point and the downstream end point of the two insertion mutations. This step can be dynamically adjusted according to the average length of the maximum variation stored to expand a certain interval. The variation feature sequence is then merged with the reference sequence to create two complete sequences to be compared. Subsequently, the two sequences to be compared are globally aligned, and the optimal matching sequence is generated by the Needleman-Wunsch (NW) two-sequence global alignment algorithm. Finally, the sequence consistency is calculated and reported using the optimal matching sequence.

[0084] In addition, the variation feature matching degree is calculated based on the matching vector, specifically,

[0085] According to the formula: , calculate the matching degree of variant features;

[0086] In the formula, is the matching degree of variant features; The value of is 1 or 0, indicating a match or no match; is the length of the matching vector. The value range is 0 to 1. In particular, when the length of the variant feature sequence is 0, it means that both reads do not contain variant features, and filtering is performed.

[0087] Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected and stored in the first cluster set and the second cluster set, respectively.

[0088] That is, when performing allele-aware variant feature clustering, first the read set within the variant region is , extract each read segment The variant feature sequences in the variant region are then R Perform allele-aware variant feature clustering. The specific steps are as follows:

[0089] The first step is to initialize the cluster read set. The goal is to find two reads with variant characteristics and differences as the initial cluster set. A and B , the specific method is as follows:

[0090] Step 1: Assemble reads within the variant region Sort by read length in descending order.

[0091] Step 2: For any read in the variant region , loop execution: for any read segment , loop execution: if and If there is a span intersection, the mutation feature matching degree is calculated based on the mutation feature sequences of the two. ;like ,Will exist Subscript in , The value is stored in the first collection ; After the inner loop ends, if If it is not empty, the outer loop ends.

[0092] Step 3: If Not empty, select the same The read segment with the largest number is added to the first cluster set A , for example join in A， Then select the same The second largest number of reads is added to the second cluster set B , for example, another join in B ;

[0093] Step 4: If Empty, All reads in A ; End the clustered read set initialization algorithm.

[0094] The existing method initializes the first clustering set A When , the first read segment is selected as the initial read segment by default, and the read segment is excluded from the subsequent initialization process, which may affect the second clustering set B The present invention is based on The clustering results are used to select the clustering set A and B The initial read fragments avoid the first clustering set A The reads selected during initialization are not representative of their regions, which leads to errors and mistakes in subsequent clustering.

[0095] Initialize clustering set A , BAfter the data is read, before clustering, the data will be loaded first. This method does not need to load all the reads of a region, but performs downsampling operations through two rounds of data retrieval while reading the data according to the coverage. The specific method is: before reading the data, the total number and total length of the variation in the region are calculated. If the current coverage is greater than the maximum coverage threshold, the first random downsampling is performed, otherwise the data of all the reads in the region are loaded; then before clustering, if the current coverage is greater than the preset coverage, the second random downsampling is performed. This method greatly reduces the resource consumption when reading data, especially in some complex areas with high coverage, the effect is more obvious.

[0096] Reads in the variant region that are not stored in the first cluster set or the second cluster set are recorded as reads to be assigned , respectively, are compared with the reads of the first cluster set and the second cluster set, where the reads of the first cluster set are recorded as the first representative reads a , the reads of the second cluster set are recorded as the second representative reads b .

[0097] According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and the read segment to be assigned is stored in the first cluster set and the second cluster set according to the variation feature matching degree, specifically:

[0098] If the read to be assigned has an intersection with both the first representative read and the second representative read, then the variation feature matching degrees of the read to be assigned, the first representative read and the second representative read are calculated respectively to obtain the first variation feature matching degree and the second variation feature matching degree, and the first variation feature matching degree is compared with the second variation feature matching degree;

[0099] If the first variation feature matching degree is not equal to the second variation feature matching degree, the read segment to be allocated is stored in a cluster set of read segments with a large variation feature matching degree;

[0100] If the read to be assigned has an intersection with the first representative read, the variation feature matching degree between the read to be assigned and the first representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the second cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the first cluster set;

[0101] If the read to be assigned has an intersection with the second representative read, the variation feature matching degree between the read to be assigned and the second representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the first cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the second cluster set.

[0102] That is, initialize the clustering setA , B Afterwards, the read segment collection R Perform allele-aware variant feature clustering as follows:

[0103] Step 1: For Not included in A or B Any read segment of , loop execution: find the corresponding of A , B The first representative read segment a , the second representative read segment b ;like and a, b All have span intersections, so calculate and a, b The matching degree of mutation characteristics , , according to the size relationship: when or When join in A , B ;when , then temporarily join C gather;

[0104] like Only with a or b There is a span intersection, calculate the representative sum of the span intersection The matching degree of mutation characteristics , Join another representative's set, and vice versa join a representative's set.

[0105] Step 2: For C Any read in , loop execution: find the corresponding of A , B The read segment represents a , b ;calculate and a , b The matching degree of mutation characteristics , ; According to the size relationship: when or When join in A , B ; End the algorithm.

[0106] The present invention uses a variety of matching index calculation methods to comprehensively analyze the variation characteristic signals and clustering methods of suspicious areas, fully utilizes the variation characteristic signals, makes the detection results more comprehensive and reliable, improves the accuracy of variation detection, and reduces false detections and missed detections, especially in terms of allele variation.

[0107] 3. After clustering, the next step is to detect heterozygous insertion or deletion mutations:

[0108] S4: The variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to the partial order alignment algorithm (abPOA) or genome assembly (wtdbg2) to construct consensus sequences, which are recorded as the first consensus sequence and the second consensus sequence.

[0109] The partial order alignment algorithm (abPOA) alignment process consumes a lot of memory, and it is impossible to estimate the memory consumption in advance. Before the abPOA work starts, the method of the present invention will decide whether to perform the partial order alignment based on the memory usage of the machine. If the used memory exceeds 80% of the current machine memory, the next abPOA work will be made to wait until the abPOA work is completed and the memory is released. This can avoid program interruption due to exceeding the memory limit.

[0110] However, since abPOA may have deviations when constructing consensus sequences for some special intervals, some variants may be missed or misdetected. For these reads, the variant interval of the read is re-detected using the genome assembly (wtdbg2). The specific method is: re-read the read set of the clustering results in the second step, including the description information and base sequence of each read, and then assemble the genome (wtdbg2) to obtain the consensus sequence.

[0111] Furthermore, in step S4, before the variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to the partial order alignment algorithm (abPOA) or the genome assembly (wtdbg2), the reads in the first cluster set and the second cluster set are smoothed to optimize the variants and smooth out the variants that have little effect on the detection results, specifically:

[0112] Based on the CIGAR data of the reads in the first cluster set and the second cluster set, if the CIGAR operand is "M", the variant feature sequence of the read is replaced with the corresponding sequence of the reference genome;

[0113] If the CIGAR operand is "I" or "D", the position is indicated as an insertion or deletion; for an insertion, the base sequence of the insertion variation is deleted from the read segment; for a deletion, the deleted base sequence is copied from the reference genome and inserted into the read segment.

[0114] S5: The first consensus sequence and the second consensus sequence are respectively aligned with the local reference genome sequence of the variant region through the whole genome alignment (minimap2) to realign the consensus sequence to the local reference sequence of the variant region, and obtain the first alignment result and the second alignment result, which include the location, size and type of the insertion or deletion variation in the variant region. In this alignment result, the heterozygous insertion or deletion variation usually produces two consensus sequences carrying the variation information.

[0115] The variant features extracted from the comparison results are then further filtered, including the following operations:

[0116] Based on the first comparison result and the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the first cluster set and the second cluster set in pairs to determine the effective variation, specifically:

[0117] Based on the first comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the first clustering set. If the first matching condition is met, it is recorded as a match; if the number of recorded matches is not less than the preset number of read segments, the variation feature is recorded as a valid variation; otherwise, it is considered that the variation feature has failed the variation detection and the data record is directly discarded.

[0118] Based on the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the second clustering set. If the first matching condition is met, it is recorded as a match; if the number of recorded matches is not less than the preset number of read segments, the variation feature is recorded as a valid variation; otherwise, it is considered that the variation feature has failed the variation detection and the data record is directly discarded.

[0119] Preferably, in step S5, the first matching condition is that the variation type is the same, the variation size is not less than a preset variation length (such as a minimum variation length minSVsize), the variation length ratio is not less than a preset variation length ratio (such as 0.7), the reference distance is not greater than a preset reference distance (such as 200bp), and the sequence consistency of the inserted variation is not less than a preset sequence consistency (such as 0.7).

[0120] S6: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

[0121] For the detected insertion or deletion variation, its genotype needs to be further detected. This process is called genotyping, also known as genotype calculation. By detecting the genotype, complex heterozygous variation can be detected more specifically to distinguish whether it is a homozygous variation or a heterozygous variation.

[0122] Specifically:

[0123] If two valid variants, respectively from the first consensus sequence and the second consensus sequence, can overlap each other after being extended to both sides by a preset length at the variant position, the two valid variants are alleles, and the variant type is determined;

[0124] If both valid variants are insertion variants and the sequence consistency is not less than the sequence consistency threshold (0.9), the two valid variants are merged and the genotype is recorded as a homozygous variant;

[0125] If both effective variants are insertion variants and the sequence identity is less than the sequence identity threshold, the genotypes of the two effective variants are recorded as heterozygous variants;

[0126] If both valid variants are deletion variants and the variant length ratio is not less than the variant length ratio threshold (0.9), the two valid variants are merged and the genotype is recorded as a homozygous variant;

[0127] If both effective variants are deletion variants and the variant size ratio is less than the variant size ratio threshold, the genotypes of the two effective variants are recorded as heterozygous variants;

[0128] If the two effective variants are of different types, the genotypes of the two effective variants are recorded as heterozygous variants.

[0129] Furthermore, the above operations can be expressed as:

[0130] When there are two consensus sequences in the variable region, before detecting the genotype, it is necessary to determine whether the two valid variants in the region belong to two alleles of the same locus. The determination method is as follows:

[0131] If the following conditions are met at the same time: a) the two effective variants are located on different consensus sequences, i.e., they correspond to different haplotypes; and b) they can overlap with each other after extending 5 bp to the left and right of the two effective variant positions, then the two effective variants are judged to be two alleles.

[0132] For genotyping of two allele variants, there are four situations:

[0133] 1. If two alleles are variants, and further satisfy the following conditions c) and d): c) the types of the two variants are the same; d) the sequence consistency of the two insertion variants is ≥ the sequence consistency threshold, such as 0.9, or the variation length ratio of the two deletion variants is ≥ the variation length ratio threshold, such as 0.9; then the two alleles are homozygous variants and are merged, and the genotype is marked as 1 / 1;

[0134] 2. Both alleles are homozygous variants and are merged, with the marker genotype being 1 / 1;

[0135] 3. If the two alleles have different mutation types, the two alleles are heterozygous mutations and the marker genotype is 0 / 1;

[0136] 4. If the two alleles have the same variation type, the insertion variation sequence consistency is less than the sequence consistency threshold, such as 0.9, or the deletion variation length ratio is less than the variation length ratio threshold, such as 0.9, then the two alleles are heterozygous variations, and the marker genotype is 0 / 1;

[0137] If the variant exists in only one consensus sequence and its allele is consistent with the reference sequence, then the variant is a heterozygous variant and the genotype is 0 / 1.

[0138] In addition, if there is only one consensus sequence in the variable region, the genotypes 1 / 1, 0 / 1, and 0 / 0 are determined based on the allele frequencies.

[0139] The operations can be as follows:

[0140] First, calculate the variation Allele frequency :

[0141] ,

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

[0143] Then, based on the allele frequency, the variation is calculated using the following formula Genotype:

[0144] ,

[0145] Among them, the genotype "1 / 1" indicates that when the allele frequency is not less than 0.75, it is recorded as 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.

[0146] The present invention extracts the insertion or deletion variation characteristics of the sequencing data at the variation position, uses a variety of matching index calculation methods, makes full use of the diploid variation characteristic information in the genome, designs a variation characteristic clustering method for different alleles, and designs an accurate insertion or deletion variation detection method based on the clustering results, thereby improving the variation detection in complex regions, especially in multi-alleles, thereby improving the variation detection accuracy. The detection method can be applied to complex genomic regions, including those with multiple alleles and medically related genes, and improves the detection capabilities in these regions.

[0147] The present invention also provides an insertion or deletion variation detection system based on allele perception, comprising:

[0148] Variation signal extraction module: used to compare sequencing data with the reference genome, and extract variation signals from the obtained comparison data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive shearing data, and excessive short insertions or deletions;

[0149] Mutation region calculation module: based on the mutation signal, calculate the interval between two adjacent normal regions to obtain the mutation region;

[0150] Clustering module: used to compare all reads in the variant region in pairs. If there is an intersection, the variant feature sequences of the two reads are extracted, and the variant feature matching degree is calculated. If the variant feature matching degree is less than the preset variant feature matching degree, the reads are stored in the first set;

[0151] Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected, and the two are stored in the first cluster set and the second cluster set respectively;

[0152] The reads in the variant region that are not stored in the first cluster set and the second cluster set are recorded as the to-be-assigned reads, and are compared with the reads of the first cluster set and the second cluster set, respectively. The reads of the first cluster set are recorded as the first representative reads, and the reads of the second cluster set are recorded as the second representative reads.

[0153] According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and according to the variation feature matching degree, the read segment to be assigned is stored in the first cluster set and the second cluster set;

[0154] Consensus sequence construction module: used for constructing consensus sequences of variant feature sequences of reads in the first cluster set and the second cluster set respectively through a partial order alignment algorithm or genome assembly, and recorded as the first consensus sequence and the second consensus sequence;

[0155] Effective variation determination module: used to compare the first consensus sequence and the second consensus sequence with the local reference genome sequence of the variation region through the whole genome, and obtain the first comparison result and the second comparison result;

[0156] Extracting variation features based on the first comparison results and the second comparison results, performing pairwise comparisons of the extracted variation features with all variations in the first cluster set and the second cluster set to determine effective variations;

[0157] Genotype determination module: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

Claims

1. A method for detecting insertion or deletion variation based on allele perception, characterized in that: The following steps are involved: S1: Sequencing data is aligned with the reference genome, and variation signals are extracted from the aligned data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive shearing data, and excessive short insertions or deletions; S2: Based on the variation signal, the interval between two adjacent normal regions is calculated to obtain the variation region; S3: performing pairwise comparisons on all reads of the variant region. If there is an intersection, extracting the variant feature sequences of the two reads, calculating the variant feature matching degree, and storing the reads in the first set if the variant feature matching degree is less than a preset variant feature matching degree; Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected, and the two are stored in the first cluster set and the second cluster set respectively; The reads in the variant region that are not stored in the first cluster set and the second cluster set are recorded as the to-be-assigned reads, and are compared with the reads of the first cluster set and the second cluster set, respectively. The reads of the first cluster set are recorded as the first representative reads, and the reads of the second cluster set are recorded as the second representative reads. According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and the read segment to be assigned is stored in the first cluster set and the second cluster set according to the variation feature matching degree, specifically: If the read to be assigned has an intersection with both the first representative read and the second representative read, then the variation feature matching degrees of the read to be assigned, the first representative read and the second representative read are calculated respectively to obtain the first variation feature matching degree and the second variation feature matching degree, and the first variation feature matching degree is compared with the second variation feature matching degree; If the first variation feature matching degree is not equal to the second variation feature matching degree, the read segment to be allocated is stored in a cluster set of read segments with a large variation feature matching degree; If the read to be assigned has only an intersection with the first representative read, the variation feature matching degree between the read to be assigned and the first representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the second cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the first cluster set; If the read to be assigned has only an intersection with the second representative read, the variation feature matching degree between the read to be assigned and the second representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the first cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the second cluster set; S4: The variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to a partial order alignment algorithm or genome assembly to construct consensus sequences, which are recorded as the first consensus sequence and the second consensus sequence; S5: After the first consensus sequence and the second consensus sequence are aligned with the local reference genome sequence of the variable region through the whole genome, a first alignment result and a second alignment result are obtained; Based on the first comparison result and the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the first cluster set and the second cluster set in pairs to determine the effective variation, specifically: Based on the first comparison result, extract the variation feature, compare the extracted variation feature with all the variations in the first cluster set in pairs, and record it as a match if the first matching condition is met; if the number of recorded matches is not less than the preset number of read segments, record the variation feature as a valid variation; Based on the second comparison result, the variation feature is extracted, and the extracted variation feature is compared with all the variations in the second cluster set in pairs. If the first matching condition is met, it is recorded as a match; if the number of recorded matches is not less than the preset number of read segments, the variation feature is recorded as a valid variation; S6: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

2. The method for detecting insertion or deletion variation based on allele perception according to claim 1, characterized in that: In step S6, based on the first consensus sequence and the second consensus sequence, allele determination is performed to determine the genotype of the effective variation, specifically: If two valid variants, respectively from the first consensus sequence and the second consensus sequence, can overlap each other after being extended to both sides by a preset length at the variant position, the two valid variants are alleles, and the variant type is determined; If both valid variants are insertion variants and the sequence consistency is not less than the sequence consistency threshold, the two valid variants are merged and the genotype is recorded as a homozygous variant; If both effective variants are insertion variants and the sequence identity is less than the sequence identity threshold, the genotypes of the two effective variants are recorded as heterozygous variants; If both valid variants are deletion variants and the variant length ratio is not less than the variant length ratio threshold, the two valid variants are merged and the genotype is recorded as a homozygous variant; If both effective variants are deletion variants and the variant size ratio is less than the variant size ratio threshold, the genotypes of the two effective variants are recorded as heterozygous variants; If the two effective variants are of different types, the genotypes of the two effective variants are recorded as heterozygous variants.

3. The method for detecting insertion or deletion variation based on allele perception according to claim 1, characterized in that: In step S3, the variation feature sequences of the two read segments are extracted and the variation feature matching degree is calculated, specifically: Extracting the variation feature sequences of the two read segments, calculating the matching scores according to the second matching condition, and constructing the variation feature comparison matrix; According to the variation feature alignment matrix, the variation feature sequences of the two read segments are aligned to obtain a matching vector, and the variation feature matching degree is calculated according to the matching vector.

4. The method for detecting insertion or deletion variation based on allele perception according to claim 3, characterized in that: The second matching condition is that the variation types are the same, the variation length ratio is not less than a preset variation length ratio, the reference distance is not greater than a reference distance threshold, and the sequence consistency of the inserted variation is not less than a preset sequence consistency.

5. The method for detecting insertion or deletion variation based on allele perception according to claim 3, characterized in that: The calculation of the variation feature matching degree according to the matching vector is specifically as follows: According to the formula: , calculate the matching degree of variant features; In the formula, is the matching degree of variant features; The value of is 1 or 0, indicating a match or a mismatch. ; is the length of the matching vector.

6. The method for detecting insertion or deletion variation based on allele perception according to claim 1, characterized in that: In step S5, the first matching condition is that the variation type is the same, the variation size is not less than the preset variation length, the variation length ratio is not less than the preset variation length ratio, the reference distance is not greater than the preset reference distance, and the sequence consistency of the inserted variation is not less than the preset sequence consistency.

7. The method for detecting insertion or deletion variation based on allele perception according to claim 1, characterized in that: In step S4, the variant feature sequences of the reads in the first cluster set and the second cluster set are respectively subjected to the partial alignment algorithm or the genome assembly, including smoothing the reads in the first cluster set and the second cluster set, specifically: Based on the CIGAR data of the reads in the first cluster set and the second cluster set, if the CIGAR operand is "M", the variant feature sequence of the read is replaced with the corresponding sequence of the reference genome; If the CIGAR operand is "I" or "D", the position is indicated as an insertion or deletion; for an insertion, the base sequence of the insertion variation is deleted from the read segment; for a deletion, the deleted base sequence is copied from the reference genome and inserted into the read segment.

8. An insertion or deletion variation detection system based on allele perception, characterized in that: include: Variation signal extraction module: used to compare sequencing data with the reference genome, and extract variation signals from the obtained comparison data, including: non-consistent bases, abnormal base coverage depth, zero coverage depth, excessive insertions or deletions, excessive shearing data, and excessive short insertions or deletions; Mutation region calculation module: based on the mutation signal, calculate the interval between two adjacent normal regions to obtain the mutation region; Clustering module: used to compare all reads in the variant region in pairs. If there is an intersection, the variant feature sequences of the two reads are extracted, and the variant feature matching degree is calculated. If the variant feature matching degree is less than the preset variant feature matching degree, the reads are stored in the first set; Based on the first set, one of the reads with the same variation feature match and the largest number and one of the reads with the same variation feature match and the second largest number are selected, and the two are stored in the first cluster set and the second cluster set respectively; The reads in the variant region that are not stored in the first cluster set and the second cluster set are recorded as the to-be-assigned reads, and are compared with the reads of the first cluster set and the second cluster set, respectively. The reads of the first cluster set are recorded as the first representative reads, and the reads of the second cluster set are recorded as the second representative reads. According to the intersection between the read segment to be assigned and the first representative read segment and the second representative read segment, the variation feature matching degree is calculated, and the read segment to be assigned is stored in the first cluster set and the second cluster set according to the variation feature matching degree, specifically: If the read to be assigned has an intersection with both the first representative read and the second representative read, then the variation feature matching degrees of the read to be assigned, the first representative read and the second representative read are calculated respectively to obtain the first variation feature matching degree and the second variation feature matching degree, and the first variation feature matching degree is compared with the second variation feature matching degree; If the first variation feature matching degree is not equal to the second variation feature matching degree, the read segment to be allocated is stored in a cluster set of read segments with a large variation feature matching degree; If the read to be assigned has only an intersection with the first representative read, the variation feature matching degree between the read to be assigned and the first representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the second cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the first cluster set; If the read to be assigned has only an intersection with the second representative read, the variation feature matching degree between the read to be assigned and the second representative read is calculated. If the obtained variation feature matching degree is less than 1, the read to be assigned is stored in the first cluster set; if the obtained variation feature matching degree is not less than 1, the read to be assigned is stored in the second cluster set; Consensus sequence construction module: used for constructing consensus sequences of variant feature sequences of reads in the first cluster set and the second cluster set respectively through a partial order alignment algorithm or genome assembly, and recorded as the first consensus sequence and the second consensus sequence; Effective variation determination module: used to compare the first consensus sequence and the second consensus sequence with the local reference genome sequence of the variation region through the whole genome, and obtain the first comparison result and the second comparison result; Based on the first comparison result and the second comparison result, the variation features are extracted, and the extracted variation features are compared with all the variations in the first cluster set and the second cluster set in pairs to determine the effective variation, specifically: Based on the first comparison result, extract the variation feature, compare the extracted variation feature with all the variations in the first cluster set in pairs, and record it as a match if the first matching condition is met; if the number of recorded matches is not less than the preset number of read segments, record the variation feature as a valid variation; Based on the second comparison result, the variation feature is extracted, and the extracted variation feature is compared with all the variations in the second cluster set in pairs. If the first matching condition is met, it is recorded as a match; if the number of recorded matches is not less than the preset number of read segments, the variation feature is recorded as a valid variation; Genotype determination module: Based on the first consensus sequence and the second consensus sequence, allele judgment is performed to determine the genotype of the effective variation.

Images