Annotation method for analyzing genome centromere region structure

Abstract

The invention discloses an annotation method for analyzing a genome centromere region structure, and belongs to the technical field of biological genes. And prerequisites are needed. The method comprises the following steps: firstly, analyzing repetitive monomers of a centromere sequence based on a heuristic strategy, and judging the centromere sequence; when the serial repetition units are identified, an adjacent matrix is obtained based on serial repetition unit information, the analysis problem of the HORs is converted into a graph theory problem of searching a specific structure in the adjacent matrix, and a sliding window strategy is combined to refine the structure analysis result of the HORs. Compared with the prior art, the method has the advantages that the length of the tandem repeat unit can be detected more accurately, and the defect that two continuous tandem repeat units are identified as a single tandem repeat unit by TRASH is effectively overcome. The method has an obvious efficiency advantage. When a human genome centromere area is analyzed, the method is 10 times faster than a current method while the precision is guaranteed.

Description

Technical Field

[0001] The present invention belongs to the field of biological gene technology, and in particular relates to an annotation method for analyzing the structure of a centromere region of a genome. Background Art

[0002] Centromeric regions are rich in satellite DNA sequences, characterized by long tandem repeats. These sequences typically follow a pattern of mmmm..., where m represents the repeating unit, called a monomer. The length of these monomers varies between organisms. For example, the monomer length in the centromeric region of humans is approximately 171 base pairs (bp), while in chicken it is 41 bp. In some organisms, the centromeric region may contain multiple types of monomers, with similarities ranging from 50% to 90%. These monomers are repeated in an organized manner to form higher-order tandem repeat units (HORs). For example, on chromosome X in the human genome, satellite DNA sequences follow a pattern similar to m1-m2-...-m11-m12-m1-m2...-m11-m12..., where m1-m2-...-m11-m12 is repeated multiple times to form a HOR consisting of 12 monomers.

[0003] To better understand centromere structure, accurate annotation is essential; it is a fundamental step in analyzing its structure. Centromere annotation focuses on identifying these monomers and HORs hidden within the centromere region. A key challenge in centromere annotation is repeat annotation, for which a range of tools have been specifically designed, such as Tandem RepeatFinder (TRF), RepeatExplorer2, and TAREAN. While these tools perform well for repeat annotation, they do not support the analysis of HORs. On the other hand, several tools, such as Alpha-CENTAURI, CentromereArchitect, HORmon, and HiCAT, are specifically designed for analyzing HORs within centromere regions. Although these tools can automate the analysis of HORs, they require known monomer templates, limiting their applicability to centromere structure analysis in certain non-model organisms. To address these challenges, TRASH and Satellite Repeat Finder (SRF) have been proposed. TRASH can automatically analyze monomers and infer HORs from centromere sequences, but requires users to manually select monomers from the inferred monomer set as templates for HOR analysis. SRF cannot provide detailed annotation results. Summary of the Invention

[0004] The purpose of the present invention is to provide an annotation method for analyzing the structure of the centromere region of a genome to make up for the deficiencies of the prior art.

[0005] There are two key challenges in designing centromere annotation methods. One is to accurately identify the repeat monomers in the centromere sequence, and the other is to accurately analyze the HORs structure in the centromere sequence. Accurate identification of the repeat monomers is an important prerequisite for analyzing the HOR structure.

[0006] The present invention first analyzes the repeat monomers of centromere sequences based on a heuristic strategy to determine their identity. Once the tandem repeat units are identified, an adjacency matrix is ​​generated based on the tandem repeat information. The HORs parsing problem is then transformed into a graph theory problem of finding specific structures within the adjacency matrix. A sliding window strategy is then used to refine the HORs' structural analysis results.

[0007] Specifically, the present invention is achieved through the following technical solutions:

[0008] A method for annotating the structure of centromeric regions of a genome (centroAnno) comprising the following steps:

[0009] S1: Based on the unknown sequence, first determine whether the sequence has tandem repeats, and screen for sequences with tandem repeats;

[0010] S2: Identify the high-frequency k-mers in the sequence and determine the distance between all high-frequency k-mers, which is used as the length of the potential tandem repeat unit;

[0011] S3: Based on the distance obtained in S2, a set of tandem repeat units is obtained by directly intercepting subsequences from the sequence obtained by screening in S1, and the sequence is decomposed based on this set; the sequence is finally decomposed into a list of blocks, each of which is represented as a tandem repeat unit;

[0012] S4: Based on the block list obtained in S3, detect the incorrectly annotated regions in the sequence: each incorrectly annotated block has a low score (less than 85);

[0013] S5: Re-infer tandem repeat units by clustering incorrectly annotated blocks and expand the current tandem repeat set; based on the expanded tandem repeat set, re-decompose the sequence to obtain more comprehensive and accurate annotation results;

[0014] S6: Convert the annotation results obtained in S5 into an adjacency matrix, and infer HORs based on the adjacency matrix to finally complete all annotations.

[0015] Furthermore, S1 includes: counting the frequency of occurrence of each k-mer in the sequence; for a sequence of length n, there are a total of n-k+1 k-mers; k-mers with a frequency exceeding 1 are called repeated k-mers; counting the proportion of repeated k-mers in all k-mers, when the proportion is greater than a threshold (the default is 0.2, which can also be selected according to actual conditions), the sequence will be judged as a tandem repeat sequence.

[0016] Furthermore, in S2, a distance list is calculated for all repeated k-mers. Two adjacent, identical k-mers correspond to one distance. Since repeated k-mers typically appear multiple times in a sequence, a distance list can be generated. If a k-mer appears n times in a sequence (n>1), the distance list generated by that k-mer contains n-1 elements. The frequency of occurrence of each distance is counted, where the distance with the highest frequency and the smallest occurrence represents the length of the tandem repeat unit. The term "high frequency" generally refers to a continuous string of length k that appears at least 10 times in the sequence.

[0017] Furthermore, the S3 includes:

[0018] S3-1: Use the distance d calculated in S2 to intercept the tandem repeat unit. Specifically, if the distance between adjacent and identical k-mers is d, then intercept the portion between the two k-mers in the sequence and use this portion as a candidate tandem repeat unit.

[0019] S3-2: Based on all tandem repeat units, a string decomposition algorithm is used to decompose the sequence into a list of blocks. Specifically, the string decomposition algorithm takes as input a sequence and a set of potential tandem repeat units, and outputs a list of blocks. Each block in the list has a label and a score. The label indicates the tandem repeat unit number, and the score indicates the confidence that the block is a specific tandem repeat unit.

[0020] Furthermore, the S4 includes:

[0021] S4-1: When the length of consecutive blocks exceeds a certain threshold (default is 5000), these blocks form a block region; the block regions with lower scores are detected, which indicate that the incomplete tandem repeat unit set can be inferred; in the default setting, when the score of a block is lower than 85, it is recorded as a block with lower score. For each block region with lower score, repeat the S2 step to infer the tandem repeat unit again. Here, the block score is the similarity between the block and the repeat unit, which is set by the similar length between the block and the repeat unit template.

[0022] S4-2: In reality, there are also blocks with lower scores scattered in various parts of the sequence, which do not form a block region with lower score. Extract these blocks separately and complete clustering (DBSCAN method can be selected); the consistent sequence generated by each cluster in the clustering result is taken as a new tandem repeat unit.

[0023] Further, the S5 includes: based on S4, some new tandem repeat units are obtained, which are taken as the input of the string decomposition algorithm, and the sequence is decomposed again, and the decomposition result is taken as the final tandem repeat unit annotation result.

[0024] Further, the above method further includes a HOR analysis step:

[0025] First step: construct an adjacency matrix. After S5 is completed, a block list will be obtained, and each block in the list will be assigned a label, which will induce a label sequence. For example, a block list corresponds to a label sequence 1-2-3-1-2-3..., where '1' implies that the first block of the sequence corresponds to the 1st tandem repeat unit. Correspondingly, the first block also corresponds to a score, which implies the confidence of the first block as the 1st tandem repeat unit. A label sequence S with a length of n corresponds to an n*n adjacency matrix. For each element A[i][j] of the matrix, 1

[0026] Step 2: Identify the specific structure in the adjacency matrix. Each diagonal in the matrix is ​​called a wave. For an n×n adjacency matrix, there are n waves. Starting from the main diagonal, from left to right, each wave has a length of n, n-1, ..., 2, 1. These waves are wave 1, wave 2, ..., wave n. For each wave i, if there are m consecutive elements with a value of 1, and m is not less than i, these consecutive elements are said to form a potential HOR region. The first i elements of these regions essentially correspond to a HOR. For example, for a simple block sequence 1-2-1-2, a potential HOR region of length 2 is obtained in the second wave of its corresponding adjacency matrix. The first two elements of this region correspond to "1-2", and "1-2" is a HOR.

[0027] Compared with the prior art, the advantages and beneficial effects of the present invention are:

[0028] The heuristic strategy proposed by the present invention effectively realizes the detection of tandem repeat units. Compared with the current mainstream centromere structure analysis tools, the present invention does not require similar prior information such as tandem repeat unit templates to realize the analysis of HORs. In addition, compared with the current tandem repeat analysis tools, the strategy proposed by the present invention is more accurate. Compared with the most advanced tool TRASH, the present invention can more accurately detect the length of the tandem repeat unit, effectively overcoming the weakness of TRASH that two consecutive tandem repeat units are identified as a single tandem repeat unit.

[0029] Compared to existing technologies, the present invention offers significant efficiency advantages. When analyzing the centromere region of the human genome, the present invention achieved a 10-fold increase in speed compared to current methods while maintaining high accuracy. This efficiency and accuracy advantage provides a strong foundation for analyzing tandem repeats directly within large amounts of sequencing data. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 This is a workflow diagram of the method of the present invention. a. Workflow diagram of centroAnno. b. Two methods for counting k-mers. c. The process of inferring candidate monomers. d. The process of inferring high-confidence monomers from candidate monomers involves using the same string decomposition algorithm as StringDecomposer. e. Bad blocks in the block list. f. The process of inferring HORs.

[0031] Figure 2 It is a brief diagram of the annotation method of the present invention.

[0032] Figure 3is a comparison of centromere analysis methods on 54 simulated datasets. a. Annotation rate of HOR structures in simulated datasets by different centromere analysis methods. b. Most frequent length of HORs inferred by all tools on simulated datasets. c. Most frequent length of HORs inferred by centroAnno on simulated datasets. d. Monomer length distribution inferred by centroAnno and TRASH on simulated datasets. e. Annotation rate of monomers by centroAnno and TRASH on simulated datasets. f. Runtime comparison of repeated annotation by centroAnno and TRASH on simulated datasets. The average runtime of each method on all datasets is shown in the plot. g. Memory usage comparison of repeated annotation by centroAnno and TRASH on simulated datasets. h. Runtime of HOR annotation by centromere analysis methods on all simulated datasets. The average runtime of each method on all datasets is shown in the plot. i. Memory usage of HOR annotation by centromere analysis methods on all simulated datasets.

[0033] Figure 4 is a comparison of centromere analysis methods on human CHM13-T2T genome (v2.0). a. Annotation rate of repeat units by centroAnno and TRASH. b. Length of repeats inferred by centroAnno and TRASH. c. Most frequent HOR length inferred by centromere analysis methods. “0” indicates that the analysis failed. HORmon-v1.0 indicates the HOR annotation result of HORmon on human CHM13-T2T genome (v1.0). d. Actual length distribution of most frequent HORs on chromosome 5 inferred by HiCAT. For convenience, only the first five HORs were used in the plotting process. e. Actual length distribution of most frequent HORs on chromosome 15 inferred by HiCAT. For convenience, only the first five HORs were used in the plotting process. f. Centromere regions containing alpha Satellite (aSat) and HumanSatellite 1 (HSat1). g. Length distribution of repeat units inferred by centroAnno on centromere regions containing aSat and HSat1. h. Length distribution of repeat units inferred by HiCAT on centromere regions containing aSat and HSat1. i. Runtime distribution of centromere analysis methods on human CHM13-T2T genome. j. Memory usage distribution of centromere analysis methods on human CHM13-T2T genome.

[0034] Figure 5This is the centromere structure analysis of four species. a. The reported monomer lengths in the four species. b. The monomer length distribution inferred by centroAnno. "MonoM" indicates a monomer with a length of M inferred by centroAnno. c. The repeat length distribution inferred by centroAnno on chromosome 3 of grape. d The RepBlock-MonoTemp identity (the similarity between the block unit obtained by the sequence decomposition algorithm and the inferred monomer template) distribution calculated by HiCAT on chromosome 3 of grape based on the reported (Mono107) and centroAnno inferred (Mono135) monomer templates. e. The repeat length distribution inferred by TRASH on the Y chromosome of fish. f. The RepBlock-MonoTemp identity distribution obtained by HiCAT on the Y chromosome of fish. HiCAT used the reported (Mono512) and centroAnno inferred (Mono524) templates. g. The RepBlock-MonoTemp identity distribution of three species obtained by centroAnno. The line graph shows the average RepBlock-MonoTemp identity of 12 rice chromosomes. h. RepBlock-MonoTemp identity distribution of grape obtained by centroAnno. i. Distribution of CentO numbers reported in rice and inferred by centroAnno. j. Estimated repeat unit ratio in rice. k. Heat map of the centromeric region of grape chromosome 2 (top) and the RepBlock-MonoTemp identity distribution inferred by centroAnno (bottom). l. Repeat unit ratio in Arabidopsis chromosomes. Monomer-x represents template x, and monomer-x' represents the reverse complementary sequence of x. m. Monomer array of Arabidopsis chromosome 1. Different colors in the monomer array represent monomers of different lengths. DETAILED DESCRIPTION

[0035] The technical solution of the present invention is further described below with reference to embodiments.

[0036] Example 1

[0037] Analyze the annotation of tandem repeat units in the centromeric region of the genome and the annotation method of HORs, such as Figure 1 、 Figure 2 As shown, specifically including the following:

[0038] 1. Theoretical support:

[0039] 1. Consider a tandem repeat sequence, particularly a centromere sequence. Its properties allow it to be decomposed into a series of blocks, each representing a tandem repeat unit. For example, a centromere sequence S = S1S2S3..., where Si∈{A, T, C, G}, can be decomposed into ABAB..., indicating that the sequence is composed of two tandem repeat units, A and B. From an evolutionary perspective, A and B evolved from the same "ancestral tandem repeat unit," with a similarity between A and B ranging from approximately 50% to 90%. In particular, the alternating appearance of "AB" in the sequence constitutes a high-order tandem repeat unit, or HOR.

[0040] 2. A provable proposition:

[0041] Proposition 1: Given two adjacent tandem repeat units A and B in a centromere sequence (containing N repeat units), assuming that the evolution rate is E and the length of the ancestral tandem repeat unit corresponding to A and B is L. Then the probability that A and B share a k-mer is not less than

[0042] 3. String decomposition algorithm

[0043] Given a string R (corresponding to a centromeric (tandem repeat) sequence) and a set of strings, Blocks, where each Block corresponds to a monomeric template, the goal of the string decomposition problem is to represent R using the Blocks alphabet. A chain is defined as any concatenation of Blocks, and the optimal chain for R is defined as the chain with the highest score when performing a global alignment of R, which is the best representation of R among all possible chains. The goal of the string decomposition problem is to find the optimal chain for R.

[0044] Before describing the string decomposition algorithm, we first define the alignment graph. Given a block b and a string R, their standard alignment graph consists of all vertices (i, j) where 0 ≤ i ≤ |b| and 0 ≤ j ≤ |R|, where |x| represents the length of string x. Vertex (i, j) is connected to vertices (i+1, j), (i, j+1), and (i+1, j+1), representing insertion, deletion, and match / unmatch operations, respectively. Informally, given a set of t blocks and a string R, their string decomposition graph (SD graph) consists of t standard alignment graphs, 'glued' together by their row 0. Vertices in this graph are indexed using three indices (b for a block, i for a position in block b, and j for a position in string R), except for vertices in row 0, which are simply indexed (0, j) because all blocks share the same 'gluing' row 0. Furthermore, edges are added connecting each vertex (b, |b|, j) in the last row (of each block b) with the vertex (0, j) in row 0. Formally, the SD graph is constructed on all vertices (b,i,j), where b is a block, 0≤i≤|b|, 0≤j≤|R|, and it is assumed that for each j, all vertices (b,0,j) form a single vertex (0,j). The edge set of the SD graph consists of the standard alignment edges in the standard alignment graph and the block switch edges described below. Vertex (b,i,j) is connected to vertices (b,i+1,j), (b,i,j+1), and (b,i+1,j+1), representing insertion, deletion, and match / unmatch operations, respectively. These transitions have scores of d (for insertion / deletion), r (for unmatch), and +1 (for match). In addition, vertex (b,|b|,j) is connected to vertex (0,j) of each block b in the block set via a zero-weight block switch edge. Each such edge models the transition from the end of alignment of one block to the beginning of alignment of the next block. Although the SD graph has directed cycles, the longest path in this graph (i.e., the path with the largest total edge weight) is well-defined because all directed cycles have negative weights. The vertex (0,0) is called the source, and each vertex (b,|b|,|R|) is called a sink of the SD graph.

[0045] The string decomposition algorithm is a wraparound dynamic programming algorithm for solving the string decomposition problem. The i-prefix of a block is defined as the string consisting of its first i symbols. Given a block b, a (b,i)-chain is defined as the chain ending at the i-prefix of b [note that there may be multiple (b,i)-chains]. The dynamic programming algorithm for solving the string decomposition problem is based on computing a variable score score(b,i,j), defined as the score of the best global alignment between all possible (b,i)-chains and the j-prefix of the string R. Alignments are scored using parameters d for the insertion / deletion penalty, r for the mismatch penalty, and +1 for the match reward. Similar to the fit alignment, which finds a substring of string v that has the highest-scoring alignment with the entire string w, the variable score score(b,i,j) is initialized to score(b,0,j) = 0 for all b and j, and score(b,i,0) = i*d for all i. To compute score(b,i,j), the string decomposition algorithm uses dynamic programming to calculate the length of the longest path to vertex (b,i,j) in the SD graph. The solution to the string decomposition problem is given by the longest path from the source to one of the sinks. The block b that maximizes this score represents the last block in the optimal chain. The remaining blocks in this chain are inferred by backtracing from the sink (b,|b|,|R|) to the source (0,0) and following the block switching edges defined by the backtracing path. Each path from a source to a sink in the string decomposition graph encodes an alignment of a set of blocks relative to the input string R (the score of this alignment is equal to the length of the path). Conversely, each such alignment is encoded by a path in the string decomposition graph. Therefore, the solution to the string decomposition problem is the longest path from the source to the sink in the string decomposition graph.

[0046] Figure 1 In the string decomposition algorithm, blocks are clustered using the edit distance metric, a density-based clustering algorithm is used, and a multiple sequence alignment algorithm is used to obtain the consensus sequence for each cluster as a high-confidence monomer. Bad blocks in the block list form bad regions. The process of inferring HORs: First, an adjacency matrix is ​​constructed for the block list. "1" indicates a "hit," and HORs are primarily determined from the diagonal containing the most "hits."

[0047] 2. Detection Algorithm for Tandem Repeat Units

[0048] Given a sequence, the detection of tandem repeat units is performed in three steps: first, a decision is made whether the sequence has tandem repeats, this scheme can quickly filter out centromere (tandem repeat) sequences when analyzing tandem repeat units directly in sequencing data. Second, the k-mer distances (k-mer distance refers to the number of bases between two identical k-mers on a linear sequence) are calculated and distances with high confidence are determined, which imply the length of the tandem repeat unit. Finally, the potential tandem repeat unit is truncated from the sequence based on the distances, the sequence is decomposed using the string decomposition algorithm, and the final candidate tandem repeat unit is determined based on the clustering scheme.

[0049] 1. Decision of tandem repeat sequence

[0050] One phenomenon is that sub-sequences in tandem repeat sequences will appear multiple times, which means that k-mers in tandem repeat sequences have a high frequency. Based on this, the present application designs a heuristic strategy to realize the decision of tandem repeat sequences. Specifically, the frequency of each k-mer in the sequence is counted, and the ratio of repeated k-mers (frequency greater than 1) is calculated. If the ratio is greater than a threshold r (default is 0.1), the sequence is preliminarily determined as a tandem repeat sequence. In subsequent analysis, the sequence will be decomposed into a series of blocks by the string decomposition algorithm, and each block has a score. The average value of all block scores is counted, and if the average value is greater than a threshold (default is 85), the sequence is finally determined as a tandem repeat sequence and the corresponding analysis result is output.

[0051] 2. Determining the length of the tandem repeat unit

[0052] A simple example is used to describe the rationality of the strategy of the present application for determining the length of the tandem repeat unit. A relatively small and frequent k-mer distance is likely to correspond to the typical length of a monomer. For example, consider a distance list [40, 170, 170, 170, 340, 340, 340, 340, 340]. Although 340 appears most frequently, 170 is more likely to be the length of a monomer. Because, the frequent appearance of 340 can be attributed to two non-adjacent monomers sharing a same k-mer, resulting in the appearance of 340 distance, and 170 is a factor of 340, further supporting the possibility that 170 is the length of a monomer. 40 appears less frequently, indicating that it is likely to be a random occurrence and is unlikely to represent the length of a monomer. Iteration starts from the first element of the distance list, and when a reliable distance is found, iteration is terminated, and this distance is inferred as the length of a monomer. Specifically, this problem is modeled as an optimization problem:

[0053]

[0054]

[0055] Among them, Loss(ad i ) measures ad i As a penalty for the length of the repeating unit, identity measures the similarity between two sequences, and f is the similarity threshold. i ), equal to i-1, as a penalty function, where penalty(ad i ) represents less than ad i The number of all elements, and min(1-p ij ,p ij ) means ad i It's ad j The cost of the factor. Define ad i Satisfy confidence condition C N,B , where N represents the minimum number of supporting k-mer distances, and B represents the maximum number of supporting elements for that distance (default is , if ad j With ad i If the difference does not exceed 5, it means that ad j Support ad i as the length of the tandem repeats).

[0056] 3. Determination of the final tandem repeat unit

[0057] After determining the length of the tandem repeat unit, the potential tandem repeat units are directly cut out from the sequence. Based on these tandem repeat units, the sequence is decomposed into a series of blocks in combination with the string decomposition algorithm. These blocks are clustered to obtain a series of clusters. The consensus sequence of each cluster is calculated, and these consensus sequences are used as the final tandem repeat units and the sequence is decomposed again by the string decomposition algorithm to obtain preliminary annotation results. Figure 1 f).

[0058] Before inferring HORs, we first check the score of each block in the sequence decomposition result. When there are a large number of blocks with low scores, we extract these blocks and obtain the consensus sequences of these blocks through clustering. These consensus sequences are used as new tandem repeat units, and the sequence is decomposed again based on the string decomposition algorithm to obtain the final annotation result. Based on this result, we analyze the HORs.

[0059] Each element in the list of blocks is assigned a monomer label, and the HORs are inferred from the label information.

[0060] The present application identifies potential HOR regions and infers from the label information. First, determine the wave with the largest 'hit' count, denoted as wave H, H∈{1,...,n-1}, which is considered to contain the most frequent HORs. Then, search for all 'hit' regions in wave H that satisfy the following conditions: 1≤m≤t, 0≤j≤j+l≤n-H-1, H≤k≤k+l≤n-1, l≥H-1.

[0061] Next, create a list L of length n to record the annotation information. When an element in the list is annotated, its value is modified to 1. Iteratively annotate a HOR region and update L. Then, use a sliding window strategy combined with the following constraints to identify individual HORs not in the HOR region and update L.

[0062] Finally, through the above steps, all HOR regions are iteratively identified and annotated from wave 1, wave 2,..., wave N (N is 50 by default).

[0063] As shown in Figure 2 , given a sequence with tandem repeats or high-order tandem repeat patterns, the first step is to obtain the sequence of the centromere region of the genome; the second step is to analyze the centromere region by heuristic algorithm; the third step is to process the data by combining graph theory and optimization theory; the fourth step is to output the analysis results as a structural annotation report. Finally, the centromere region in the genomic fragment is highlighted.

[0064] Example 2:

[0065] Based on the centromere structure analysis algorithm provided in Example 1, an experiment was designed to comprehensively verify the algorithm.

[0066] The annotation method centroAnno provided by the present application is superior to existing methods on simulated data sets: 54 simulated data sets were compared between centroAnno and other centromere analysis tools Figure 3). In the absence of template information, centroAnno's HOR annotation rate is comparable to HiCAT and HORmon. When template information is provided, centroAnno-mono (meaning providing a template to centroAnno) performs better on most datasets. Although centroAnno has a HOR annotation rate of 0 on some datasets, its length recognition is accurate. Compared with TRASH, centroAnno is more reliable in identifying and annotating repetitive regions, and has obvious advantages in speed and memory usage. Overall, centroAnno performs well on both simulated and real datasets, and is an efficient, accurate and low-resource utilization tool for centromere analysis.

[0067] centroAnno outperforms current methods on real human genome centromere regions: We extracted the centromere regions of all chromosomes from the human CHM13-T2T genome (v2.0) and used these regions to compare the performance of centroAnno with other centromere analysis tools ( Figure 4 CENM represents the centromere region of chromosome M, for example, CEN1 represents the centromere region of chromosome 1. CentroAnno performed excellently in repeat annotation, with a repeat annotation rate exceeding 0.9 for all centromere regions, while TRASH's annotation rate for CEN15 was less than 0.8. Furthermore, centroAnno accurately inferred repeat lengths, with most regions consistent with the actual length of 171 bp. In contrast, TRASH misinterpreted repeat lengths in some regions.

[0068] There are some differences between centroAnno and other tools in terms of HOR annotation, particularly in the slightly different HOR lengths inferred for CEN5, CEN10, and CEN15. It's worth noting that the most frequent HOR length inferred by HiCAT may not always be the actual length. Sometimes, we encounter situations where pure HOR regions cannot be obtained, including unexpected repeat regions such as HSat1. In such cases, centroAnno successfully detects all repeat types, while HORmon is unable to process them, and HiCAT misses the HSat1 repeat.

[0069] CentroAnno completed centromere analysis in an average of 190 seconds, 12 to 44 times faster than other methods. CentroAnno also used only 0.3GB of memory on average, 0.03 to 0.6 times more than other methods. CentroAnno-mono also outperformed centroAnno in both time and memory. These results demonstrate that centroAnno excels in centromere analysis, offering both efficiency and accuracy that surpasses other tools.

[0070] centroAnno can effectively analyze the centromere regions of the genomes of multiple species: in addition to the human genome, the T2T genomes of four species (Arabidopsis thaliana, rice, grape, and the Y chromosome of fish) were collected and the centromere regions of these genomes were extracted for analysis ( Figure 5 ). First, the performance of centroAnno for repeat annotation was evaluated, and the results showed that its repeat annotation rate in all centromeric regions exceeded 0.9, while TRASH performed slightly worse in some regions. Next, centroAnno was used to infer the monomers of chromosome 1 and Y chromosome of these species, and the results were consistent with previous reports. However, the monomer length inferred on grape chromosome 3 was inconsistent with previous reports. centroAnno inferred its length to be 135bp, and TRASH also inferred it to be 135bp. These two monomer lengths were input into HiCAT for further analysis, and the results showed that the RepBlock-MonoTemp identity generated by the 135bp monomer template was higher than the 107bp template, suggesting that the monomer length of grape chromosome 3 is more likely to be 135bp. In addition, the monomer length of the fish Y chromosome inferred by centroAnno was 512bp, which was different from the previous 524bp, but this difference was also supported by TRASH and HiCAT.

[0071] This example also analyzes the distribution of RepBlock-MonoTemp identities on chromosomes of various species and finds that centroAnno exhibits high RepBlock-MonoTemp identities on Arabidopsis and fish, while the values ​​on specific chromosomes of rice and grape are less than 0.9. In addition, centroAnno was used to determine the high-confidence blocks in the centromere region of each rice chromosome and to calculate the tandem repeats of each centromere region. It was also found that centromere regions with lower RepBlock-MonoTemp identities typically show lower annotation rates and lower self-similarity. Through further analysis, it was found that some centromere regions may lack self-similarity, such as the region on grape chromosome 2. The distribution of monomers in Arabidopsis and fish was also explored, and it was found that the centromere region of Arabidopsis is typically composed of 2 to 3 monomers, while the Y chromosome of fish is mainly composed of one monomer. Overall, centroAnno performs well in repeat annotation and monomer inference, and provides valuable information on the self-similarity and monomer distribution of centromere regions.

[0072] On the basis of the above embodiments, the present invention continues to describe in detail the technical features involved therein and the functions and roles played by the technical features in the present invention, so as to help technicians in this field fully understand the technical solutions of the present invention and reproduce them.

[0073] Finally, although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A method for annotating the centromere region structure of a genome, characterized in that: The annotation method comprises the following steps: S1: Based on the unknown sequence, first determine whether the sequence has tandem repeats, and screen for sequences with tandem repeats; S2: Identify the high-frequency k-mers in the sequence and determine the distance between all high-frequency k-mers, which is used as the length of the potential tandem repeat unit; S3: Based on the distance obtained in S2, a set of tandem repeat units is obtained by intercepting subsequences, and the sequence is decomposed based on this set; the sequence is finally decomposed into a list of blocks, each of which is represented as a tandem repeat unit; S4: Based on the block list obtained in S3, detect the incorrectly annotated regions in the sequence; S5: Re-infer the tandem repeat units by clustering the incorrectly annotated blocks and expand the current tandem repeat set; based on the expanded tandem repeat set, re-decompose the sequence to obtain the annotation results of the tandem repeat units; S6: Convert the annotation results obtained in S5 into an adjacency matrix, and infer HORs based on the adjacency matrix to finally complete all annotations.

2. The annotation method according to claim 1, wherein: The S1 includes: counting the frequency of occurrence of each k-mer in the sequence; for a sequence of length n, there are a total of n-k+1 k-mers; k-mers with a frequency exceeding 1 are called repeated k-mers; counting the proportion of repeated k-mers in all k-mers, when the proportion is greater than a threshold, the sequence will be judged as a tandem repeat sequence.

3. The annotation method according to claim 1, wherein: In S2, a distance list of all repeated k-mers is calculated; two adjacent and identical k-mers correspond to one distance, and since repeated k-mers appear multiple times in a sequence, a distance list can be generated; if a k-mer appears n times in a sequence, then the distance list generated by the k-mer contains n-1 elements; the frequency of occurrence of each distance is counted, where the distance with the highest frequency and the smallest frequency represents the length of the tandem repeat unit; the high frequency refers to a continuous string of length k that appears at least 10 times in the sequence.

4. The annotation method according to claim 1, wherein: The S3 includes: S3-1: The distance d calculated by S2 is used to intercept the tandem repeat unit; if the distance between adjacent and identical k-mers is d, the part between the two k-mers in the sequence is intercepted and used as the candidate tandem repeat unit; S3-2: Based on all tandem repeat units, a string decomposition algorithm is used to decompose the sequence into a block list; the string decomposition algorithm input is a sequence and a set of potential tandem repeat units, and the output is a block list. Each block in the list has a label and a score, where the label implies the tandem repeat unit number and the score represents the confidence that the block is a specific tandem repeat unit.

5. The annotation method according to claim 1, wherein: The S4 includes: S4-1: When the length of consecutive blocks exceeds a threshold, these blocks form a block region; detect the block regions with lower detection scores. For each block region with a lower score, repeat step S2 to infer tandem repeat units again. S4-2: The blocks that do not form block regions with lower scores need to be extracted separately and clustering is completed; the consensus sequences generated by each cluster in the clustering result are used as new tandem repeat units.

6. The annotation method according to claim 1, wherein: S5 includes: Based on the new tandem repeat units obtained in S4, use these tandem repeat units as the input of the string decomposition algorithm and decompose the sequence again. The decomposition result at this time is used as the final tandem repeat unit annotation result.

7. The annotation method according to claim 1, wherein: In S6, inferring HORs includes: The first step: Construct an adjacency matrix; each block in the block list obtained in S5 will be assigned a label, which can generate a label sequence; a label sequence S with a length of n corresponds to an n×n adjacency matrix. For each element A[i][j] of the matrix, where 1 < i, j ≤ n, when S[i] = S[j], A[i][j] = 1, otherwise it is 0. The second step: Determine the specific structure in the adjacency matrix; each diagonal line in the matrix is called a wave. For an n×n adjacency matrix, it has n waves; starting from the main diagonal line, the length of each wave from left to right is n, n - 1,..., 2, 1, and these waves are wave 1, wave 2,..., wave n in turn; for each wave i, there are m consecutive elements with values all equal to 1, and m is not less than i. It is said that these consecutive elements form a potential HOR region, and the first i elements corresponding to this region correspond to HOR, that is, the inference of HORs is completed.

Images