Gene Sequence Variation Detection Method Based on Robust Statistics and Multi-Strategy Fusion

Abstract

The present invention relates to the technical field of gene mutation detection, and specifically belongs to a gene sequence mutation detection method based on robust statistics and multi-strategy fusion, including: aligning to generate a sorted BAM file; preprocessing missing values and N positions in the data, replacing missing values with a 0 filling strategy, and removing data at N positions; extracting RD signals and MQ signals as feature values; performing GC bias correction on the extracted RD signals; performing smoothing and noise reduction processing on the RD signals and MQ signals; adopting a two-step progressive segmentation strategy to identify continuous segments with highly consistent RD values; constructing a two-dimensional profile of the RD signals and MQ signals and performing normalization processing; tandem repeat detection and calculation of anomaly scores; setting a threshold using the box plot method to determine outliers; and refining the tandem repeat region. The present invention effectively solves the problem of tandem repeat mutation detection under low-coverage sequencing data, and significantly improves the accuracy and reliability of detection.

Description

Technical Field

[0001] The present invention relates to the technical field of gene mutation detection in bioinformatics, and specifically relates to a gene sequence mutation 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, they are divided into single nucleotide polymorphisms and structural variations. Among them, structural variation is one of the main types of genetic variation, which usually manifests as a sequence change of more than 50bp in the biological genome, covering forms such as large fragment deletions, insertions, duplications, inversions, translocations, and copy number variations in the genomic sequence. In the genomes of the normal population, structural variations account for about 13%, and at least 240 genes show homozygous deletion polymorphisms. Due to its rich types and high proportion, it plays a key role in explaining the diversity characteristics of organisms and populations.

[0003] Tandem repeat, as a specific form of structural variation, refers to a genomic region composed of adjacent repeat units, that is, a DNA sequence continuously and closely repeats 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 intergenic regions. This distribution characteristic enables them to have a profound impact on the structure and function of DNA, RNA, and proteins. Research has confirmed that tandem repeats are closely related to a variety of diseases, such as inducing amplification diseases and causing gene silencing, and are associated with characteristics of different cancers and acute myeloid leukemia. Therefore, it is particularly necessary to accurately detect tandem repeats in gene testing, and the continuous development and gradual improvement of next-generation sequencing technology provide a solid data foundation for the detection of tandem repeats.

[0004] With the development of next-generation sequencing technology, more and more detection strategies based on next-generation sequencing data have been developed, mainly divided into four categories: detection methods based on paired-end mapping, detection methods based on split reads, detection methods based on read depth, and detection methods based on de novo assembly. Detection methods based on paired-end mapping and detection methods based on split reads are only applicable to paired-end sequencing data. The detection method based on paired-end mapping has low resolution for detecting variant sites, while the detection method based on split reads can accurately detect breakpoints but has relatively high requirements for data quality and computing resources. The detection method based on read depth is the mainstream algorithm for detecting tandem repeat variations, applicable to single and paired-end sequencing data, and judges repeat or deletion regions according to mapping depth and read distribution, but cannot accurately detect the boundary sites of variant regions.

[0005] Second-generation sequencing technology, with its relatively low sequencing cost and high sequencing throughput, has become an important tool in genomics research 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 performing well under specific conditions, all have certain limitations.

[0006] BreakDancer identifies structural variations by analyzing abnormal mapping positions of paired ends in high-throughput sequencing data. However, in low-coverage sequencing data, its detection accuracy may decrease. Pindel-TD adapts to the characteristics of tandem repeat short-read alignment and can obtain single-nucleotide resolution at breakpoints. However, in low-coverage cases, it is difficult to obtain sufficient 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 detected structural variation features and has advantages in scenarios requiring high-precision structural variation detection. However, low-coverage samples may result in insufficient feature collection or poor clustering effect, affecting the detection performance. ScanITD accurately identifies internal tandem repeats using a split-read realignment strategy but may be limited when detecting special types of tandem repeats, and its detection accuracy may be affected in low-coverage sequencing data.

[0007] In summary, most existing methods for detecting tandem repeats can only show ideal effects in the environment of high-coverage sequencing data. In low-coverage cases, due to high sensitivity to noise, they are easily interfered by noise factors such as sequencing errors, mapping errors, and GC content biases, resulting in an increase in the false positive rate and a significant reduction in detection performance. Therefore, detecting tandem repeats in low-coverage cases 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 is used to detect tandem repeat sequences in the human genome, effectively solves the problem of detecting tandem repeat variations in 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] A gene sequence variation detection method based on robust statistics and multi-strategy fusion provided by the present invention is characterized by including the following steps:

[0010] Step 1, input a reference genome and a sequencing sample, and generate a sorted BAM file by alignment;

[0011] Step 2: Preprocess the missing values and N positions in the data. Replace the missing values using the 0 filling strategy and remove the data at the N positions, where the N positions refer to the positions in the sequence where the specific bases cannot be determined due to sequencing technology limitations or data quality issues;

[0012] Step 3: Align the genome in the sorted BAM file to the reference genome, and extract the RD signal and MQ signal as feature values. Among them, the RD signal reflects the read depth, the MQ signal reflects the mapping quality, the RD signal corresponds to the numerical representation of the RD value, and the MQ signal corresponds to the numerical representation of the MQ value;

[0013] Step 4: Perform GC bias correction on the extracted RD signal;

[0014] Step 5: Perform smoothing and noise reduction processing on the RD signal and MQ signal;

[0015] Step 6: Adopt a two-step progressive segmentation strategy to identify continuous segments with highly consistent RD values;

[0016] Step 7: Construct a two-dimensional profile of the RD signal and MQ signal and perform standardization processing;

[0017] Step 8: Apply the minimum covariance determinant method to tandem repeat detection and calculate the anomaly score using the binary combination weighted evaluation method;

[0018] Step 9: Use the box plot method to set a reasonable threshold to determine outliers;

[0019] Step 10: Implement the SR-PEM collaborative optimization method to refine the tandem repeat region.

[0020] Further, in Step 1, use the BWA-MEM tool in the BWA software to perform alignment analysis on the reference genome and the sequencing sample to generate an initial BAM file. Use the SAMtools software to sort the initial BAM file in order according to the reference genome position, and screen out paired reads and split reads to obtain a sorted BAM file containing specific reads.

[0021] Further, in Step 2, perform data preprocessing operations on the sorted BAM file, identify the missing values in the data, and uniformly fill them with 0; the operation of removing the data at the N positions 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, it is determined that the window is an abnormal window containing the N position and is removed, otherwise it is determined as a normal window, and the corresponding RD value, as well as the start position and end position of the window, are recorded.

[0022] Further, in step three, extract samples from the BAM file and calculate the read count at each position, denoted as RC. Extract the RD value at each position in the genome from the sorted BAM file, i.e.,

[0023]

[0024] where, represents the RD value of the th window; represents the RC value at the th position in the th window; represents the length of the genomic window in the sorted BAM file;

[0025] Extract the MQ value at each position in the genome from the sorted BAM file, i.e.,

[0026]

[0027] where, represents the MQ value of the th window; represents the Mq value at the th position in the th window; represents the length of the genomic window in the sorted BAM file.

[0028] Further, in step four, perform GC bias correction on the extracted RD signals, i.e.,

[0029]

[0030] where, represents the RD value of the th window after correction; represents the mean value of the RD signals of all windows; represents the mean RD value of the window with a high degree of GC content consistency with the th window; represents the RD value of the th window.

[0031] Further, in step five, use the total variation model to smooth and denoise the RD signals and MQ signals. The operation process includes calculating the gradients between adjacent points in the RD signals and MQ signals to quantify the changes in their respective signals; adjusting the gradient values to the minimum to reduce the noise in the signals; and during the optimization process, implementing adaptive weight allocation for the RD signals and MQ signals and finely adjusting the weights of each point to effectively retain the structural features of the signals.

[0032] Further, in step six, the two-step progressive segmentation strategy includes using a sliding window segmentation method to locally process the genome in the sorted BAM file with a sliding window of a fixed width, splitting the genome in the sorted BAM file into consecutive non-overlapping windows of the same length; applying a circular binary segmentation method to segment the genome in the sorted BAM file in units of windows according to the amplitude of the RD signal into multiple genome segments, and calculating the RD mean of each segment; performing a statistical significance test on the RD means of adjacent segments, and if a significant difference is detected, confirming the segmentation point and locating the variant break point, otherwise merging the segments without difference to simplify the data; recursively executing the above segmentation and testing processes until no further segmentation is possible. Specifically, in step six, the two-step progressive segmentation strategy includes using a sliding window segmentation method to locally process the genome in the sorted BAM file with a sliding window of a fixed width, splitting the genome in the sorted BAM file into consecutive non-overlapping windows of the same length; applying a circular binary segmentation method to segment the genome in the sorted BAM file in units of windows according to the amplitude of the RD signal into multiple genome segments, and calculating the RD mean of each segment; performing a statistical significance test on the RD means of adjacent segments, and if a significant difference is detected, confirming the segmentation point and locating the variant break point, otherwise merging the segments without difference to simplify the data; recursively executing the above segmentation and testing processes until no further segmentation is possible.

[0033] Further, in step seven, transform the RD signal in one-dimensional space into a two-dimensional profile combining the RD signal and the MQ signal, where the RD signal reflects the read depth and the MQ signal maps the quality, that is, where the RD signal reflects the read depth and the MQ signal maps the quality, that is,

[0034]

[0035] where, represents the RD value of the th window; represents the MQ value of the th window; represents the number of genomic windows in the sorted BAM file generated for all regions;

[0036] At the same time, standardize the RD signal and the MQ signal, that is,

[0037] where, represents the standardized RD signal; represents the RD value after GC correction and smoothing denoising; represents the mean of the RD signals of all windows; represents the standard deviation of the RD signal;

[0038]

[0039] where, represents the standardized MQ signal; represents the MQ value after GC correction and smoothing denoising; represents the mean of the MQ signals of all windows; represents the standard deviation of the MQ signal.

[0040] Further, in step eight, the MCD anomaly detection method is used to identify the anomaly points 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 by 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 Mahalanobis distance of each data point in the dataset is calculated using the optimized covariance matrix and mean, and it is determined whether each data point is an outlier, and the calculation process is as follows.

[0041] The standardized RD signal and MQ signal are used as eigenvalues and input into the MCD anomaly detection method to calculate the mean vector ,

[0042]

[0043] where represents the mean vector of the two-dimensional profile of the RD signal and MQ signal; represents the mean of the RD signals of all windows; represents the mean of the MQ signals of all windows; represents the number of genomic windows in the sorted BAM files generated in all regions; represents the th window's RD value; represents the th window's MQ value;

[0044] Calculate the covariance matrix ,

[0045]

[0046] where represents the covariance matrix of the two-dimensional profile of the RD signal and MQ signal; represents the variance of the RD signal; represents the variance of the MQ signal; represents the covariance of the RD signal and MQ signal;

[0047] Calculate the Mahalanobis distance ,

[0048]

[0049] where represents the Mahalanobis distance of the th window; represents the eigenvector of the th window, including and , that is ; 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.

[0050] Furthermore, using the binary combination weighted evaluation method, after combining the original anomaly score and the normalized anomaly score with weights, a comprehensive anomaly score is obtained, that is,

[0051]

[0052] where, represents the comprehensive anomaly score; represents the original anomaly score, that is, the Mahalanobis distance ; represents the normalized anomaly score; and represent two different weights, and .

[0053] Furthermore, in step nine, the box plot method is used to effectively identify and mark outliers,

[0054] Calculate the interquartile range, which is the difference between the upper quartile and the lower quartile, that is,

[0055]

[0056] where, represents the interquartile range; represents the upper quartile; represents the lower quartile;

[0057] Introduce the parameter , and calculate the upper limit value and the lower limit value, that is,

[0058]

[0059] where, represents the upper limit value; represents the upper quartile; represents the parameter; represents the interquartile range;

[0060]

[0061] where, represents the lower limit value; represents the upper quartile; represents the parameter; represents the interquartile range;

[0062] Data points falling above the upper limit value or below the lower limit value are determined as outliers.

[0063] Further, in step ten, the SR-PEM collaborative optimization method includes two stages. In the first stage, a detection method based on split reads is used to parse the CIGAR field in the BAM file. The CIGAR field, as a key information carrier, is used to represent the alignment status after the sequencing reads are aligned with the reference genome. The combination form of its symbols and numbers is as follows:

[0064]

[0065] Among them, The part representing the length of the number before this symbol matches the reference position; Indicates a mismatch and soft clipping; Indicates a mismatch and hard skipping; Indicates the number of bases where the sequencing read is completely matched with the reference genome; Indicates the number of bases where the sequencing read does not match the reference genome; Indicates the length of the short sequencing read;

[0066] In the second stage, on the basis of completing the parsing of the CIGAR field, a detection method based on paired-end mapping is used. The operation steps are as follows: accurately extract the 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 forming inconsistent reads; 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 a refined tandem repeat region.

[0067] The present invention proposes a gene sequence variation detection method based on robust statistics and multi-strategy fusion, namely MCD-TD, aiming to solve the problems of large differences in the proportions of normal and abnormal genes in samples and difficulty in accurately determining the positions of variant breakpoints. MCD-TD uses the sequencing data generated by the second-generation sequencing technology to detect tandem repeat variations for a single sequencing sample, and has the following technical effects:

[0068] (1) The present invention uses the MCD anomaly detection method for tandem repeat detection. Utilizing its sensitivity to changes in characteristic signals, MCD-TD can more effectively detect tandem repeats with low coverage;

[0069] (2) The present invention not only extracts the RD signal as a characteristic value, but also extracts the MQ signal as a characteristic value, which can effectively characterize the anomalies in each region, thereby improving the detection ability of tandem repeats;

[0070] (3) The present invention integrates three strategies of the second-generation sequencing technology, namely, the detection method based on paired-end mapping, the detection method based on split reads, and the detection method based on read depth, to filter false positives and reduce boundary deviation, showing more superior performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 is the implementation flowchart of the present invention;

[0072] Figure 2 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.3;

[0073] Figure 3 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.4;

[0074] Figure 4 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.5;

[0075] Figure 5 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.6;

[0076] Figure 6 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.7;

[0077] Figure 7 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 4X and the tumor purity is 0.8;

[0078] Figure 8 is the contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the condition that the coverage of the simulated sequencing sample is 6X and the tumor purity is 0.3;

[0079] Figure 9This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 6X and a tumor purity of 0.4;

[0080] Figure 10 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 6X and a tumor purity of 0.5;

[0081] Figure 11 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 6X and a tumor purity of 0.6;

[0082] Figure 12 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 6X and a tumor purity of 0.7;

[0083] Figure 13 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 6X and a tumor purity of 0.8;

[0084] Figure 14 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.3;

[0085] Figure 15 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.4;

[0086] Figure 16 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.5;

[0087] Figure 17 This is the contour map for comparing the detection results of the present invention with the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.6;

[0088] Figure 18This is a contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.7;

[0089] Figure 19 This is a contour map of the comparison of the detection results between the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD under the conditions of a simulated sequencing sample coverage of 8X and a tumor purity of 0.8;

[0090] Figure 20 This is a bar chart of 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 under the condition of a simulated sequencing sample coverage of 4X;

[0091] Figure 21 This is a bar chart of 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 under the condition of a simulated sequencing sample coverage of 6X;

[0092] Figure 22 This is a bar chart of 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 under the condition of a simulated sequencing sample coverage of 8X;

[0093] Figure 23 This is a chord diagram of the overlapping tandem repeat distribution detected by the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD in the real sequencing sample NA19238;

[0094] Figure 24 This is a chord diagram of the overlapping tandem repeat distribution detected by the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD in the real sequencing sample NA19239;

[0095] Figure 25 This is a chord diagram of the overlapping tandem repeat distribution detected by the present invention and the existing detection methods Pindel-TD, SVIM, and ScanITD in the real sequencing sample NA19240. Detailed implementation manners

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

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

[0098] Step 2: To ensure accurate and reliable results in bioinformatics analysis, it is first necessary to identify missing values in the data and uniformly fill them with 0. However, in the reference genome, it is usually composed of four bases, namely A, T, G, C, and there is also a special character N, indicating that the base at this position has not 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 (e.g., -9999) for distinction; when calculating the RD value of each window, if the RD value is less than 0, determine that the window is an abnormal window containing the N position and remove it, otherwise determine it as a normal window, and record the corresponding RD value, as well as the start position and end position of the window. After such processing, the analysis error can be reduced, and the accuracy and reliability of data analysis can be improved.

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

[0100]

[0101] where, represents the RD value of the th window; represents the RC value of the th window at the th position; represents the length of the genomic window in the sorted BAM file, which is set to 1000bp in the present invention;

[0102] Extract the MQ value at each position in the genome of the sorted BAM file, that is,

[0103]

[0104] where, represents the MQ value of the th window; represents the Mq value of the th window at the th position; Indicates the length of the genomic window in the sorted BAM file, which is set to 1000bp in the present invention.

[0105] Step four, the GC content refers to the proportion of G and C in a DNA or RNA sequence, and the GC content bias refers to the phenomenon that the GC content distribution in different regions of the genome or sequencing data in the sorted BAM file is uneven. Since the GC content can affect the sequencing bias in DNA sequencing, that is, regions rich in GC may cause amplified bias during the sequencing process, resulting in a relatively high sequencing depth in that region; regions poor in GC may lead to a lower sequencing depth. Therefore, it is necessary to correct the GC bias of the extracted RD signal to reduce the impact of GC bias on the results, that is,

[0106]

[0107] where, represents the RD value of the th window after correction; represents the mean value of the RD signals of all windows; represents the mean RD value of the windows with GC content similar to that of the th window; represents the RD value of the th window.

[0108] Step five, due to sample contamination and sequencing technology limitations during the sequencing process, the original sequencing data has strong noise. However, this noise may be misidentified as tandem repeats, resulting in false positive results. At the same time, tandem repeat regions may be masked by noise, leading to the failure to correctly detect true variations. Therefore, the total variation model is used to denoise the RD signal and the MQ signal. The total variation model is a denoising method based on the total variation of the signal, which specifically includes the following steps: First, calculate the gradient between adjacent points in the RD signal and the MQ signal to quantify the change of each signal; Second, adjust the gradient value to the minimum to reduce the noise in the signal; Finally, during the optimization process, perform adaptive weight allocation on the RD signal and the MQ signal, and finely adjust the weights of each point, so as 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 can effectively avoid destroying the information expressed by the data while reducing noise through the above steps.

[0109] Step 6: The present invention adopts a two-step progressive segmentation strategy to process the genomic data in the sorted BAM file. The specific implementation techniques are as follows. The first step is the sliding window segmentation method, which processes the genome in the sorted BAM file from a local perspective using a sliding window of fixed width, and splits the genome in the sorted BAM file into consecutive non-overlapping windows of the same length. The second step is to use the circular binary segmentation method. Based on the amplitude of the RD signal, the genome in the sorted BAM file in units of windows is segmented into multiple genomic fragments, and the RD mean value of each fragment is calculated. Then, a statistical significance test is performed on the RD mean values of adjacent fragments. If a significant difference is detected, the segmentation point is confirmed and the mutation breakpoint is located; otherwise, the segments without difference are merged to simplify the data. The above segmentation and testing processes are recursively executed until no further segmentation is possible.

[0110] Step 7: Transform the RD signal in one-dimensional space into a two-dimensional profile combining the RD signal and the MQ signal , where the RD signal reflects the read depth and the MQ signal maps the quality, that is,

[0111]

[0112] where, represents the RD value of the th window; represents the MQ value of the th window; represents the number of genomic windows in the sorted BAM file generated for all regions;

[0113] Since there is a large difference in the data range between the two features of RD and MQ, there is a problem of data imbalance. Therefore, in order to eliminate the dimensional difference and data range difference 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,

[0114]

[0115] where, represents the standardized RD signal; represents the RD value after GC correction and smoothing noise reduction; represents the mean value of the RD signal for all windows; represents the standard deviation of the RD signal;

[0116]

[0117] where, represents the standardized MQ signal; represents the MQ value after GC correction and smoothing noise reduction; The mean value of the MQ signal representing all windows; The standard deviation of the MQ signal.

[0118] Step 8: Use the MCD anomaly detection method to identify the anomaly points in the data. Determine the outliers by calculating the covariance matrix of the data, including 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 determinant of the covariance matrix 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 dataset and determine whether each data point is an anomaly point.

[0119] Subsequently, use the standardized RD signal and MQ signal as eigenvalues and input them into the MCD anomaly detection method.

[0120] First, calculate the mean vector , that is,

[0121]

[0122] Among them, Represents the mean vector of the two-dimensional profile of the RD signal and the MQ signal; Represents the mean value of the RD signal of all windows; Represents the mean value of the MQ signal of all windows; Represents the number of genomic windows in the sorted BAM files generated in all regions; Represents the RD value of the th window; Represents the MQ value of the

[0123] Secondly, calculate the covariance matrix , that is,

[0124]

[0125] Among them, Represents the covariance matrix of the two-dimensional profile 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 the RD signal and the MQ signal;

[0126] Finally, calculate the Mahalanobis distance , that is,

[0127]

[0128] Among them, The Mahalanobis distance of the \(i\)th window; denotes the \(i\)th eigenvector of the window, including and , that is ; denotes the mean vector of the two-dimensional profiles of the RD signal and the MQ signal; denotes the covariance matrix of the two-dimensional profiles of the RD signal and the MQ signal,

[0129] The calculated Mahalanobis distance itself is a metric value and can be converted into an original anomaly score, thus more intuitively representing the anomaly degree of data points.

[0130] Adopting the binary combination weighted evaluation method, after weighting and combining the original anomaly score and the normalized anomaly score, the comprehensive anomaly score is obtained, that is

[0131]

[0132] where denotes the comprehensive anomaly score; denotes the original anomaly score, that is the Mahalanobis distance ; denotes the normalized anomaly score; and denote two different weights, and in the present invention is set to 0.8, is set to 0.2.

[0133] Step Nine. In order to accurately declare tandem repeats according to the outliers, a reasonable statistical model needs to be established to provide a cut-off point for the outliers. Therefore, the box plot method is adopted to effectively identify and label the outliers, and there are no restrictive requirements on the data (such as the data needs to follow a certain distribution).

[0134] Using the interquartile range of the box plot method to detect the outliers, the interquartile range is the difference between the upper quartile and the lower quartile, which contains half of all the data, that is,

[0135]

[0136] where denotes the interquartile range; denotes the upper quartile; denotes the lower quartile;

[0137] Secondly, introduce the parameter , and calculate the upper limit value and the lower limit value, that is,

[0138]

[0139] Among them, represents the upper limit value; represents the upper quartile; represents a parameter, in the present invention is set to 0.5; represents the interquartile range;

[0140]

[0141] Among them, represents the lower limit value; represents the upper quartile; represents a parameter, in the present invention the default value is 0.5; represents the interquartile range;

[0142] Finally, objects falling above the upper limit value or below the lower limit value of the box plot are determined as outliers. In tandem repeats, a higher outlier indicates a greater probability of being an outlier. Therefore, the present invention focuses on objects above the upper limit value of the box plot and initially locates them as rough tandem repeat regions.

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

[0144] In the first stage, a detection method based on split reads is used to parse the CIGAR field in the BAM file; the CIGAR field, as a key information carrier, is used to represent the alignment status after the sequencing reads are aligned with the reference genome. The combination form of its symbols and numbers is as follows:

[0145]

[0146] Among them, represents the part with the length of the number before this symbol that matches the reference position; represents a mismatch and soft clipping; represents a mismatch and hard skipping; represents the number of base pairs where the sequencing read completely matches the reference genome; represents the number of base pairs where the sequencing read does not match the reference genome; represents the length of the short sequencing read.

[0147] By scanning the split reads in the alignment results, potential breakpoint sets are identified. It is summarized that the split reads generated by tandem repeat variations come from reads spanning the junction of repeat sequences, and these reads show post-alignment at the start position of the repeat fragment in the reference genome and pre-alignment at the end position.

[0148] In the second stage, affected by low coverage and low tumor purity, the regions that can be refined by split reads are limited. However, after processing the rough tandem repeat region boundaries by the split-read-based detection method, there are still some regions that are not refined. Therefore, on the basis of completing the CIGAR field parsing, the split-read-uncovered boundaries are further processed by using a detection method based on paired-end mapping. The specific operation steps are as follows: First, the paired-end reads spanning the tandem repeat junction are accurately extracted from the BAM file; Second, the extracted paired-end reads are strictly aligned with the reference sequence to accurately identify and mark the regions that form inconsistent reads; Then, all inconsistent reads are collected as important clues for parsing the structural features 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 region.

[0149] To verify the effectiveness and fairness of the theoretical results proposed by the present invention, during the experiment, 50 samples are generated for each configuration to reduce the randomness of the experiment, and it is stipulated that only when the detection result covers more than half of the true tandem repeat region can it be recorded as a true positive. For other comparison methods, such as Pindel-TD, SVIM, and ScanITD, experiments are carried out using their default parameters.

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

[0151]

[0152] Among them, TP represents the number of tandem repeats correctly predicted by the model, P represents the total number of tandem repeats that are actually tandem repeats, and FP represents the number of tandem repeats wrongly predicted by the model.

[0153] The experimental results are as Figures 2 - 19As 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 precision 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 rates and tumor purities. In contrast, due to the noise interference of low coverage and low tumor purity, the detection results of ScanITD were not ideal. SVIM and Pindel-TD were more effective in medium and high tumor purity data. Especially in low tumor purity data, MCD-TD demonstrated significant superiority, with both sensitivity and F1 score exceeding those of the three methods of SVIM, Pindel-TD, and ScanITD. For example, when the coverage was 4X, 6X, and 8X, and the tumor purity was 0.3, the F1 scores of MCD-TD were 0.70, 0.69, and 0.69 respectively, while the F1 scores of ScanITD were 0.09, 0.12, and 0.21 respectively, and the F1 scores 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.

[0154] In addition, the present invention also compared the number of tandem repeats that MCD-TD, SVIM, Pindel-TD, and ScanITD could detect and correctly judge. The experimental results are as Figures 20 to 22 shown. Among them, the first column of each group, namely All, represents the total number of tandem repeats present. As the coverage depth and tumor purity increased, the number of tandem repeats that all algorithms could detect increased significantly. Specifically, the MCD-TD of the present invention demonstrated superior detection performance compared to other methods under almost all coverage levels and tumor purity conditions. For example, when the coverage was 4X and the tumor purity was 0.6, the number of tandem repeats correctly detected by MCD-TD was 18, while that of SVIM was 12, that of Pindel-TD was 8, and that of ScanITD was 1. When the coverage was 6X and the tumor purity was 0.7, the number of tandem repeats correctly detected by MCD-TD was 18, while that of SVIM was 15, that of Pindel-TD was 13, and that of ScanITD was 6. Similarly, when the coverage was 8X and the tumor purity was 0.8, the number of tandem repeats correctly detected by MCD-TD was 19, while that of SVIM was 16, that of Pindel-TD was 14, and that of ScanITD was 10. It can be seen that as the coverage and tumor purity increase, the detection performance of MCD-TD continues to improve, further demonstrating the superiority of the MCD-TD method.

[0155] Simulation experiments can be an effective way to evaluate the performance of algorithms, but they cannot fully capture the complexity and diversity of genomes in real sequencing samples. To more comprehensively evaluate the true detection effect of the MCD-TD method, it was applied to analyze three real sequencing samples, NA19238, NA19239, and NA19240, from the 1000 Genomes Project. By using whole-genome sequencing data to predict tandem repeats in individual tumor samples and using the overlap density score to verify the feasibility of the present invention and existing detection methods on real data. The higher the value of the overlap density score, the better the performance of the algorithm, as shown below,

[0156]

[0157] wherein, represents the average number of overlapping events of one algorithm with other algorithms, represents the ratio of the average number of overlapping events of the algorithm to the number of predicted events.

[0158] The experimental results are as shown in Figures 23 - 25 The upper half of the circle is divided into four parts, representing four detection methods, MCD-TD, ScanITD, Pindel-TD, and SVIM. The lower half of the circle is divided into 22 parts, representing chromosomes 1 to 22. Among the three real sequencing samples, the total number of overlapping tandem repeat distributions detected by MCD-TD is the highest in all cases, demonstrating excellent detection ability. However, for the three existing methods, SVIM, Pindel-TD, and ScanITD, the total number of overlapping tandem repeat distributions detected in different samples is ranked as follows: in Figure 23 , the total number detected by SVIM is the largest, followed by ScanITD, and Pindel-TD is the least; in Figure 24 , the total number detected by ScanITD is the largest, followed by SVIM, and Pindel-TD is still the least; in Figure 25 , the total number detected by ScanITD is the largest, followed by Pindel-TD, and SVIM is the least. It should be noted that the tandem repeat structures of different chromosomes usually vary, which may be one of the reasons for the differences in the detection results of each method. The specific values of the overlap density scores are shown in Table 1,

[0159] Table 1 Overlap Density Score Table

[0160] True sequencing sample 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

[0161] The above results indicate that the gene sequence variation detection method based on robust statistics and multi-strategy fusion proposed by the present invention effectively solves the problem of tandem repeat variation detection under low coverage. Under the premise of ensuring fairness, this experiment tested 900 simulated data sets and 3 real data sets. Among them, the 900 simulated data sets were derived from the combination of three different coverages 4X, 6X, 8X and six tumor purities 0.3, 0.4, 0.5, 0.6, 0.7, 0.8. Each combination contained 50 data samples. Specifically, by combining three coverages with six tumor purities, 18 different experimental conditions were formed, and 50 data samples were generated under each condition, that is, 3 coverages × 6 tumor purities × 50 samples = 900 data sets, thus obtaining a total of 900 data sets. In addition, the 3 real data sets corresponded to the real sequencing samples NA19238, NA19239 and NA19240 respectively. The test results showed that the MCD-TD proposed by the present invention achieved the best balance in terms of sensitivity, precision and F1 score, and was also superior to the three existing detection methods SVIM, Pindel-TD and ScanITD in terms of overlapping density score.

[0162] 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, and 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. The 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; Represents 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