Gene sequence variation detection method based on robust statistics and multi-strategy fusion

Abstract

The invention relates to the technical field of gene variation detection, in particular to a gene sequence variation detection method based on robust statistics and multi-strategy fusion, which comprises the following steps: comparing to generate a sorted BAM file; preprocessing a missing value and an N position in the data, replacing the missing value by adopting a 0 filling strategy, and removing the data at the N position; extracting an RD signal and an MQ signal as characteristic values; performing GC deviation correction on the extracted RD signal; performing smooth noise reduction processing on the RD signal and the MQ signal; a two-step progressive segmentation strategy is adopted, and continuous fragments with highly consistent RD values are identified; constructing two-dimensional profiles of the RD signal and the MQ signal, and performing standardization processing; carrying out serial repeated detection and calculating an abnormal score; setting a threshold value by using a box plot method, and judging an abnormal value; and refining a tandem repeated region. According to the method, the problem of tandem repeat variation detection under low-coverage sequencing data is effectively solved, and the detection accuracy and reliability are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of gene variation detection in bioinformatics, and in particular to a gene sequence variation detection method based on robust statistics and multi-strategy fusion. Background Art

[0002] In the field of genetics, a large number of studies have shown that there are rich and diverse genetic variations in the human genome, and these variations mainly originate from gene mutations. According to the difference in the number of mutated bases, it is divided into single nucleotide polymorphisms and structural variations. Among them, structural variation is one of the main types of genetic variation. It is usually manifested as sequence changes of more than 50bp in the biological genome, covering large fragment deletions, insertions, duplications, inversions, translocations and copy number variations of genomic sequences. In the genome of the normal population, structural variation accounts for about 13%, and at least 240 genes are homozygous deletion polymorphisms. Because of its rich types and high proportion, it plays a key role in explaining the diversity characteristics of organisms and populations.

[0003] Tandem repeats, as a specific form of structural variation, refer to genomic regions composed of adjacent repeating units, that is, a DNA sequence that is repeated continuously and closely in the genome. There are more than 1 million different tandem repeat sequences in the human genome, which are widely distributed within genes and in gene intergenic regions. This distribution characteristic enables them to have a profound impact on the structure and function of DNA, RNA, and proteins. Studies have confirmed that tandem repeats are closely related to a variety of diseases, such as inducing amplification diseases and causing gene silencing. Different cancers and acute myeloid leukemia have related characteristics. Therefore, it is particularly necessary to accurately detect tandem repeats in genetic testing, and the continuous development and gradual improvement of second-generation sequencing technology provide a solid data foundation for the detection of tandem repeats.

[0004] With the development of second-generation sequencing technology, more and more detection strategies based on second-generation sequencing data have been studied, which are mainly divided into four categories: detection methods based on double-end mapping, detection methods based on split reads, detection methods based on read depth, and detection methods based on re-joining. Detection methods based on double-end mapping and detection methods based on split reads are only applicable to double-end sequencing data. The detection methods based on double-end mapping have low resolution in detecting variant sites, while the detection methods based on split reads can accurately detect breakpoints but have relatively high requirements for data quality and computing resources. The detection method based on read depth is the mainstream algorithm for detecting tandem repeat variants. It is suitable for single-end and double-end sequencing data. It determines the duplication or deletion region based on the mapping depth and read distribution, but cannot accurately detect the boundary sites of the variant region.

[0005] The second generation sequencing technology has become an important tool for genomic research due to its low sequencing cost and high sequencing throughput, and has promoted the development of a large number of algorithms for detecting tandem repeat variations. However, existing detection methods such as BreakDancer, Pindel-TD, SVIM and ScanITD, although they perform well under specific conditions, all have certain limitations.

[0006] BreakDancer identifies structural variations by analyzing the abnormal mapping positions of paired ends in high-throughput sequencing data, but its detection accuracy may be reduced in low-coverage sequencing data. Pindel-TD adapts to the characteristics of tandem repeat short read alignment and can obtain single nucleotide resolution at breakpoints, but in low-coverage cases, it is difficult to obtain enough information to accurately identify variant breakpoints, which may lead to false positive or false negative results. SVIM uses a graph-based method and a new distance metric to cluster the detected structural variation features, which has advantages in scenarios requiring high-precision structural variation detection, but low-coverage samples may lead to insufficient feature collection or poor clustering effects, affecting detection performance. ScanITD uses a realignment strategy for split reads to accurately identify internal tandem repeats, but may be limited in detecting special types of tandem repeats, and detection accuracy may be affected in low-coverage sequencing data.

[0007] In summary, most existing methods for detecting tandem repeats can only show ideal results in a high-coverage sequencing data environment. In the case of low coverage, due to the high sensitivity to noise, it is easily interfered by noise factors such as sequencing errors, mapping errors and GC content deviations, resulting in an increase in false positive rates and a significant decrease in detection performance. Therefore, the detection of tandem repeats in low coverage conditions still faces many challenges, and there is an urgent need to develop more robust and efficient detection methods. Summary of the invention

[0008] The present invention discloses a gene sequence variation detection method based on robust statistics and multi-strategy fusion, which detects tandem repeat sequences in the human genome, effectively solves the problem of tandem repeat variation detection under low-coverage sequencing data, significantly improves the accuracy and reliability of detection, reduces the false positive rate, and provides a more accurate tool for genomics research.

[0009] The present invention provides a gene sequence variation detection method based on robust statistics and multi-strategy fusion, characterized in that it comprises the following steps: Step 1: Input the reference genome and sequencing samples, and generate a sorted BAM file through alignment; Step 2: Preprocess the missing values ​​and N positions in the data, use a 0-filling strategy to replace the missing values, and remove the N position data, where the N position refers to the position in the sequence where the specific base cannot be determined due to sequencing technology limitations or data quality issues; Step 3, aligning the genome in the sorted BAM file with the reference genome, extracting RD signal and MQ signal as feature values, wherein the RD signal reflects the reading depth, the MQ signal reflects the mapping quality, the RD signal corresponds to the numerical expression of the RD value, and the MQ signal corresponds to the numerical expression of the MQ value; Step 4, performing GC bias correction on the extracted RD signal; Step 5, smoothing and noise reduction processing is performed on the RD signal and the MQ signal; Step 6: A two-step progressive segmentation strategy is adopted to identify continuous segments with highly consistent RD values; Step 7, construct the two-dimensional profiles of the RD signal and the MQ signal, and perform standardization processing; Step 8, applying the minimum covariance determinant method to tandem duplication detection, and using a binary combination weighted evaluation method to calculate the anomaly score; Step nine, use the box plot method to set a reasonable threshold and determine the outliers; Step 10: Implement the SR-PEM collaborative optimization method to refine the tandem repeat region.

[0010] Furthermore, in step one, the BWA-MEM tool in the BWA software is used to compare and analyze the reference genome and the sequenced samples to generate an initial BAM file, and the SAMtools software is used to arrange the initial BAM file in order according to the reference genome position, and to screen out paired reads and split reads to obtain a sorted BAM file containing specific reads.

[0011] Furthermore, in step 2, data preprocessing operations are performed on the sorted BAM file to identify missing values ​​in the data and fill them with 0 uniformly; the operation of removing the N-position data is to define the window containing the N positions as an abnormal window, and set the corresponding RD value to a specific negative number for distinction; when calculating the RD value of each window, if the RD value is less than 0, the window is determined to be an abnormal window containing N positions and is removed, otherwise it is determined to be a normal window, and the corresponding RD value as well as the starting and ending positions of the window are recorded.

[0012] Furthermore, in step 3, samples are extracted from the BAM file and the read counts at each position are calculated, denoted as RC, and the RD value of each position in the genome in the sorted BAM file is extracted, i.e.,

[0013] in, Indicates RD value of a window; Indicates In the window RC value of each position; Indicates the length of the genomic window in the sequenced BAM file; Extract the MQ value for each position in the genome in the sorted BAM file, i.e.,

[0014] in, Indicates MQ value of a window; Indicates In the window The Mq value of each position; Indicates the length of the genomic window in the sequenced BAM file.

[0015] Furthermore, in step 4, the extracted RD signal is corrected for GC deviation, that is,

[0016] in, Indicates the corrected RD value of a window; represents the mean value of RD signal in all windows; Indicates The RD mean of the windows with highly consistent GC content; Indicates The RD value of the window.

[0017] Furthermore, in step five, the total variation model is used to smooth and reduce the noise of the RD signal and the MQ signal. The operation process includes calculating the gradient between adjacent points in the RD signal and the MQ signal to quantify the changes of each signal; adjusting the gradient value to the minimum to reduce the noise in the signal; in the optimization process, implementing adaptive weight allocation for the RD signal and the MQ signal, and finely adjusting the weight of each point, so as to effectively retain the structural characteristics of the signal.

[0018] Furthermore, in step 6, the two-step progressive segmentation strategy includes adopting a sliding window segmentation method to segment the sorted BAM files using a sliding window of fixed width. The genome in the sorted BAM file is locally processed, and the genome in the sorted BAM file is split into continuous non-overlapping windows of the same length; the genome in the sorted BAM file in windows is split into multiple genome fragments according to the amplitude of the RD signal using the cyclic binary segmentation method, and the RD mean of each fragment is calculated; the RD means of adjacent fragments are statistically tested for significance, and if a significant difference is detected, the segmentation point is confirmed and the variation breakpoint is located, otherwise the fragments with no difference are merged to simplify the data; the above segmentation and testing process is recursively executed until no further segmentation is possible.

[0019] Furthermore, in step 7, the RD signal in the one-dimensional space is transformed into a two-dimensional profile combining the RD signal and the MQ signal. , RD signal reflects the depth of read segments, and MQ signal mapping quality, i.e.,

[0020] in, Indicates RD value of a window; Indicates MQ value of a window; represents the number of genomic windows in the sorted BAM file generated for all regions; At the same time, the RD signal and the MQ signal are normalized, that is,

[0021] in, represents the normalized RD signal; It represents the RD value after GC correction and smoothing noise reduction; represents the mean value of RD signal in all windows; represents the standard deviation of the RD signal;

[0022] in, represents the standardized MQ signal; It represents the MQ value after GC correction and smoothing noise reduction; represents the mean of MQ signals in all windows; Indicates the standard deviation of the MQ signal.

[0023] Furthermore, in step eight, the MCD anomaly detection method is used to identify anomalies in the data, and the outliers are determined by calculating the covariance matrix of the data. In the initialization stage, an initial subset is selected and the covariance matrix and mean vector of the initial subset are calculated; in the iterative optimization stage, the subset is optimized by an iterative method, and the subset members are repeatedly adjusted until the covariance matrix determinant of the subset reaches the minimum value; in the anomaly detection stage, the Mahalanobis distance of each data point in the data set is calculated using the optimized covariance matrix and mean to determine whether each data point is The calculation process of outliers is as follows: The standardized RD and MQ signals are used as feature values ​​and input into the MCD anomaly detection method to calculate the mean vector ,

[0024]

[0025] in, represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; represents the mean value of RD signal in all windows; represents the mean of MQ signals in all windows; represents the number of genomic windows in the sorted BAM file generated for all regions; Indicates RD value of a window; Indicates MQ value of a window; Calculate the covariance matrix ,

[0026] in, represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal; represents the variance of the RD signal; represents the variance of the MQ signal; represents the covariance of RD signal and MQ signal; Calculate Mahalanobis distance ,

[0027] in, Indicates The Mahalanobis distance of the window; Indicates The feature vector of the window includes and ,Right now ; represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; Represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal.

[0028] Furthermore, the binary combination weighted evaluation method is used to weight the original anomaly score and the normalized anomaly score to obtain the comprehensive anomaly score, that is,

[0029] in, represents the comprehensive anomaly score; Represents the raw anomaly score, i.e., the Mahalanobis distance ; represents the normalized anomaly score; and represents two different weights, and .

[0030] Furthermore, in step nine, the box plot method is used to effectively identify and mark outliers. The interquartile range is calculated as the difference between the upper quartile and the lower quartile, that is,

[0031] in, represents the interquartile range; Indicates the upper four digits of a fraction; Indicates the next four digits of the fraction; Introducing parameters , and calculate the upper and lower limits, that is,

[0032] in, Indicates the upper limit value; Indicates the upper four digits of a fraction; Indicates parameters; represents the interquartile range;

[0033] in, Indicates the lower limit value; Indicates the upper four digits of a fraction; Indicates parameters; represents the interquartile range; Data points that fall above the upper limit or below the lower limit are considered outliers.

[0034] Furthermore, in step 10, the SR-PEM collaborative optimization method includes two stages. In the first stage, the CIGAR field in the BAM file is parsed using a split reads-based detection method; the CIGAR field is used as a key information carrier to indicate the alignment status of the sequencing read segment after alignment with the reference genome. The combination of its symbols and numbers is shown below:

[0035] in, The portion indicating the length of the digits preceding the symbol matches the reference position; Indicates mismatch and performs soft clipping; Indicates mismatch and hard skip; Indicates the number of bases that the sequencing reads completely match with the reference genome; Indicates the number of bases that the sequencing reads do not match the reference genome; represents the length of short sequencing reads; In the second stage, based on the completion of CIGAR field analysis, a detection method based on paired-end mapping was used. The operation steps are as follows: accurately extract paired-end reads spanning the junctions of tandem repeats from the BAM file; strictly align the extracted paired-end reads with the reference sequence to accurately identify and mark the regions where inconsistent reads are formed; collect all inconsistent reads as important clues for analyzing the structural characteristics of tandem repeat sequences; based on the collected clues, further accurately identify potential tandem repeat boundaries, thereby obtaining refined tandem repeat regions.

[0036] The present invention proposes a gene sequence variation detection method based on robust statistics and multi-strategy fusion, namely MCD-TD, which aims to solve the problem that the ratio of normal genes to abnormal genes in samples is large and the location of variation breakpoints is difficult to accurately determine. MCD-TD uses sequencing data generated by second-generation sequencing technology to detect tandem repeat variations in a single sequencing sample, and has the following technical effects: (1) The present invention adopts the MCD anomaly detection method to detect tandem repeats. By utilizing its sensitivity to characteristic signal changes, MCD-TD can more effectively detect tandem repeats with low coverage; (2) The present invention not only extracts the RD signal as a feature value, but also extracts the MQ signal as a feature value, which can effectively characterize the abnormality of each region, thereby improving the detection ability of tandem duplications; (3) The present invention integrates three strategies of the second-generation sequencing technology, namely, a detection method based on double-end mapping, a detection method based on split reads, and a detection method based on read depth, to filter out false positives and reduce boundary deviations, thereby showing superior performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a flow chart for realizing the present invention; Figure 2 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD when the simulated sequencing sample coverage is 4X and the tumor purity is 0.3; Figure 3 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 4X and tumor purity of 0.4; Figure 4 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 4X and tumor purity of 0.5; Figure 5 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 4X and tumor purity of 0.6; Figure 6 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 4X and tumor purity of 0.7; Figure 7 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 4X and tumor purity of 0.8; Figure 8 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.3; Fig. 9 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.4; Fig.10 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.5; Fig.11 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.6; Fig.12 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.7; Fig.13It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 6X and tumor purity of 0.8; Fig.14 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.3; Fig.15 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.4; Fig.16 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.5; Fig.17 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.6; Fig.18 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.7; Fig.19 It is a contour map comparing the detection results of the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD under the conditions of simulated sequencing sample coverage of 8X and tumor purity of 0.8; Fig. 20 A histogram showing the number of correctly predicted tandem repeats that can be detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD when the simulated sequencing sample coverage is 4X; Fig.21 A histogram showing the number of correctly predicted tandem repeats that can be detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD when the simulated sequencing sample coverage is 6X; Fig. 22 A histogram showing the number of correctly predicted tandem repeats that can be detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD when the simulated sequencing sample coverage is 8X; Fig.23A chord diagram of the distribution of overlapping tandem repeats detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD in the case of the real sequencing sample NA19238; Fig.24 A chord diagram of the distribution of overlapping tandem repeats detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD in the case of the real sequencing sample NA19239; Fig.25 It is a chord diagram of the distribution of overlapping tandem repeats detected by the present invention and the existing detection methods Pindel-TD, SVIM and ScanITD in the real sequencing sample NA19240. DETAILED DESCRIPTION

[0038] like Figure 1 As shown, the present invention provides a tandem repeat detection method based on robust statistics and multi-strategy fusion, which is mainly achieved through the following steps.

[0039] Step 1: Use the BWA-MEM tool in the BWA software to compare and analyze the reference genome and sequencing samples to generate an initial BAM file. Then use the SAMtools software to sort the initial BAM file according to the reference genome position, and filter out paired reads and split reads to obtain a sorted BAM file containing specific reads.

[0040] Step 2: In order to ensure accurate and reliable results in bioinformatics analysis, it is necessary to first identify missing values ​​in the data and fill them with 0. However, in the reference genome, it is usually composed of four bases, namely A, T, G, and C, and there is a special character N, which indicates that the base at this position has not yet been determined and may be any of the above four bases. In this case, it is necessary to perform the operation of removing the data at the N position. The specific steps are as follows: define the window containing the N position as an abnormal window, and set the corresponding RD value to a specific negative number (such as -9999) to distinguish it; when calculating the RD value of each window, if the RD value is less than 0, the window is determined to be an abnormal window containing the N position and is removed, otherwise it is determined to be a normal window, and the corresponding RD value and the starting and ending positions of the window are recorded. After such processing, the error of the analysis can be reduced and the accuracy and credibility of the data analysis can be improved.

[0041] Step 3: Extract samples from the BAM file and calculate the read count at each position, denoted as RC, and extract the RD value of each position in the genome in the sorted BAM file, that is,

[0042] in, Indicates RD value of a window; Indicates In the window RC value of each position; represents the length of the genomic window in the sorted BAM file, which is set to 1000 bp in the present invention; Extract the MQ value for each position in the genome in the sorted BAM file, i.e.,

[0043] in, Indicates MQ value of a window; Indicates In the window The Mq value of each position; It represents the length of the genome window in the sorted BAM file, and in the present invention, the length is set to 1000 bp.

[0044] Step 4: GC content refers to the ratio of G and C in DNA or RNA sequences, while GC content deviation refers to the uneven distribution of GC content in different regions of the genome or sequencing data in the sorted BAM file. Since GC content can affect sequencing deviation in DNA sequencing, that is, GC-rich regions may cause amplification deviation during the sequencing process, resulting in a relatively high sequencing depth in the region; GC-poor regions may result in a lower sequencing depth. Therefore, it is necessary to perform GC deviation correction on the extracted RD signal to reduce the impact of GC deviation on the results, that is,

[0045] in, Indicates the corrected RD value of a window; represents the mean value of RD signal in all windows; Indicates The RD mean of windows with similar GC content; Indicates The RD value of the window.

[0046] Step 5: Due to sample contamination and sequencing technology limitations during sequencing, the original sequencing data has strong noise. However, these noises may be mistakenly identified as tandem repeats, resulting in false positive results. At the same time, the tandem repeat region may be masked by noise, resulting in the inability to correctly detect the true variation. Therefore, the total variation model is used to reduce the noise of the RD signal and the MQ signal. The total variation model is a noise reduction method based on the total change of the signal, which specifically includes the following steps: first, the gradient between adjacent points in the RD signal and the MQ signal is calculated to quantify the change of each signal; second, the gradient value is adjusted to the minimum to reduce the noise in the signal; finally, in the optimization process, the RD signal and the MQ signal are adaptively weighted, and the weight of each point is finely adjusted to effectively retain the structural characteristics of the signal. Finally, compared with simple denoising methods such as median denoising and linear smoothing, the total variation model effectively avoids destroying the information expressed by the data through the above steps while reducing noise.

[0047] Step six, the present invention adopts a two-step progressive segmentation strategy to process the genome data in the sorted BAM file. The specific implementation technology is as follows: the first step is a sliding window segmentation method, which uses a sliding window of fixed width to process the genome in the sorted BAM file from a local perspective, and splits the genome in the sorted BAM file into continuous non-overlapping windows of the same length; the second step is to use a cyclic binary segmentation method, based on the amplitude of the RD signal, to divide the genome in the sorted BAM file in windows into multiple genome fragments, and calculate the RD mean of each fragment; then, the RD means of adjacent fragments are statistically tested for significance. If a significant difference is detected, the segmentation point is confirmed and the variation breakpoint is located, otherwise the indifferent sub-segments are merged to simplify the data; the above segmentation and testing process is recursively performed until further segmentation is impossible.

[0048] Step 7: transform the RD signal in one-dimensional space into a two-dimensional profile combining the RD signal and the MQ signal , RD signal reflects the depth of read segments, and MQ signal mapping quality, i.e.,

[0049] in, Indicates RD value of a window; Indicates MQ value of a window; represents the number of genomic windows in the sorted BAM file generated for all regions; Since the data ranges of the two features RD and MQ are quite different, there is a data imbalance problem. Therefore, in order to eliminate the dimension differences and data range differences between different features, the RD signal and the MQ signal are standardized so that the data can be compared and analyzed on the same scale, that is,

[0050] in, represents the normalized RD signal; It represents the RD value after GC correction and smoothing noise reduction; represents the mean value of RD signal in all windows; represents the standard deviation of the RD signal;

[0051] in, represents the standardized MQ signal; It represents the MQ value after GC correction and smoothing noise reduction; represents the mean of MQ signals in all windows; Indicates the standard deviation of the MQ signal.

[0052] Step eight, use the MCD anomaly detection method to identify anomalies in the data, and determine the outliers by calculating the covariance matrix of the data. It includes three stages: the initialization stage is to select an initial subset and calculate the covariance matrix and mean vector of the initial subset; the iterative optimization stage is to optimize the subset through an iterative method. The specific process includes repeatedly adjusting the subset members until the covariance matrix determinant of the subset reaches the minimum value; the anomaly detection stage is to use the optimized covariance matrix and mean to calculate the Mahalanobis distance of each data point in the data set to determine whether each data point is an anomaly.

[0053] Then, the standardized RD and MQ signals are used as feature values ​​and input into the MCD anomaly detection method. First, calculate the mean vector ,Right now,

[0054] in, represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; represents the mean value of RD signal in all windows; represents the mean of MQ signals in all windows; represents the number of genomic windows in the sorted BAM file generated for all regions; Indicates RD value of a window; Indicates MQ value of a window; Next, calculate the covariance matrix ,Right now,

[0055] in, represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal; represents the variance of the RD signal; represents the variance of the MQ signal; represents the covariance of RD signal and MQ signal; Finally, calculate the Mahalanobis distance ,Right now,

[0056] in, represents the Mahalanobis distance of the th window; Indicates The feature vector of the window includes and ,Right now ; represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal, The Mahalanobis distance calculated above is itself a metric value, which can be converted into a raw anomaly score to more intuitively represent the degree of anomaly of a data point.

[0057] The binary combination weighted evaluation method is used to weight the original anomaly score and the normalized anomaly score to obtain the comprehensive anomaly score, that is,

[0058] in, represents the comprehensive anomaly score; Represents the raw anomaly score, i.e., the Mahalanobis distance ; represents the normalized anomaly score; and represents two different weights, and In the present invention Set to 0.8, Set to 0.2.

[0059] Step 9: In order to accurately declare tandem duplications based on outliers, a reasonable statistical model needs to be established to provide a cutoff point for outliers. Therefore, the box plot method is used to effectively identify and mark outliers, and there are no restrictive requirements on the data (such as the data must follow a certain distribution).

[0060] The interquartile range of the box plot method is used to detect outliers. The interquartile range is the difference between the upper quartile and the lower quartile, which includes half of all the data, that is,

[0061] in, represents the interquartile range; Indicates the upper four digits of a fraction; Indicates the next four digits of the fraction; Secondly, introduce the parameters , and calculate the upper and lower limits, that is,

[0062] in, Indicates the upper limit value; Indicates the upper four digits of a fraction; Indicates parameters, in the present invention Set to 0.5; represents the interquartile range;

[0063] in, Indicates the lower limit value; Indicates the upper four digits of a fraction; Indicates parameters, in the present invention The default value of is 0.5; represents the interquartile range; Finally, objects that fall above the upper limit of the box plot or below the lower limit are determined to be outliers. In tandem repeats, higher outliers indicate a greater probability of being outliers, so the present invention focuses on objects above the upper limit of the box plot and preliminarily locates them as rough tandem repeat regions.

[0064] Step 10, based on the RD strategy, the detection of tandem repeats is relatively easy to implement, and does not require complex processing of sequencing data, but the boundaries of its detection results cannot be accurate to the base pair level. In addition, affected by low coverage, its detection results have a high false positive rate. Therefore, the present invention uses the SR-PEM collaborative optimization method to further process the detection results to obtain a refined tandem repeat region.

[0065] In the first stage, the split reads-based detection method is used to parse the CIGAR field in the BAM file; the CIGAR field is used as a key information carrier to indicate the alignment status of the sequencing read segment and the reference genome. The combination of its symbols and numbers is shown below:

[0066] in, The portion indicating the length of the digits preceding the symbol matches the reference position; Indicates mismatch and performs soft clipping; Indicates mismatch and hard skip; Indicates the number of bases that the sequencing reads completely match with the reference genome; Indicates the number of bases that the sequencing reads do not match the reference genome; Represents the length of short sequencing reads.

[0067] By scanning the split reads in the alignment results, potential breakpoint sets were found, and it was concluded that the split reads generated by tandem repeat variations came from reads spanning the junction of repeated sequences, and these reads were post-aligned at the start position of the repeated fragment of the reference genome and pre-aligned at the end position.

[0068] In the second stage, due to the low coverage and low tumor purity, the areas that can be refined by split reads are limited. However, after the rough tandem repeat region boundaries are processed by the split reads-based detection method, some areas are still not refined. Therefore, on the basis of completing the CIGAR field analysis, the detection method based on double-end mapping is used to further process the boundaries not covered by the split reads. The specific operation steps are as follows: First, the double-end reads spanning the junction of tandem repeats are accurately extracted from the BAM file; secondly, the extracted double-end reads are strictly compared with the reference sequence to accurately identify and mark the areas where inconsistent reads are formed; then, all inconsistent reads are collected as important clues for parsing the structural characteristics of tandem repeat sequences; finally, based on the collected clues, the potential tandem repeat boundaries are further accurately identified to obtain the refined tandem repeat regions.

[0069] In order to verify the validity and fairness of the theoretical results proposed by the present invention, 50 samples were generated for each configuration during the experiment to reduce the randomness of the experiment, and it was stipulated that only when the detection result covered more than half of the true tandem repeat region, it could be recorded as a true positive. For other comparison methods, such as Pindel-TD, SVIM and ScanITD, their default parameters were used for the experiment.

[0070] In the simulation data test, MCD-TD is compared with Pindel-TD, SVIM and ScanITD in the existing detection methods in terms of accuracy, sensitivity and F1 score, namely,

[0071] Among them, TP represents the number of tandem repeats correctly predicted by the model, P represents the total number of actual tandem repeats, and FP represents the number of tandem repeats incorrectly predicted by the model.

[0072] The experimental results are as follows Figure 2-Figure 19 As shown, in each case, the MCD-TD of the present invention obtained the highest F1 score, followed by SVIM, Pindel-TD and ScanITD. In terms of accuracy and sensitivity, MCD-TD performed best in all data, followed by SVIM and Pindel-TD. It is worth noting that MCD-TD is the only method suitable for detecting all coverage and tumor purity. In contrast, the detection results of ScanITD are not ideal due to noise interference from low coverage and low tumor purity. SVIM and Pindel-TD are more effective in medium and high tumor purity data. Especially in low tumor purity data, MCD-TD showed significant superiority, with sensitivity and F1 scores exceeding those of SVIM, Pindel-TD and ScanITD. For example, at coverage of 4X, 6X, and 8X, and a tumor purity of 0.3, the F1 scores of MCD-TD were 0.70, 0.69, and 0.69, respectively, while those of ScanITD were 0.09, 0.12, and 0.21, and those of Pindel-TD were 0.21, 0.40, and 0.50, respectively. Although SVIM performed better than other methods at low tumor purity, it was still insufficient compared to MCD-TD.

[0073] In addition, the present invention also compares the number of tandem repeat sequences that can be detected and judged correctly by MCD-TD, SVIM, Pindel-TD and ScanITD. The experimental results are as follows: Figure 20 to Figure 22As shown. Among them, the first column of each group, i.e. All, represents the total number of tandem repeat sequences. With the increase of coverage depth and tumor purity, the number of tandem repeat sequences that can be detected by all algorithms increases significantly. Specifically, the MCD-TD of the present invention exhibits detection performance superior to other methods under almost all coverage levels and tumor purity conditions. For example, when the coverage is 4X and the tumor purity is 0.6, the number of tandem repeat sequences correctly detected by MCD-TD is 18, while SVIM is 12, Pindel-TD is 8, and ScanITD is 1. When the coverage is 6X and the tumor purity is 0.7, the number of tandem repeat sequences correctly detected by MCD-TD is 18, while SVIM is 15, Pindel-TD is 13, and ScanITD is 6. Similarly, when the coverage was 8X and the tumor purity was 0.8, the number of tandem repeat sequences correctly detected by MCD-TD was 19, while SVIM was 16, Pindel-TD was 14, and ScanITD was 10. It can be seen that with the increase of coverage and tumor purity, the detection performance of MCD-TD continued to improve, further proving the superiority of the MCD-TD method.

[0074] Simulation experiments can be used as an effective method to evaluate algorithm performance, but they cannot fully capture the complexity and diversity of the genome in real sequencing samples. In order to more comprehensively evaluate the real detection effect of the MCD-TD method, it was applied to analyze three real sequencing samples NA19238, NA19239, and NA19240 from the 1000 Genome Project. The tandem repeats of a single tumor sample were predicted by using whole genome sequencing data, and the overlap density score was used to verify the feasibility of the present invention and the existing detection method on real data. The higher the value of the overlap density score, the better the performance of the algorithm, as shown below,

[0075] in, represents the average number of overlapping events of one algorithm with other algorithms, It represents the ratio of the average number of overlapping events of the algorithm to its predicted number of events.

[0076] The experimental results are as follows Figure 23-Figure 25As shown, the upper half of the circle is divided into four parts, representing the four detection methods of MCD-TD, ScanITD, Pindel-TD and SVIM, and the lower half of the circle is divided into 22 parts, representing chromosomes 1 to 22. In the three real sequencing samples, the total number of overlapping tandem repeats detected by MCD-TD is the highest, showing excellent detection ability. However, for the three existing methods of SVIM, Pindel-TD and ScanITD, the total number of overlapping tandem repeats detected in different samples is ranked as follows: Fig.23 Among them, SVIM detected the most total number, followed by ScanITD, and Pindel-TD the least; Fig.24 Among them, ScanITD detected the most, followed by SVIM, and Pindel-TD still had the least; Fig.25 Among them, ScanITD detected the most total number, followed by Pindel-TD, and SVIM had the least. It is worth noting that the tandem repeat structures of different chromosomes are usually different, which may be one of the reasons for the differences in the detection results of each method. The specific overlap density scores are shown in Table 1. Table 1 Overlap density score table Real sequencing samples MCD-TD ScanITD Pindel-TD SVIM NA19238 2787.62 195.97 103.07 940.41 NA19239 380.15 167.02 52.17 106.85 NA19240 379.94 103.49 92.85 13.49 The above results show that the gene sequence variation detection method based on robust statistics and multi-strategy fusion proposed in the present invention effectively solves the problem of tandem repeat variation detection under low coverage. This experiment tested 900 simulated data sets and 3 real data sets under the premise of ensuring fairness. Among them, the 900 simulated data sets are derived from three different coverages of 4X, 6X, 8X and six tumor purities of 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8, and each combination contains 50 groups of data samples. Specifically, by combining three coverages with six tumor purities, 18 different experimental conditions are formed, and 50 groups of data samples are generated under each condition, that is, 3 coverages × 6 tumor purities × 50 groups of samples = 900 data sets, thereby obtaining a total of 900 data sets. In addition, the three real data sets correspond to real sequencing samples NA19238, NA19239, and NA19240, respectively. The test results show that the MCD-TD proposed in the present invention achieves the best balance in sensitivity, accuracy and F1 score, and is also superior to the three existing detection methods, SVIM, Pindel-TD and ScanITD, in terms of overlap density score.

[0077] In summary, it can be proved that the present invention is a reliable tool for detecting tandem repeats from second-generation sequencing data, especially in the case of low coverage where factors such as noise have a greater impact, the detection effect is still effective and reliable.

Claims

1. A gene sequence variation detection method based on robust statistics and multi-strategy fusion, characterized in that: The following steps are included: Step 1: Input the reference genome and sequencing samples, and generate a sorted BAM file through alignment; Step 2: Preprocess the missing values ​​and N positions in the data, use a 0-filling strategy to replace the missing values, and remove the N position data, where the N position refers to the position in the sequence where the specific base cannot be determined due to sequencing technology limitations or data quality issues; Step 3, aligning the genome in the sorted BAM file with the reference genome, extracting RD signal and MQ signal as feature values, wherein the RD signal reflects the reading depth, the MQ signal reflects the mapping quality, the RD signal corresponds to the numerical expression of the RD value, and the MQ signal corresponds to the numerical expression of the MQ value; Step 4, performing GC bias correction on the extracted RD signal; Step 5, smoothing and noise reduction processing is performed on the RD signal and the MQ signal; Step 6: A two-step progressive segmentation strategy is adopted to identify continuous segments with highly consistent RD values; Step 7, construct the two-dimensional profiles of the RD signal and the MQ signal, and perform standardization processing; Step 8, applying the minimum covariance determinant method to tandem duplication detection, and using a binary combination weighted evaluation method to calculate the anomaly score; Step nine, use the box plot method to set a reasonable threshold and determine the outliers; Step 10: Implement the SR-PEM collaborative optimization method to refine the tandem repeat region.

2. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 1, characterized in that: The BWA-MEM tool in the BWA software was used to compare and analyze the reference genome and sequenced samples to generate an initial BAM file. The SAMtools software was used to arrange the initial BAM file in order according to the reference genome position, and paired reads and split reads were screened to obtain a sorted BAM file containing specific reads. Data preprocessing operations were performed on the sorted BAM file to identify missing values ​​in the data and fill them with 0. The operation of removing the N-position data was to define the window containing the N position as an abnormal window and set the corresponding RD value to a specific negative number to distinguish it. When calculating the RD value of each window, if the RD value is less than 0, the window is determined to be an abnormal window containing the N position and is removed. Otherwise, it is determined to be a normal window, and the corresponding RD value as well as the starting and ending positions of the window are recorded.

3. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 2, characterized in that: Extract samples from the BAM file and calculate the read count at each position, denoted as RC, and extract the RD value of each position in the genome in the sorted BAM file, that is, in, Indicates RD value of a window; Indicates In the window RC value of each position; Indicates the length of the genomic window in the sequenced BAM file; Extract the MQ value for each position in the genome in the sorted BAM file, i.e., in, Indicates MQ value of a window; Indicates In the window The Mq value of each position; Indicates the length of the genomic window in the sequenced BAM file.

4. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 3, characterized in that: The extracted RD signal is corrected for GC bias, that is, in, Indicates the corrected RD value of a window; represents the mean value of RD signal in all windows; Indicates The RD mean of the windows with highly consistent GC content; Indicates The RD value of the window.

5. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 4, characterized in that: The total variation model is used to smooth and reduce noise of the RD signal and the MQ signal. The operation process includes calculating the gradient between adjacent points in the RD signal and the MQ signal to quantify the changes of each signal; adjusting the gradient value to the minimum to reduce the noise in the signal; During the optimization process, adaptive weight allocation is implemented for RD signals and MQ signals, and the weight of each point is finely adjusted to effectively retain the structural characteristics of the signals.

6. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 5, characterized in that: The two-step progressive segmentation strategy includes: using a sliding window segmentation method, using a sliding window of fixed width to perform local processing on the genome in the sorted BAM file, and splitting the genome in the sorted BAM file into continuous non-overlapping windows of the same length; using a cyclic binary segmentation method, according to the amplitude of the RD signal, the genome in the sorted BAM file in units of windows is segmented into multiple genome fragments, and the RD mean of each fragment is calculated; performing a statistical significance test on the RD means of adjacent fragments, if a significant difference is detected, confirming the segmentation point and locating the variation breakpoint, otherwise merging the indifferent fragments to simplify the data; recursively executing the above segmentation and testing process until further segmentation is impossible.

7. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 6, characterized in that: Transform the RD signal in one-dimensional space into a two-dimensional profile combining the RD signal and the MQ signal , RD signal reflects the depth of read segments, and MQ signal mapping quality, i.e., in, Indicates RD value of a window; Indicates MQ value of a window; represents the number of genomic windows in the sorted BAM file generated for all regions; At the same time, the RD signal and the MQ signal are normalized, that is, in, represents the normalized RD signal; It represents the RD value after GC correction and smoothing noise reduction; represents the mean value of RD signal in all windows; represents the standard deviation of the RD signal; in, represents the standardized MQ signal; It represents the MQ value after GC correction and smoothing noise reduction; represents the mean of MQ signals in all windows; Indicates the standard deviation of the MQ signal.

8. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 7, characterized in that: The MCD anomaly detection method is used to identify outliers in the data. The outliers are determined by calculating the covariance matrix of the data. In the initialization stage, an initial subset is selected and the covariance matrix and mean vector of the initial subset are calculated. In the iterative optimization stage, the subset is optimized through an iterative method, and the subset members are repeatedly adjusted until the determinant of the covariance matrix of the subset reaches the minimum value; In the anomaly detection stage, the optimized covariance matrix and mean are used to calculate the Mahalanobis distance of each data point in the data set to determine whether each data point is an anomaly. The calculation process is as follows: The standardized RD and MQ signals are used as feature values ​​and input into the MCD anomaly detection method to calculate the mean vector , in, represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; represents the mean value of RD signal in all windows; represents the mean of MQ signals in all windows; represents the number of genomic windows in the sorted BAM file generated for all regions; Indicates RD value of a window; Indicates MQ value of the window; calculate the covariance matrix , in, represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal; represents the variance of the RD signal; represents the variance of the MQ signal; represents the covariance of RD signal and MQ signal; Calculate Mahalanobis distance , in, Indicates The Mahalanobis distance of the window; Indicates The feature vector of the window includes and ,Right now ; represents the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; Represents the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal; the binary combination weighted evaluation method is used to weight the original anomaly score and the normalized anomaly score to obtain the comprehensive anomaly score, that is, in, represents the comprehensive anomaly score; Represents the raw anomaly score, i.e., the Mahalanobis distance ; represents the normalized anomaly score; and represents two different weights, and .

9. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 8, characterized in that: The box plot method is used to effectively identify and mark outliers, and the interquartile range is calculated, which is the difference between the upper quartile and the lower quartile, that is, in, represents the interquartile range; Indicates the upper four digits of a fraction; Represents the next four digits of the fraction; introduces parameters , and calculate the upper and lower limits, that is, in, Indicates the upper limit value; Indicates the upper four digits of a fraction; Represents parameters; represents the interquartile range; in, Indicates the lower limit value; Indicates the upper four digits of a fraction; Indicates parameters; represents the interquartile range; Data points that fall above the upper limit or below the lower limit are considered outliers.

10. A gene sequence variation detection method based on robust statistics and multi-strategy fusion according to claim 9, characterized in that: The SR-PEM collaborative optimization method includes two stages. In the first stage, the CIGAR field in the BAM file is parsed using a split reads-based detection method. The CIGAR field is used as a key information carrier to indicate the alignment status of the sequencing read segment after alignment with the reference genome. The combination of its symbols and numbers is shown below: in, The portion indicating the length of the digits preceding the symbol matches the reference position; Indicates mismatch and performs soft clipping; Indicates mismatch and hard skip; Indicates the number of bases that the sequencing reads completely match with the reference genome; Indicates the number of bases that the sequencing reads do not match the reference genome; Indicates the length of the short sequencing read; in the second stage, based on the completion of CIGAR field analysis, a detection method based on paired-end mapping is used. The operation steps are as follows: accurately extract paired-end reads spanning the junction of tandem repeats from the BAM file; strictly align the extracted paired-end reads with the reference sequence to accurately identify and mark the regions where inconsistent reads are formed; collect all inconsistent reads as important clues for analyzing the structural characteristics of tandem repeat sequences; based on the collected clues, further accurately identify potential tandem repeat boundaries, thereby obtaining refined tandem repeat regions.

Images