A method, apparatus, device, and storage medium for detecting tandem repeating regions.

Abstract

A method, apparatus, device, and storage medium for detecting tandem repeat regions are disclosed. The method includes: extracting read pair depth and mapping quality from a genomic sequence after processing for missing values ​​and N positions; correcting for GC bias in the read pair depth, and then subjecting the read pair depth and mapping quality to smoothing, denoising, and normalization; classifying outliers to obtain coarse tandem repeat regions; and processing the coarse tandem repeat regions using a split read pairing method and a paired end mapping method to obtain refined tandem repeat regions. This method can eliminate the influence of noise caused by sequencing errors, mapping errors, and GC content bias, resulting in a balanced distribution of read pair depth and mapping quality. Furthermore, it can further refine tandem repeat regions under conditions of low coverage and low tumor purity, thereby obtaining tandem repeat regions with higher accuracy.

Description

Technical Field

[0001] This invention relates to the field of genome structure detection technology, specifically to a method, apparatus, device, and storage medium for detecting tandem repeat regions. Background Technology

[0002] Structural variations (SVs) are one of the major types of genetic variation besides single nucleotide variants and small nucleotide variants. SVs affect approximately 0.5% of the genome of a given individual and are more capable of explaining the characteristics of diversity in organisms and populations. SVs include deletions, insertions, duplications, inversions, and translocations of DNA segments. With the decreasing cost of sequencing, SV research has gradually become a hot topic in the field of genomics.

[0003] Tandem duplications (TDs) are regions of genomic DNA containing multiple adjacent copies of a sequence, comprising approximately 10% of the human genome. Due to replication errors occurring during cell division, TDs are highly variable among individuals and exhibit a high mutation rate across the genome. They can contribute to amplification diseases, gene silencing, and rapid morphological variations. Studies have shown that TDs may be an important mechanism for cancer gene activation. Furthermore, related expansions of TD genes play a significant role in the genetic etiology and phenotypic complexity of autism spectrum disorder (ASD). Therefore, TDs play a crucial role in human disease mechanisms, necessitating accurate detection of TDs to effectively identify pathogenic genes and develop targeted therapies.

[0004] Currently, there are five strategies for detecting structural variations (SVs) using next-generation sequencing (NGS) technology. These mainly include paired-end mapping (PEM), split reads (SR), read depth (RD), de novo assembly (AS), and combinations of the above strategies (CB). PEM uses the salient features of inconsistent read pairs to infer the presence of SVs, but it is not suitable for cases where the inserted sequence is larger than the average insertion, and it cannot detect fragment duplications in SVs in low-complexity regions. SR pairs cannot perfectly match reads on the reference genome for split alignment; the detection resolution can reach the base pair (bp) level, but its performance is poor for detecting large inserted sequences. RD highlights duplication or deletion regions by studying the difference between mapping depth and read distribution, but it cannot accurately detect the boundaries of SV regions. AS constructs short DNA fragments into longer continuous sequences through sequence alignment and sequence merging calculations; however, its computational cost is too high and it is not suitable for detecting long repetitive regions and heterozygous regions. The CB strategy integrates the advantages of the above four strategies, which helps to improve the accuracy and sensitivity of the test results.

[0005] In recent years, many methods for detecting TDs from NGS data have been developed. VNTRseek maps sequencing reads to a set of reference TDs and then identifies the putative TDs based on the difference between the copy number of the references and their mapped reads, but it is only suitable for detecting microsatellite TDs. ScanITD performs a stepwise seeding and rearrangement process, using string rotation to determine whether the inserted sequence is a repetitive genomic sequence or a novel sequence, accurately predicting the region of the ITD using VAF. It is often used to detect small to medium-sized repetitions (50-300 bp), but its performance is poor when the sample sequencing coverage is low. SVIM uses a graph-based clustering method and a novel SVs signature distance metric to cluster the detected SVs signatures. Although it achieves good results across all coverage ranges, its detection performance is poor at low tumor purity. TIDDIT uses inconsistency pairs and split reads (SRs) to detect the genomic location of SVs with RD signals and uses read depth information to classify variants and assess quality; however, it has a high false positive rate when sequencing data coverage is low. TARDIS integrates multiple sequence signatures and identifies and clusters potential SV regions under the maximum parsimony assumption. It then uses a probabilistic likelihood model to distinguish specific TDs by calculating the likelihood score of each SV. However, it is not suitable for detecting sequencing data with low coverage and low tumor purity, and the accuracy of the detection results is relatively low.

[0006] The above analysis shows that accurately detecting TD regions from sequencing samples with low coverage and low tumor purity is the main challenge in current research. Existing technologies face several challenges: (1) Due to the influence of noise such as sequencing errors, mapping errors, and GC content (the ratio of guanine and cytosine in the genome) bias in NGS, the actual distribution of RD signals is uncertain. Therefore, it is inefficient to detect variant regions by simply observing the variance of RD in the genome sequence; (2) Sequencing data from low coverage depth and low tumor purity are particularly sensitive to noise, and the false positive rate of their detection results will be high; (3) Gene variant regions usually account for a very small proportion of the genome to be tested, that is, the sample category distribution is unbalanced. Summary of the Invention

[0007] The purpose of this invention is to provide a method, apparatus, device, and storage medium for detecting tandem repeating regions.

[0008] The technical solution of this invention is as follows:

[0009] A method for detecting tandem repeat regions includes: extracting read pair depth and mapping quality from a genomic sequence after processing for missing values ​​and N positions; correcting for GC bias by processing the read pair depth, and then performing smoothing, denoising, and normalization processing on the mapping quality; obtaining outliers after classification processing, screening the outliers to obtain coarse tandem repeat regions; and processing the coarse tandem repeat regions using a split read pair method and a paired end mapping method to obtain refined tandem repeat regions.

[0010] As described above, the method for detecting contiguous repeating regions, specifically the operation of obtaining outliers after classification processing, involves: using the K-nearest neighbor algorithm to process the read pair depth and mapping quality after smoothing, denoising, and standardization; finding the k nearest neighbors of any sample point in the read pair depth and mapping quality; and calculating an anomaly score based on the arbitrary sample point and the k nearest neighbors. When the anomaly score is greater than the average anomaly score of all sample points, the sample point corresponding to the anomaly score is an anomaly.

[0011] The method for detecting serial repeating regions as described above, wherein the coarse serial repeating regions are processed by the split read pair method and the paired end mapping method, specifically: the read pairs containing split read pairs are matched with the coarse serial repeating regions to obtain the refined serial repeating regions; then the read pairs containing end mappings are matched with the refined serial repeating regions to obtain the refined serial repeating regions.

[0012] As described above, the method for detecting serial repetition regions includes the following steps: Filtering outliers using a box plot method; when an outlier falls above the upper limit of the box plot, it is considered an outlier; the outliers are counted to obtain the coarse serial repetition region.

[0013] As described above, in a method for detecting contiguous repeating regions, the smoothing and denoising process specifically involves: processing the mapping quality and the read pair depth after correcting GC bias using a total variational model, as shown in the formula: x represents the noisy signal, and y represents the smoothed signal. Both x and y are vectors. and T denotes matrix transpose, N r λ represents the total number of genomic regions, and λ represents the penalty parameter, which ranges from [0.15, 0.30].

[0014] In the method for detecting contiguous repeating regions as described above, the standardization process specifically involves: using a standardization method to process the read pair depth and mapping quality after smoothing and denoising, with the formula r′=(r-μ) / σ, where r represents the mapping quality or read pair depth, r′ is the standardized sample, μ represents the mean of all sample data, and σ represents the standard deviation of all sample data.

[0015] As described above, in a method for detecting tandem repeat regions, the operation of correcting GC bias for read pair depth specifically involves: calculating the GC score of each segment in the genome sequence and performing GC correction on the read pair depth of each segment, using the following formula:

[0016]

[0017] RD′ i RD represents the correction value for the read depth of the i-th segment. i This represents the original value of the read depth of the i-th segment. This represents the mean depth of all read pairs. This represents the average depth of read pairs similar to the GC score of the i-th bin.

[0018] The present invention also provides an apparatus for detecting tandem repeating regions, comprising:

[0019] The extraction module is used to process the genome sequence to obtain the read pair depth and mapping quality;

[0020] The preprocessing module is used to process the read pair depth and mapping quality to maintain data balance;

[0021] The classification processing module is used to process the read pair depth and mapping quality to obtain a coarse concatenated repeating region;

[0022] The refinement module processes the coarsely defined series repeating region to obtain a refined series repeating region.

[0023] The present invention also provides an apparatus for detecting serial repeating regions, comprising a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the above-described method for detecting serial repeating regions.

[0024] The present invention also provides a computer-readable storage medium for storing a computer program, wherein the computer program, when executed by a processor, implements the above-described method for detecting serial repeating regions.

[0025] The beneficial effects of this invention are as follows:

[0026] (1) This invention provides a method for detecting tandem repeat regions. The method corrects GC bias processing RD, and performs smoothing and noise reduction processing and standardization processing of RD and MQ. This can eliminate the influence of noise caused by sequencing error, mapping error and GC content bias, so that RD and TDs regions follow a linear relationship and RD and MQ are evenly distributed, laying a data foundation for obtaining more accurate TDs regions in the future.

[0027] (2) This invention provides a method for detecting serial repeating regions. Based on the K-nearest neighbor algorithm, RD and MQ are processed. Anomaly scores are obtained based on the sample points and their k nearest neighbors. Then, the anomaly scores are compared with the average of the anomaly scores of all sample points to classify the sample points and obtain anomaly points. This method can find out anomaly points that are easily masked. Then, anomaly points are filtered based on the box plot program to obtain TDs regions with higher accuracy.

[0028] (3) The present invention provides a method for detecting tandem repeat regions, which further refines the TDs region based on SR and PEM, and can further refine the TDs region under low coverage and low tumor purity conditions, thereby obtaining a TDs region with higher accuracy.

[0029] (4) The present invention provides a method for detecting serial repeating regions, which has good detection effect in processing and analyzing both simulated and real data. Attached Figure Description

[0030] The solutions and advantages of this application will become clear to those skilled in the art upon reading the following detailed description of preferred embodiments. The accompanying drawings are for illustrative purposes only and are not intended to limit the scope of the invention.

[0031] In the attached diagram:

[0032] Figure 1 This is a flowchart illustrating the detection method in the embodiment;

[0033] Figure 2 This is an example image of the TDs region detected based on the SR method in the embodiment;

[0034] Figure 3 This is an example image of the TDs region detected by the PEM method in the embodiment;

[0035] Figure 4 This is a summary chart of the results of KNNTD, SVIM, TARDIS, and TIDDIT on the dataset simulation in the examples;

[0036] Figure 5 This is a summary diagram of the boundary deviations of KNNTD with SVIM, TARDIS, and TIDDIT in the embodiments;

[0037] Figure 6 This is a statistical chart showing the average number of overlap events per chromosome between KNNTD and SVIM, TARDIS, and TIDDIT in real samples, as illustrated in the example.

[0038] Figure 7 This is a distribution overview of overlapping TD events detected by KNNTD, SVIM, TARDIS, and TIDDIT in NA19238, NA19239, NA19240, and HG00266 in the embodiments.

[0039] Figure 8 This is a schematic diagram of the detection device in the embodiment;

[0040] Figure 9 This is a schematic diagram of the detection device in the embodiment. Detailed Implementation

[0041] Exemplary embodiments of this disclosure will now be described in more detail with reference to the accompanying drawings.

[0042] Example

[0043] This embodiment provides a method for detecting serial repeating regions; see [link to documentation]. Figure 1 ,include:

[0044] Step 1: Extract RD and MQ from the genome sequence after processing for missing values ​​and N positions;

[0045] Step 2: After correcting the GC bias and processing RD, the MQ is smoothed, denoised, and normalized.

[0046] Step 3: After classification, outliers are obtained. These outliers are then filtered to obtain a rough TDs region.

[0047] Step 4: The coarse TDs region is processed by the SR method and the PEM method to obtain the refined TDs region.

[0048] 1. Extract RD and MQ from the genome sequence after processing for missing values ​​and N positions.

[0049] In this embodiment, the applicant of the present invention uses the reference genome (Fasta) and sequencing samples (Fastq) as source data, and extracts RD and MQ from these two data.

[0050] First, Fasta and Fastq were aligned using BWA software. Then, SAMtools software was used to extract BAM format information configuration files (BAM files) from the alignment results and sort them according to genomic location. Next, Sambamba software was used to remove duplicates from the BAM files. Finally, SAMtools software was used to extract the required inconsistent read pairs from the sorted BAM files and index them.

[0051] In the process of retrieving RD and MQ from the BAM file, missing values ​​and N positions in the data must first be processed before retrieving RD and MQ.

[0052] 1.1 Handling Missing Values ​​and N Positions

[0053] Fasta typically consists of four bases: "A", "T", "G", and "C". An additional special character, "N", indicates that the base at that position is undetermined and could be any of the four bases mentioned above. It has been observed that existing Fasta datasets contain missing positions and "N" positions. These issues can lead to mapping errors, meaning a discrepancy between the read count and the actual read count. To provide a more reasonable read count, this embodiment fills missing positions with zeros and employs a removal strategy for the "N" positions. Specific techniques can be found in the articles "A local outlier factor-based detection of copy number variations from NGS data" by YUAN et al., and "A Local Kernel Density-Based Approach to the Detection of Copy Number Variations by Using Next-Generation Sequencing Data" by LIU et al.

[0054] 1.2 Obtaining RD and MQ

[0055] Due to noise such as sequencing errors, mapping errors, and GC content bias, the variance of reads (RDs) between bins in different genomic sequences does not follow a linear relationship with the variance of read-transformation regions (TDs). Therefore, simply observing the variance of RDs in the genomic sequence to detect variant regions is inefficient. Read quotients (MQs), on the other hand, indicate the reliability of aligning each read pair to a specific position in the reference sequence. They can be used in conjunction with other features for detection.

[0056] To address the aforementioned issues, this embodiment incorporates MQ as part of the features used to detect TDs. Considering that this embodiment primarily relies on low coverage for detection and is sensitive to coverage depth, low-quality readings will not be filtered out in actual operation to avoid reducing coverage depth.

[0057] To facilitate the extraction of RD and MQ, the genome sequence can be divided into contiguous, non-overlapping bins using a fixed-width sliding window. Since the "N" position in the reference sequence does not clearly indicate whether TD events occurred in that region, bins containing "N" are filtered out.

[0058] The filtering process specifically involves setting the length of the genome bins to L. b The number of genome bins generated in all regions is N. r Subsequently, the average read count (RC) at each position within each bin is calculated as the read depth signal for that segment, which is RD. The formula for calculating RD is as follows:

[0059]

[0060] Among them, RD i Represents the RD value of the i-th bin, RC ij L represents the RC value at the j-th position in the i-th bin; b This indicates the length of the genome bins, for example, 1000 bp.

[0061] The MQ value is the average of the Mq values ​​at each position in each bin. Its calculation method is similar to RD, and the formula is as follows:

[0062]

[0063] Among them, MQ i Mq represents the MQ value of the i-th bin. ij This represents the MQ value at the j-th position in the i-th bin.

[0064] 2. After correcting the GC bias and processing RD, the MQ is smoothed, denoised, and normalized.

[0065] To eliminate noise and imbalance in the data of RD and MQ, RD and MQ need to be preprocessed to lay the foundation for obtaining TDs regions with higher accuracy in the future.

[0066] 2.1 Correction of GC Bias Processing RD

[0067] First, the RD (Reference Scale) needs to be corrected for GC bias. NGS can introduce GC content bias during PCR amplification, which affects the RD and consequently interferes with the prediction of subsequent TD (True Dimensions) regions. To eliminate GC content bias, we calculate the GC score for each bin and correct the RD for each bin using the following formula:

[0068]

[0069] Among them, RD' i and RD i These represent the corrected value and the original value of RD for the i-th bin, respectively; This represents the mean of all RD values; This represents the average RD value that is similar to the GC score of the i-th bin.

[0070] 2.2 Smoothing and Denoising Processing for MQ and RD

[0071] Because the RD and MQ of adjacent bins have a natural correlation. However, due to sequencing and mapping errors, the RD and MQ of adjacent bins may fluctuate randomly, generating noise. Therefore, it is necessary to smooth and denoise the MQ and correct the RD after GC bias processing. Thus, a regularization term needs to be added to maintain smoothness.

[0072] This embodiment utilizes the Total Variation (TV) model, or TV model for short, to smooth and denoise the MQ and RD after correcting for GC content bias. Compared with simple denoising methods such as median denoising and linear smoothing, the TV model can remove noise, preserve the boundary information of bins, and does not damage the information expressed by the data. The TV model is expressed as the solution to a minimization problem, as shown in the following formula:

[0073]

[0074] Where x represents the noisy signal; y represents the signal obtained after smoothing using the above formula; both x and y are vectors, i.e. and N r Represents the total number of bins; T represents the matrix transpose; This represents the fitting error between x and y; Represents the total variation penalty; λ represents the penalty parameter, used for... and The trade-off between these factors is considered, and the parameter λ is selected within the range of [0.15, 0.30].

[0075] 2.3 Standardized processing of RD and MQ

[0076] Because MQ and RD may suffer from data imbalance (MQ value is much larger than RD value), when searching for nearest neighbors in the feature space during computation, MQ contributes more than RD, leading to unfairness in the detection computation. Therefore, standardization is needed to smooth and denoise RD and MQ. To reduce the impact of the difference in MQ and RD ranges on subsequent calculations, a standardization method is used to process the smoothed and denoised RD and MQ. The standardization formula is as follows:

[0077] r′=(r-μ) / σ (5)

[0078] Where r represents RD or MQ; r′ is the standardized sample; μ represents the mean of all sample data; and σ represents the standard deviation of all sample data.

[0079] After classifying and processing RD and MQ, coarse TDs regions are obtained.

[0080] To obtain a coarse TDs region, the RD and MQ after classification, denoising, and normalization processes can be used to obtain a coarse TDs region.

[0081] In this embodiment, the inventors of this application use the K-nearest neighbor (KNN) algorithm to process the denoised and normalized RD and MQ, then find the k nearest neighbors of any sample point in RD and MQ, and calculate the anomaly score based on the arbitrary sample point and the k nearest neighbors. When the anomaly score is greater than the average anomaly score of all sample points, the sample point corresponding to the anomaly score is an anomaly point.

[0082] Specifically, in finding the k nearest neighbors of any sample point in RD and MQ, a high-dimensional indexed tree data structure method (k-dimensional tree, KDTree) is used to find the k nearest neighbors.

[0083] In the outlier screening process, a box plot method is used to screen outliers. When an outlier falls above the upper limit of the box plot, it is regarded as an outlier. Outliers are counted to obtain a rough TDs region.

[0084] 3.1 The KNN algorithm processes RD and MQ to identify outliers.

[0085] To address the problem that the large size of bins hinders rapid computation, the applicant of this invention utilizes the Circular Binary Segmentation (CBS) algorithm, or CBS algorithm for short, to divide bins into regions of different sizes, thereby reducing the computational complexity of the KNN algorithm.

[0086] When using the KNN algorithm to process the two feature values ​​RD and MQ after region segmentation and standardization, RD is combined with its corresponding alignment quality to reflect the magnitude of the mutation and the quality of the mapping, respectively. RD in one-dimensional space is transformed into a two-dimensional profile D combining RD and MQ, as shown in formula (6). Analyzing the mutation interval from two perspectives further improves the detection accuracy. The formula is as follows:

[0087] D = {(RD)} i MQ i )|, i∈[1,N r ]} (6)

[0088] Since the KNN algorithm is a classification method based on a certain distance metric, its core idea is to find the k nearest neighbors of any sample point in RD and MQ, and then calculate the anomaly score of the sample point based on the sample point and the k nearest neighbors. By comparing the anomaly score of the sample point with the average anomaly score of all sample points, it is determined whether it is an anomaly.

[0089] To facilitate understanding of the k nearest neighbors in the KNN algorithm by those skilled in the art, in this embodiment, the distance between the outlier and its k nearest neighbors is relatively large, and several related definitions are introduced, namely the k-distance of object p, the k-nearest neighborhood, and the outlier score.

[0090] Definition 1: For any positive integer k, the k-distance of object p is defined as equation (7).

[0091] k_dist(p)=d(p,q) (7)

[0092] Where k represents the number of nearest neighbors for any object; d(p, q) represents the distance between object p and object q, where p, q ∈ D; and object q is the k-th nearest object to p among all objects. There are three common distance metrics for calculating d(p, q). Here, we choose the Euclidean distance, as shown in the formula below.

[0093]

[0094] Where l is the dimension of the object, which is 2 here.

[0095] Definition 2: Given the k-distance of p, the k-nearest neighbor of object p contains all objects whose distance to p is not greater than k_dist(p). The set of these objects is called the k nearest neighbors of p, defined by formula (9).

[0096] N k_dist(p) ={t|d(p,t)≤k_dist(p),t∈D\{p}} (9)

[0097] Therefore, the k nearest neighbors of object p with other objects can be represented as {d(p, t) | t∈N} k_dist(p)}

[0098] To reduce computational costs and accelerate the search for the k nearest neighbors, a high-dimensional index tree data structure is used to search for the k nearest neighbors of any sample point on a large scale in the data space, with a computational cost of O(N log N).

[0099] KDTree employs a divide-and-conquer approach, alternating between the two eigenvalues ​​RD and MQ to construct a binary tree, thus dividing the entire feature space D into multiple smaller parts. Each node in a KDTree is a k-dimensional binary tree. All non-leaf nodes can be viewed as hyperplanes that divide a space into two half-spaces.

[0100] The specific workflow of KDTree is as follows:

[0101] a. Starting from the root node, recursively descend according to the segmentation dimension of that point. If the input point is smaller than the hyperplane in the segmentation dimension, proceed to the left child node; otherwise, proceed to the right child node.

[0102] b. When descending to a leaf node, this node is taken as the "current best point" by default. The distance between the current point and the target point is calculated using Euclidean distance, as shown in formula (8).

[0103] c. Backtrack forward and perform the following operations on each node visited:

[0104] c1. If the current point is closer to the input point than the "current best point", then change it to the "current best point".

[0105] c2. Check if there is a closer point in the other subtree. If so, search downwards from that node.

[0106] d. After the root node search is completed, the nearest neighbor search is performed.

[0107] Definition 3: Given the k nearest neighbors of object p, the anomaly score S of object p can be obtained according to formula (10). p As shown in formula (10).

[0108]

[0109] Where, |N k_dist(p) | represents set N k_dist(p) The cardinality, 0≤|N k_dist(p) |≤k. Furthermore, regarding d(p,t) and N k_dist(p) For more information, please see Definition 1 and Definition 2.

[0110] To reduce the impact of the value of k, this embodiment specifies that the integer k is 20% of the sample size.

[0111] Based on the above method, the outlier scores of all sample points in RD and MQ are calculated. Then, the average of all outlier scores is calculated and used as the average outlier score. When the outlier score of a sample point is greater than the average outlier score, it is considered to be outlier and is regarded as an outlier. After comparing and analyzing the outlier scores of all sample points, all outliers can be obtained.

[0112] 3.2 Box plot method for filtering outliers to obtain coarse TDs region

[0113] From the perspective of TDs (True-Valued Regions), adjacent regions in the genome are locationally related. After performing the KNN algorithm and KDtree processing, the resulting outliers can be considered as a set of potential TD-containing regions. These outliers can be merged to obtain TD-containing regions. However, because the KNN algorithm is quite sensitive to outliers, it is prone to generating many outliers that do not possess TDs during the detection process. Therefore, the aforementioned outlier screening is necessary.

[0114] In the filtering process, a box plot method was chosen for screening, specifically using the interquartile range (IQR) of the box plot to statistically detect outliers. The interquartile range is the difference between the upper quartile (QU) and the lower quartile (QL), encompassing half of all data. A parameter θ was introduced, with upper and lower limits defined as QU + θ * IQR and QL - θ * IQR, respectively. When an outlier falls above the upper limit of the box plot, it is considered an outlier and is classified as a desired Data Point (TD) object. By summing all TD objects, a coarse TD region can be obtained.

[0115] The 4SR method and the PEM method process the coarse TDs region to obtain a refined TDs region.

[0116] 4.1 Refining the TDs region based on the SR method

[0117] Since the KNN algorithm only uses RD and MQ, the detected TDs region boundaries are relatively coarse. To improve the boundary accuracy of coarse TDs regions to the base level, this embodiment uses the SR method to process the coarse TDs regions. The core idea of ​​the SR method is to perform breakpoint matching between read pairs containing SR and the coarse TDs regions to obtain a slightly refined TDs region.

[0118] Specifically, firstly, a local alignment algorithm is used to match Fastq and Fasta. If a mismatch occurs, the mismatch location is the breakpoint, which is stored in the SR (Search Request) file and saved in the BAM file. Then, when the SR is used to refine the boundary of the coarse TDs region, two reads, R1 and R2, are randomly selected from the BAM file containing the SR. R1 and R2 are matched with the coarse TDs region. When the breakpoints in R1 and R2 have matching breakpoints at both ends of the coarse TDs region (including the regions with step sizes S before and after them), the matching position of the breakpoint is the start and end position of the refined coarse TDs region.

[0119] In the BAM file, the CIGAR field indicates the alignment status of reads after comparison with Fasta. It is usually composed of numbers, M, and S. SR can be determined by whether the CIGAR field contains the "S" cutoff symbol.

[0120] In this embodiment, CIGAR mainly analyzes two types of front alignment and back alignment, as represented by formula (11):

[0121]

[0122] Wherein, "M" indicates that the portion of the number preceding the symbol matches the reference position, "S" indicates that it does not match and has been soft-trimmed; "xMyS" indicates front alignment; "ySxM" indicates back alignment; R l This represents the length of a short sequencing read, typically set to 100 bp. For example, if a paired-end short sequencing read is 100 bp long, and one read has a CIGAR field of "20S80M", this indicates that the read is back-aligned, the first 20 bp does not match the reference position, while the remaining 80 bp does.

[0123] To facilitate understanding of the above method by those skilled in the art, the applicant of this invention provides an example, see [link to example]. Figure 2 Using the length of m bins as the effective search step size S, the starting position a and ending position b of the detected coarse TDs region are taken as the search center, and the search is extended to the left and right by S. The detection range at both ends of the coarse TDs region is 2S. Within this range, R1 and R2 are matched at breakpoints. During the search, it was found that R1 and R2 both exhibited the SR phenomenon within the 2S search range near positions a and b, i.e. Figure 2The breakpoints in R1 and R2 are back aligned (ySxM) at position A and front aligned (xMyS) at position B. Therefore, the starting position for refining the coarse TDs region is A, and the ending position for refining the coarse TDs region is B. Thus, the boundary [a,b] of the coarse TDs region can be refined to [A,B], that is, [A,B] is the sub-refined TDs region obtained by using SR.

[0124] 4.2 Refining the TDs Region Based on the PEM Method

[0125] Under conditions of low coverage and low tumor purity, the SR method can only refine a limited area because after processing the coarse TDs region [a,b], some areas will remain unrefined, such as [a,B], [A,b], and [a,b].

[0126] To address this issue, the applicant of this invention further processes the boundaries not covered by SR based on the PEM method, that is, it refines the sub-refined TDs region based on PEM, and its processing operation is similar to the core idea of ​​the SR method. First, the Fastq and Fasta are matched using a local alignment algorithm. If a mismatch occurs, the mismatch location is the breakpoint, which is stored in the PEM and saved in the BAM file. Then, when processing the sub-refined TDs region using PEM, four reads are randomly selected from the BAM file containing the PEM, namely R1', R2', R3', and R4'. R1', R2', R3', and R4' are matched with the sub-refined TDs region. When the breakpoints in R1', R2', R3', and R4' have matching breakpoints at both ends of the sub-refined TDs region (including the regions before and after its step size S), the matching position of the breakpoint is the start and end position of the sub-refined TDs region after refinement.

[0127] Among them, the extracted R1', R2', R3' and R4' are within the search range of 2S of point A, because this range may contain real breakpoint information; and the absolute value of the length of R1', R2', R3' and R4' does not exceed 4 times the length of the sub-refined TDs region.

[0128] To facilitate understanding of the above method by those skilled in the art, the applicant of this invention provides an example, see [link to example]. Figure 3 The starting position A and ending position B of the detected sub-refined TDs region are used as the search center, and the region is expanded to the left and right by S. The detection range at both ends of the sub-refined TDs region is 2S. Within this range, breakpoint matching is performed on R1', R2', R3' and R4'. During the search, it was found that PEM phenomenon appeared at positions A and B of R1', R2', R3' and R4'. Therefore, the refined TDs region obtained after PEM processing is [A,B].

[0129] 5. Experiments using simulated and real data

[0130] To verify the effectiveness of the proposed method in this embodiment—the AK-nearest neighbor detection method based on tandem duplications (KNNTD)—the applicant of this invention tested it using simulated data and real data.

[0131] In simulated data experiments, KNNTD was compared with three publicly available algorithms: SVIM (structural variant identification using mapped long reads), TARDIS (targeted digital sequencing), and TIDDIT (translocations, inversions, deletions, interspersed duplications, insertions, and tandem duplications). Performance was evaluated in four areas: sensitivity, accuracy, F1-scores (harmonic averages of sensitivity and accuracy), and boundary deviation.

[0132] In real-world data, four real-world sequencing samples from the 1000 Genomes Project (http: / / www.1000genomes.org) were used to detect KNNTD. Overlap density score (ODS) was used to evaluate the performance of each method.

[0133] During the experiment, the length L of the genome bins b The sample genome length is 1000 bp, divided by L. b The number of genome bins is N, which is then rounded up. r The parameter λ of the total variation model is set to 0.25 by default. In the KNN algorithm, the integer k is taken as 20% of the sample size, i.e., 0.2N. r In the box plot, the parameter θ is set to 0.8 by default. The effective search step size in SR is S = m·L. b , where m = 5.

[0134] 5.1 Simulated Data Research Experiment

[0135] Datasets with varying tumor purity and coverage depth were generated using SinC and seqtk (1-footnotes). Tumor purity ranged from 0.2 to 0.6, increasing by 0.2 for each increment. Coverage depths were 4x, 6x, and 8x. A total of 26 TD regions, ranging in size from 10kbp to 50kbp, were embedded in each simulation configuration. A chromosome (e.g., Chr21) from hg19 was selected as the reference genome, and each configuration was repeated 50 times to reduce experimental randomness.

[0136] KNNTD and three existing methods were tested using a simulated dataset. Here, a true positive is recorded when more than half of a real TD region is covered by the detection result. To ensure method stability, sensitivity, accuracy, and F1-score are the average of 50 replicate samples. Table 1-3 shows the specific F1-scores for different configurations. Figure 4 The specific sensitivity and accuracy under different configurations are given. Figure 4 The F1 score is represented by a black curve, ranging from 0.1 to 0.9, and the coverage depth is abbreviated as Cov.

[0137] Table 1: F1-scores of the four methods on a 4X coverage sample.

[0138]

[0139] Table 2: F1-scores of the four methods on a 6X coverage sample.

[0140]

[0141] Table 3: F1-scores of the four methods on 8X coverage samples

[0142]

[0143] From Table 1-3 and Figure 4The results show that the F1-score and sensitivity of all four methods improve with increasing tumor purity. For example, at a coverage depth of 4X, KNNTD's F1-score increases from a minimum of 0.82 to 0.98 with increasing tumor purity. In terms of F1-score, in samples with a coverage depth of 4X and tumor purity of 0.2, the maximum value of 0.82 is 0.68 higher than the minimum value of 0.14; in samples with a coverage depth of 8X and tumor purity of 0.6, the maximum value of 0.97 is 0.14 higher than the minimum value of 0.83. Regarding sensitivity, KNNTD achieves the highest value across all configurations. Meanwhile, KNNTD achieves a high accuracy of 1. Considering both sensitivity and accuracy, KNNTD ranks first in F1-score. Therefore, KNNTD achieves a better balance between accuracy and sensitivity under different coverage depths and tumor purity configurations.

[0144] To better demonstrate the data distribution and improve the detection accuracy of TDs boundaries, boundary bias was plotted for 50 sets of experiments under each configuration. Boundary bias was defined as the average number (bp) of base pairs (the difference between the recalled TDs and the actual TDs boundary bias in each experiment). A smaller boundary bias indicates a more accurate TD detection method. The comparison results of the boundary biases for the four methods are shown below. Figure 5 As shown, from Figure 5 It can be seen that as the tumor purity and coverage depth increase, the boundary error of the four methods decreases accordingly. Combined with the F1-score, the KNNTD method is superior.

[0145] In summary, the detection method (KNNTD) provided in this embodiment has high detection accuracy, sensitivity, and F1-score.

[0146] 5.2 Real-world dataset research

[0147] Simulations cannot reflect the algorithm's real-world detection performance on biological data. To validate the effectiveness of KNNTD, the applicant also applied it to analyze four real sequencing samples (NA19238, NA19239, NA19240, and HG00266) from the 1000 Genomes Project. The first three samples came from the Yoruba family's YRI trio, representing the mother, father, and daughter, respectively. The last sample was from a Finnish woman of European descent. Each of the four samples was read at a length of 150 bp, and whole-genome sequencing was performed on these samples.

[0148] KNNTD was compared with three other methods (SVIM, TARDIS, and TIDDIT). Due to the lack of ground truth files containing real data, their sensitivity, accuracy, F1-score, and boundary bias could not be calculated. Therefore, we analyzed the detection performance of the four methods by calculating the ODS for each method. The calculation of ODS is given in formula (12):

[0149] ODS=Over_lap·Over_rate, (12)

[0150] Here, Over_lap represents the average number of overlapping events between one algorithm and other algorithms, and Over_rate represents the ratio of the average number of overlapping events to the number of predicted events for that algorithm. This paper assumes that overlapping events between different algorithms are true positives; therefore, Over_lap can be considered as sensitivity, and Over_rate as accuracy. Thus, the higher the ODS value, the better the algorithm's performance.

[0151] The ODS calculation results for these four methods are shown in Table 4. As can be seen from Table 4, KNNTD achieved the highest ODS in three of the real samples.

[0152] Table 4: Comparison of ODS in four whole-genome samples (22 autosomes) using four methods

[0153]

[0154] Figure 6 The average number of overlap events per chromosome for four methods was recorded across three samples. A higher number of overlap events indicates higher sensitivity. In the NA19239 sample, TARDIS achieved the highest ODS, followed by KNNTD, TIDDIT, and SVIM. (See Table 4 for details.) Figure 6 Analysis shows that while KNNTD has a lower number of overlapping events than TIDDIT, it boasts higher accuracy. Therefore, it can be concluded that the KNNTD method is relatively reliable for application to real-world data.

[0155] Figure 7 The distribution of TD overlap events detected by the four methods is shown in a chord diagram. The upper half of the circle is divided into four parts, representing KNNTD, SVIM, TARDIS, and TIDDIT, respectively. The lower half of the circle is divided into 22 parts, representing autosomes from 1 to 22. The four methods detected the highest total number of overlaps on chromosome 1, and the best results were achieved on chromosome 10. It can also be seen that KNNTD also has good detection performance, indicating that the detection method proposed in this embodiment has good detection performance.

[0156] This embodiment discloses an image classification and recognition device, see [link to relevant documentation]. Figure 8 ,include:

[0157] The extraction module is used to process the genome sequence to obtain RD and MQ;

[0158] The preprocessing module is used to process RD and MQ data and maintain data balance.

[0159] The classification processing module is used to process RD and MQ to obtain a rough TDs region;

[0160] The refinement module is used to process the coarse TDs region to obtain the refined TDs region.

[0161] The preprocessing module includes:

[0162] The GC bias correction module is used to correct the data in RD.

[0163] The noise reduction module is used to remove noise from RD and MQ and extract useful data;

[0164] The standardization module is used to standardize data in RD and MQ to maintain data balance.

[0165] The classification processing module includes:

[0166] The K-nearest neighbor algorithm module is used to calculate the k nearest neighbors and anomaly scores of RD and MQ to obtain anomaly points;

[0167] The filtering module is used to filter out outliers and obtain a rough TDs region.

[0168] The detailed modules include:

[0169] The SR method module is used to process the coarse TDs region to obtain a more refined TDs region;

[0170] The PEM method module is used to process the sub-refined TDs region to obtain the refined TDs region.

[0171] In addition, this embodiment also discloses a device for detecting tandem repeating regions, see [link to relevant documentation]. Figure 9 The system includes a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the above-described method for detecting serial repeating regions.

[0172] In addition, this embodiment also discloses a computer-readable storage medium for storing a computer program, wherein the computer program, when executed by a processor, implements the above-described method for detecting serial repeating regions.

[0173] This invention provides a method for detecting tandem repeat regions. This method corrects for GC bias in the RD (Resolution-Oriented) processing, and performs smoothing, denoising, and normalization on RD and MQ (Molecular Quotient). This eliminates the influence of noise caused by sequencing errors, mapping errors, and GC content bias, ensuring a linear relationship between RD and TDs (Tandem Repeat) regions and a balanced distribution of RD and MQ, laying a data foundation for obtaining more accurate TDs regions. Furthermore, it uses the K-nearest neighbor algorithm to process RD and MQ, obtaining anomaly scores based on sample points and their k nearest neighbors. These anomaly scores are then compared with the average anomaly scores of all sample points to classify the sample points and identify anomalies. This method can identify easily masked anomalies, and then uses a box plot procedure to filter anomalies, obtaining more accurate TDs regions. Simultaneously, it further refines the TDs regions using SR (Sequencing Error) and PEM (Protection-Oriented Mechanism), enabling further refinement of TDs regions even under low coverage and low tumor purity conditions, thus obtaining more accurate TDs regions. In addition, this method shows good detection performance in both simulated and real data processing.

[0174] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for detecting tandem repeating regions, characterized in that, include: Extract read pair depth and mapping quality from the genome sequence after processing for missing values ​​and N positions; The operations for handling missing values ​​and the N position are as follows: fill the missing positions with zeros, and apply a removal strategy to the N position; After correcting the GC bias and processing the read pair depth, the mapping quality is smoothed, denoised, and normalized. The operation to correct read pair depth for GC bias is as follows: calculate the GC score for each segment in the genome sequence, and perform GC correction on the read pair depth in each segment. The formula is: RD′ i RD represents the correction value for the read depth of the i-th segment. i This represents the original value of the read depth of the i-th segment. This represents the mean depth of all read pairs. This represents the average depth of read pairs similar to the GC score of the i-th bin; After classification, outliers are obtained. By filtering out outliers, a rough contiguous repeating region is obtained. The operation to obtain outliers is as follows: use the K-nearest neighbor algorithm to process the read pair depth and mapping quality after denoising and standardization, find the k nearest neighbors of any sample point in the read pair depth and mapping quality, and calculate the outlier score based on the arbitrary sample point and the k nearest neighbors. When the outlier score is greater than the average outlier score of all sample points, the sample point corresponding to the outlier score is an outlier. The outlier screening process is as follows: use the box plot method to screen outliers. When an outlier falls above the upper limit of the box plot, it is considered an outlier. Count the outliers to obtain a rough contiguous repeating region. The coarse serial repeating region is processed by the split read pair method and the paired end mapping method to obtain the refined serial repeating region. Specifically: the Fastq and Fasta are matched using a local alignment algorithm. If a mismatch occurs, the mismatch position is the breakpoint, indicating the existence of a split read pair, which is saved in the BAM file. Two reads, R1 and R2, are randomly selected from the BAM file containing the split read pair. R1 and R2 are matched with the coarse serial repeating region. When the breakpoints in R1 and R2 have matching breakpoints at both ends of the coarse serial repeating region, the matching position is the start and end position of the refined serial repeating region, resulting in a sub-refined serial repeating region. The read pair containing the end mapping is then matched with the sub-refined serial repeating region to obtain the refined serial repeating region.

2. The method according to claim 1, characterized in that, The smoothing and denoising process specifically involves: The total variational model is used to process the mapping quality and the read pair depth after correcting for GC bias, as shown in the formula: x represents the noisy signal, and y represents the smoothed signal. Both x and y are vectors. and T denotes matrix transpose, N r λ represents the total number of genomic regions, and λ represents the penalty parameter, which ranges from [0.15, 0.30].

3. The method according to claim 1, characterized in that, The standardization process specifically involves: The normalization method is used to process the read pair depth and mapping quality after smoothing and denoising, and the formula is r. ′ = (r-μ) / σ, where r represents the mapping quality or read pair depth. ′ It is the standardized sample, where μ represents the mean of all sample data and σ represents the standard deviation of all sample data.

4. A device for detecting serially repeating regions, characterized in that, include: The extraction module is used to extract read pair depth and mapping quality from the genome sequence after processing for missing values ​​and N positions; The operations for handling missing values ​​and the N position are as follows: fill the missing positions with zeros, and apply a removal strategy to the N position; The preprocessing module is used to correct GC bias after processing read pair depth, and then perform smoothing, noise reduction and normalization processing on the mapping quality. The GC bias correction process involves calculating the GC score for each segment in the genome sequence and performing GC correction on the read pair depth of each segment, using the following formula: RD′ i RD represents the correction value for the read depth of the i-th segment. i This represents the original value of the read depth of the i-th segment. This represents the mean depth of all read pairs. This represents the average depth of read pairs similar to the GC score of the i-th bin; The classification module is used to identify outliers after classification, filter outliers, and obtain a coarse list of repeating regions. The outlier identification process involves: using the K-nearest neighbor algorithm to process the read pair depth and mapping quality after denoising and standardization; finding the k nearest neighbors for any sample point in the read pair depth and mapping quality; and calculating an anomaly score based on the arbitrary sample point and its k nearest neighbors. If the anomaly score is greater than the average anomaly score of all sample points, then the sample point corresponding to that score is considered an outlier. The outlier filtering process involves using a box plot method to filter outliers. When an outlier falls above the upper limit of the box plot, it is considered an outlier. Outliers are then counted to obtain a coarse list of repeating regions. The refinement module processes the coarse serial repeating region using the split read pair method and the paired end mapping method to obtain a refined serial repeating region. Specifically, it uses a local alignment algorithm to match Fastq and Fasta. If a mismatch occurs, the mismatch location is the breakpoint, indicating the existence of a split read pair, which is then saved in the BAM file. Two reads, R1 and R2, are randomly selected from the BAM file containing the split read pair. R1 and R2 are matched with the coarse serial repeating region. When a breakpoint in R1 or R2 matches a breakpoint at either end of the coarse serial repeating region, the matched position is the start and end position of the refined serial repeating region, resulting in a sub-refined serial repeating region. Finally, the read pair containing the end mapping is matched with the sub-refined serial repeating region to obtain the refined serial repeating region.

5. A device for detecting tandem repeating regions, characterized in that, The system includes a processor and a memory, wherein the processor executes a computer program stored in the memory to implement a method for detecting a series repeating region as described in any one of claims 1-3.

6. A computer-readable storage medium, characterized in that, Used to store a computer program, wherein the computer program, when executed by a processor, implements a method for detecting a series repeating region according to any one of claims 1-3.

Images