A method, system, device and storage medium for comparing long genome sequences

Abstract

The invention relates to the technical field of gene sequence alignment, and in particular to a method, system, device and storage medium for aligning long genome sequences. To solve the problem of poor accuracy of alignment methods in the prior art, the invention firstly performs weighted syncmer sampling processing on a reference genome sequence and a query sequence respectively. Then, the starting position and the ending position of a seq sequence in a sampled k-mer sequence in a query sequence and a reference genome index are clarified respectively, and an accurate position alignment interval is provided. Then, a candidate position set formed based on a query ending position and a reference ending position is subjected to chaining processing to obtain an aligned gene sequence position chain, which reflects the change of gene region alignment results. Finally, different alignment results are output according to requirements. The method is applied to aligning long genome sequences, and can improve alignment accuracy, efficiency and flexibility.

Description

Technical Field

[0001] The present invention relates to the technical field of gene sequence comparison, and in particular to a method, system, device and storage medium for comparing long genome sequences. Background Art

[0002] As a key step in research based on genome sequencing technology, the sensitivity and accuracy of genome sequence alignment will directly affect the results of subsequent research. With the widespread application of third-generation sequencing technology, the genome long sequence alignment system for third-generation sequencing technology has become one of the key tools in genomic research. Improving the sensitivity and accuracy of the genome long sequence alignment system, and then improving the effect of downstream analysis such as variant detection, has become a hot topic in current research.

[0003] Currently, the Minimap2 method is commonly used for genome sequencing alignment. The Minimap2 method uses minimizer as a subsequence (k-mer) sampling method and combines it with a seed-extension strategy for sequence alignment, which greatly reduces the index storage space and improves the sequence alignment speed. However, this method is affected by adjacent sequences and may cause the minimizer sequence sampled by the index to be discarded during the query, thereby reducing the alignment sensitivity. In addition, this method discards highly repetitive minimizer sequences to improve alignment efficiency, which will affect the alignment sensitivity and accuracy in highly repetitive regions of the genome. Summary of the invention

[0004] The purpose of the present invention is to provide a method, system, device and storage medium for aligning long genome sequences.

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

[0006] A method for comparing long genome sequences, characterized by comprising the following operations:

[0007] S1. The reference genome sequence is processed by weighted syncmer sampling to obtain a number of first sampling k-mer sequences; the first sampling k-mer sequences and their corresponding positions in the reference genome sequence are converted by hash values ​​to obtain a reference genome index;

[0008] S2. The query sequence is processed by weighted syncmer sampling to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; all query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set;

[0009] S3, the candidate position set is processed by the chain formation method to obtain a chain of aligned gene sequence positions;

[0010] S4. Based on comparing the gene sequence position chain with the query start position, query end position, reference start position, and reference end position, a genome long sequence alignment result is obtained.

[0011] The specific operation of weighted syncmer sampling processing in S1 is as follows: obtain the k-mer sequence of each position of the reference genome sequence, and determine whether there is a high-frequency s-mer sequence in each k-mer sequence; if so, perform hash value conversion processing on all s-mer sequences in the current k-mer sequence and assign weights, and increase the weight of the high-frequency s-mer sequence to the corresponding weight threshold to obtain the updated weight of each s-mer sequence; after the updated weight of each s-mer sequence is converted by hash value, it is sorted from large to small according to the hash value; if the sorting position of the last sorted s-mer sequence is the same as the position offset, the current k-mer sequence is the first sampled k-mer sequence; if not, perform hash value conversion processing on all s-mer sequences in the current k-mer sequence, and if the s-mer sequence with the smallest hash value is the same as the position offset in the current k-mer sequence, the current k-mer sequence is the first sampled k-mer sequence.

[0012] A high-frequency s-mer sequence is an s-mer sequence whose occurrence frequency in the reference genome sequence is greater than a frequency threshold; the frequency threshold is based on the total frequency of the s-mer sequence in the reference genome sequence.

[0013] The specific operations of the chaining method processing in S3 are as follows: based on the candidate position set, construct a comparison position map with the query end position as the horizontal axis and the reference end position as the vertical axis; use the anchor point of the query end position that is the last in the query sequence in the comparison position map as the first chaining starting point; within the first chaining area, use the anchor point with the highest chaining score between the first chaining starting point and the neighborhood as the anchor point to be chained, and connect it with the first chaining starting point to obtain a first sequence position chain; use the anchor point to be chained as the starting point of the first sequence position chain to obtain a second chaining starting point, and perform a connection operation based on the chaining score in the second chaining area to obtain a second sequence position chain; and so on, after completing the connection of the last anchor point, obtain a compared gene sequence position chain.

[0014] The chaining score is obtained by the following formula:

[0015] ,

[0016] f(i) For the i The chain score of the anchor point is j The starting point of the chain i The chaining scores between anchor points are equal. f(j) For the j The chain score of the chain starting point is j The chaining score of the chaining starting point is equal to that of the previous chaining starting point. α(j,i) For the j The starting point of the chain i The number of matching bases between anchor points, β(j,i) For the vacancy penalty, k is the length of the k-mer sequence.

[0017] The method for obtaining the number of matching bases is as follows: obtain the query end position difference and the query end position difference between the j-th chain starting point and the i-th anchor point, and take the minimum value of the query end position difference, the query end position difference and the length of the k-mer sequence as the number of matching bases between the j-th chain starting point and the i-th anchor point.

[0018] If there are multiple aligned gene sequence position chains, the aligned gene sequence position chains whose total position chain scores are greater than the total score threshold are sequentially subjected to the operations in S4 in descending order of the total position chain scores; the total position chain score is the chain score of the last chained anchor point in the aligned gene sequence position chain.

[0019] A genome long sequence alignment system, used to implement the above-mentioned genome long sequence alignment method, comprising:

[0020] A reference genome index generation module is used to obtain a number of first sampling k-mer sequences by weighted syncmer sampling processing of the reference genome sequence; a number of first sampling k-mer sequences and corresponding positions in the reference genome sequence are converted into hash values ​​to obtain a reference genome index;

[0021] The candidate position set generation module is used to query the sequence through weighted syncmer sampling processing to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; all query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set;

[0022] A module for generating a position chain of an aligned gene sequence, which is used to process a candidate position set through a chain formation method to obtain a position chain of an aligned gene sequence;

[0023] The genome long sequence comparison result generation module is used to obtain the genome long sequence comparison result based on the comparison gene sequence position chain with the query start position, query end position, reference start position and reference end position.

[0024] A long genome sequence comparison device comprises a processor and a memory, wherein the processor implements the above-mentioned long genome sequence comparison method when executing a computer program stored in the memory.

[0025] A computer-readable storage medium is used to store a computer program, wherein the computer program implements the above-mentioned method for comparing long genome sequences when executed by a processor.

[0026] The beneficial effects of the present invention are:

[0027] The present invention provides a method for aligning a long genome sequence. First, a weighted syncmer sampling process is performed on a reference genome sequence and a query sequence respectively, and a more valuable sampled k-mer sequence is screened out for enhancing sensitivity. Then, a query start position and a query end position of a seq sequence in the sampled k-mer sequence in the query sequence are obtained to define a region that needs to be focused on in the query sequence, and a reference start position and a reference end position of the seq sequence in a reference genome index are obtained, so as to provide an accurate alignment interval for aligning the query sequence with the reference genome sequence. Next, a chaining process is performed on a candidate position set formed based on the query end position and the corresponding reference end position to obtain an aligned gene sequence position chain, which can reflect the change of the alignment result in a continuous gene region and improve the accuracy and intuitiveness of the alignment result. Finally, according to different requirements for output results, corresponding alignment result output is performed on the information related to the aligned gene sequence position chain to improve the flexibility of the alignment method. The method is applied to the alignment of long genome sequences and can improve the alignment accuracy, efficiency and flexibility. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] By reading the detailed description of the preferred embodiment below, the scheme and advantages of the present application will become clear to those skilled in the art. The accompanying drawings are only for the purpose of illustrating the preferred embodiment and are not to be considered as limiting the present invention.

[0029] In the attached picture:

[0030] Figure 1 Schematic diagram of the flow of the comparison method of this embodiment in the embodiment;

[0031] Figure 2 2 is a graph showing the sampling density of human chromosome X using the method of this embodiment and the prior art. DETAILED DESCRIPTION

[0032] Exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings.

[0033] This example provides a method for comparing long genome sequences. Figure 1 , including the following operations:

[0034] S1. The reference genome sequence is processed by weighted syncmer sampling to obtain a number of first sampling k-mer sequences; the first sampling k-mer sequences and their corresponding positions in the reference genome sequence are converted by hash values ​​to obtain a reference genome index;

[0035] S2. The query sequence is processed by weighted syncmer sampling to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position;

[0036] S3, all query end positions and corresponding reference end positions, as well as k-mer sequence lengths, form a candidate position set; the candidate position set is processed by the chain formation method to obtain a aligned gene sequence position chain;

[0037] S4. Based on comparing the gene sequence position chain with the query start position, query end position, reference start position, and reference end position, a genome long sequence alignment result is obtained.

[0038] S1. The reference genome sequence is processed by weighted syncmer sampling to obtain a number of first sampling k-mer sequences; the first sampling k-mer sequences and their corresponding positions in the reference genome sequence are converted by hash values ​​to obtain a reference genome index.

[0039] By performing weighted syncmer sampling on the reference genome sequence, several more valuable first sampling k-mer sequences are screened out. These first sampling k-mer sequences can better reflect the essential characteristics of the genome and help improve accuracy in alignment. At the same time, these first sampling k-mer sequences and their positions are converted into hash values, which can significantly reduce storage space, speed up alignment, and improve alignment efficiency.

[0040] First, the reference genome sequence Ref_seq in Fastq format is read from the hard disk. The reference genome sequence is processed by weighted syncmer sampling to obtain several first sampling k-mer sequences. The k-mer sequence refers to a subsequence of length k continuously extracted from the reference genome sequence, which consists of four bases (A, T, C, G).

[0041] The operation of weighted syncmer sampling processing is as follows: obtain the k-mer sequence of each position of the reference genome sequence, and determine whether there is a high-frequency s-mer sequence in each k-mer sequence (the s-mer sequence is a subsequence of length s from the reference genome sequence); if so, perform hash value conversion processing on all s-mer sequences in the current k-mer sequence and assign weights, and increase the weight of the high-frequency s-mer sequence to the corresponding weight threshold to obtain the updated weight of each s-mer sequence; after the updated weight of each s-mer sequence is converted by hash value, it is sorted from large to small according to the hash value. If the sorting position of the last sorted s-mer sequence is the same as the position offset, the current k-mer sequence is the first sampled k-mer sequence; if not, perform hash value conversion processing on all s-mer sequences in the current k-mer sequence. If the s-mer sequence with the smallest hash value is in the same position as the position offset in the current k-mer sequence, the current k-mer sequence is the first sampled k-mer sequence.

[0042] The above-mentioned high-frequency s-mer sequence is an s-mer sequence whose frequency of occurrence in the reference genome sequence is greater than the frequency threshold; the frequency threshold is obtained based on the total frequency of the s-mer sequence in the reference genome sequence, and is preferably 0.02% of the total frequency. There are many types of s-mer sequences. Taking the s-mer sequences with high occurrence frequency as high-frequency s-mer sequences for weight adjustment can filter out the highly repeated syncmer sequences in the first sampled k-mer sequence, reduce the repeatability of the syncmer sequence, and thus improve the sensitivity and accuracy of long sequence alignment.

[0043] Then, several first sampled k-mer sequences and their corresponding positions in the reference genome sequence are converted into hash values, that is, several first sampled k-mer sequences and their corresponding positions are stored in a hash table to obtain a reference genome index. The key value in the reference genome index is the base sequence of the first sampled k-mer sequence (which can be called a weighted syncmer), and the value is the set of positions corresponding to the first sampled k-mer sequence (weighted syncmer) in the reference genome sequence. At the same time, in order to provide base-level sequence alignment results, the reference genome index retains the binary coding sequence of the reference genome sequence, and the four bases A, C, G, and T are encoded as 00, 01, 10, and 11, respectively. For example, the sequence AACGTTA is encoded as 00000110111100.

[0044] The specific processing process of S1 can refer to the logical algorithm flow in Table 1.

[0045] Table 1 Logical algorithm flow for obtaining reference genome index

[0046] .

[0047] S2. The query sequence is processed by weighted syncmer sampling to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; all query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set.

[0048] After the query sequence is processed by weighted syncmer sampling, the query start position and query end position are obtained based on the seq sequence in the second sampling k-mer sequence, the area that needs to be focused on in the query sequence is accurately defined, and the reference start position and reference end position of the seq sequence in the second sampling k-mer sequence in the reference genome index are obtained, which provides an accurate comparison interval for the comparison between the query sequence and the reference genome sequence, thereby improving the comparison accuracy and efficiency.

[0049] First, a query sequence Read_seq in Fastq format is read from the hard disk, and the query sequence is processed by weighted syncmer sampling to obtain a number of second sampling k-mer sequences. The operation of weighted syncmer sampling processing is the same as that in S1, and will not be repeated here to save space.

[0050] Next, traverse the seq sequence in each second sampled k-mer sequence, and use the position of the seq sequence in the query sequence as the query start position; the sum of the query start position and the k-mer sequence length as the query end position; at the same time, use the position of the above seq sequence in the reference genome index as the reference start position; and the sum of the reference start position and the k-mer sequence length as the reference end position. The query end position, the corresponding reference end position and the k-mer sequence length k form a candidate position, which is a triple ( ref_end , read_end , k ), ref_end , read_end are the reference end position and the query end position respectively. All query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set, that is, all candidate positions form a candidate position set.

[0051] The specific processing process of S2 can refer to the logical algorithm flow in Table 2.

[0052] Table 2 Logical algorithm flow for obtaining candidate position sets

[0053] .

[0054] S3. The candidate position set is processed by the chain formation method to obtain a chain of aligned gene sequence positions.

[0055] The candidate position set is processed by the chaining method, and the candidate positions are connected into chains to obtain the aligned gene sequence position chain, which can reflect the changes in the alignment results in the continuous gene region.

[0056] The specific operations of the chaining method are as follows: based on the candidate position set, a comparison position diagram is constructed with the horizontal axis being the query end position and the vertical axis being the reference end position, so that the distribution trend of the end position of the sequence in the query sequence and the reference genome sequence during the entire comparison process can be intuitively observed; then the anchor point in the comparison position diagram of the query end position that is the last in the query sequence (the anchor point that is farthest from the origin) is used as the first chaining starting point; within the first chaining region, the anchor point with the highest chaining score between the first chaining starting point and the neighborhood range is used as the anchor point to be chained, and is connected to the first chaining starting point to obtain the first sequence position chain; the anchor point to be chained is used as the starting point of the first sequence position chain to obtain the second chaining starting point, and a connection operation based on the chaining score is performed in the second chaining region to obtain the second sequence position chain; and so on, after the connection of the last anchor point is completed, the comparison gene sequence position chain is obtained, thereby improving the accuracy and intuitiveness of the comparison result.

[0057] The first chained region is a region whose horizontal coordinate does not exceed the horizontal coordinate of the first chained starting point, and whose vertical coordinate does not exceed the vertical coordinate of the first chained starting point; and so on, other chained regions are obtained.

[0058] At the same time, the chain score is obtained by the following formula:

[0059] ,

[0060] f(i) For the i The chain score of the anchor point is j The starting point of the chain i The chaining scores between anchor points are equal. f(j) For the j The chain score of the chain starting point is j The chaining score of the chaining starting point is equal to that of the previous chaining starting point. α(j,i) For the j The starting point of the chain i The number of matching bases between anchor points, β(j,i) For the vacancy penalty,k is the length of the k-mer sequence.

[0061] The method for obtaining the above-mentioned number of matching bases is: obtain the query end position difference and the query end position difference between the j-th chain starting point and the i-th anchor point, and take the minimum value of the query end position difference, the query end position difference and the length of the k-mer sequence as the number of matching bases between the j-th chain starting point and the i-th anchor point.

[0062] Above β(j,i) The gap penalty is greater than 0. The method for obtaining the gap penalty is as follows: respectively obtain the query end position difference and the query end position difference between the j-th chain starting point and the i-th anchor point, and multiply the difference between the query end position difference and the query end position difference by the penalty coefficient as the gap penalty.

[0063] The penalty coefficient can be obtained by the following formula:

[0064] ,

[0065] γ c (l) The length is l The penalty coefficient for sequence position chaining, k is the k-mer sequence length.

[0066] Since the distance from the origin and the chain score may be the same, there may be multiple aligned gene sequence position chains. If there are multiple aligned gene sequence position chains, the aligned gene sequence position chains with a total position chain score greater than the total score threshold are sequentially subjected to the operations in S4 in descending order of the total position chain score. The total position chain score is the chain score of the last linked anchor point in the aligned gene sequence position chain.

[0067] S4. Based on comparing the gene sequence position chain with the query start position, query end position, reference start position, and reference end position, a genome long sequence alignment result is obtained.

[0068] After obtaining the alignment gene sequence position chain, the corresponding alignment result output can be performed according to different requirements of the output results to improve the flexibility of the alignment method.

[0069] If the output format is the PAF alignment format, the query start position, query end position, reference start position, reference end position in the alignment gene sequence position chain, as well as the reference genome sequence and the query sequence can be output as a PAF format alignment result.

[0070] If the output format is SAM alignment format, use the SSE-based double sequence alignment algorithm (call the ksw2 library) to calculate the cigar value of the reference genome sequence and the query sequence, and then output the cigar value, the query start position, query end position, reference start position, reference end position in the alignment gene sequence position chain, as well as the reference genome sequence and the query sequence as SAM format alignment results.

[0071] To verify the effect of the comparison method of this embodiment, the following experiment was conducted.

[0072] Experiment 1: A comparative experiment between the comparison method of this embodiment and the advanced prior art Minimap2 method.

[0073] Experiment 1 was run on a T640 server with an Intel Xeon Gold 6226r processor (2.9 GHZ, 16 cores), 128 GB memory, and an operating system of Ubuntu 20.04. The specific experimental steps are as follows: (1) The simulation sequence tool Pbsim3 was used to simulate a sequencing sequence with a depth of 20 times for chromosome 14 of the human reference genome GRch38 with default parameters (the simulation mode was PacBio RS II CLR, and the error rate was about 10%), generating a total of 116,186 sequences; (2) The simulated sequences were aligned to chromosome 14 of the reference genome GRch38 using the alignment method of this embodiment and the Minimap2 method with default parameters; (3) Pbsim3 could output the correct alignment result (answer position) of the sequencing sequence when simulating the sequence; in order to evaluate the alignment effect of the alignment method of this embodiment and the Minimap2 method, the difference between the output alignment result position and the answer position was statistically analyzed. If the difference was less than 30 (30 bases), it was considered to be a correct alignment, otherwise it was an alignment error. The comparison results are shown in Table 3, where t is the position offset. It can be seen from Table 3 that under the same parameters, the comparison method of this embodiment has more correct comparisons than the Minimap2 method and has a higher accuracy rate.

[0074] Table 3 Simulation data results of the comparison method of this embodiment and the Minimap2 method

[0075] .

[0076] Experiment 2: A comparative experiment between the comparison method of this embodiment and the advanced existing method Winnowmap method.

[0077] Experiment 2 was run on a T640 server with an Intel Xeon Gold 6226r processor (2.9 GHZ, 16 cores), 128 GB memory, and Ubuntu 20.04 operating system. The specific experimental steps are as follows: (1) The simulation sequence tool Pbsim3 was used to simulate the sequencing sequence of chromosome 14 of the human reference genome GRch38 with a depth of 20 times using default parameters (the simulation mode was PacBio RS II CLR, and the error rate was about 10%), generating a total of 116,186 sequences; (2) The simulated sequences were aligned to chromosome 14 of the reference genome GRch38 using the alignment method of this embodiment and the Winnowmap method using default parameters; (3) Pbsim3 could output the correct alignment result (answer position) of the sequencing sequence when simulating the sequence; in order to evaluate the alignment effect of the alignment method of this embodiment and the Winnowmap method, the difference between the output alignment result position and the answer position was statistically analyzed. If the difference was less than 30 (30 bases), it was considered to be a correct alignment, otherwise it was an alignment error.

[0078] The comparison results are shown in Table 4. It can be seen from Table 4 that under the same parameters, the comparison method of this embodiment has more correct comparisons and a higher accuracy rate than the Winnowmap method.

[0079] Table 4 Simulation data results of the comparison method of this embodiment and Winnowmap method

[0080] .

[0081] Experiment 3: Comparative experiment of different sampling methods.

[0082] In order to verify that the weighted syncmer sampling results in the alignment method of this embodiment are distributed more evenly, all sampling methods are sampled and density maps are drawn. The specific steps are as follows: (1) The X chromosome of CHM13 is sampled using the minimizer method, the weighted minimizer method and the weighted syncmer sampling method, respectively, where the k-mer length is 19 and the window length is 17, to ensure that the three methods have similar sampling numbers; (2) The density maps are drawn for the three methods, with the horizontal axis being the X chromosome coordinate of CHM13 and the vertical axis being the sampling density of the region, as shown in FIG. Figure 2 As shown; "minCHX" is a minimizer method, "winCHX" is a weighted minimizer method, and "this embodiment" is a weighted syncmer sampling method; from Figure 2 It can be seen that the three methods achieved nearly uniform density, but the first two had significant differences in the highly repetitive centromere region ( Figure 2A significant density drop was observed in the region around the 6.0e+07 position in the image, while the weighted syncmer sampling method of this embodiment slightly increased the density in this region, indicating that the sampling of this embodiment is more uniform and less affected by highly repetitive sequences.

[0083] This embodiment also provides a genome long sequence alignment system, which is used to implement the above-mentioned genome long sequence alignment method, including:

[0084] A reference genome index generation module is used to obtain a number of first sampling k-mer sequences by weighted syncmer sampling processing of the reference genome sequence; a number of first sampling k-mer sequences and corresponding positions in the reference genome sequence are converted into hash values ​​to obtain a reference genome index;

[0085] The candidate position set generation module is used to query the sequence through weighted syncmer sampling processing to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; all query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set;

[0086] A module for generating a position chain of an aligned gene sequence, which is used to process a candidate position set through a chain formation method to obtain a position chain of an aligned gene sequence;

[0087] The genome long sequence comparison result generation module is used to obtain the genome long sequence comparison result based on the comparison gene sequence position chain with the query start position, query end position, reference start position and reference end position.

[0088] This embodiment also provides a long genome sequence comparison device, including a processor and a memory, wherein the processor implements the above-mentioned long genome sequence comparison method when executing a computer program stored in the memory.

[0089] This embodiment also provides a computer-readable storage medium for storing a computer program, wherein the computer program implements the above-mentioned method for comparing long genomic sequences when executed by a processor.

[0090] The present embodiment provides a method for aligning a long genome sequence. First, a weighted syncmer sampling process is performed on a reference genome sequence and a query sequence respectively, and a more valuable sampled k-mer sequence is screened out for enhancing sensitivity. Then, the query start position and the query end position of a seq sequence in the sampled k-mer sequence in the query sequence are obtained to define the region that needs to be focused on in the query sequence, and the reference start position and the reference end position of the seq sequence in the reference genome index are obtained, so as to provide an accurate alignment interval for the alignment between the query sequence and the reference genome sequence. Next, a chaining process is performed on a candidate position set formed based on the query end position and the corresponding reference end position to obtain an aligned gene sequence position chain, which can reflect the change of the alignment result in a continuous gene region and improve the accuracy and intuitiveness of the alignment result. Finally, according to different requirements for the output result, the information related to the aligned gene sequence position chain is outputted as a corresponding alignment result to improve the flexibility of the alignment method. The method is applied to the alignment of long genome sequences and can improve the alignment accuracy, efficiency and flexibility.

Claims

1. A method for comparing long genome sequences, characterized in that: The following operations are included: S1. The reference genome sequence is processed by weighted syncmer sampling to obtain a number of first sampling k-mer sequences; the first sampling k-mer sequences and their corresponding positions in the reference genome sequence are converted by hash values ​​to obtain a reference genome index; S2, the query sequence is processed by weighted syncmer sampling to obtain a number of second sampling k-mer sequences; The position of the seq sequence in the query sequence in the second sampled k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; All query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set; S3, the candidate position set is processed by a chain formation method to obtain a chain of aligned gene sequence positions; S4. Based on comparing the gene sequence position chain with the query start position, query end position, reference start position, and reference end position, a genome long sequence alignment result is obtained.

2. The method for comparing long genome sequences according to claim 1, characterized in that: The weighted syncmer sampling process in S1 is specifically performed as follows: Obtain the k-mer sequence of each position of the reference genome sequence, and determine whether there is a high-frequency s-mer sequence in each k-mer sequence; If it exists, all s-mer sequences in the current k-mer sequence are converted to hash values ​​and then assigned weights, and the weights of high-frequency s-mer sequences are increased to the corresponding weight threshold to obtain the updated weight of each s-mer sequence; after the updated weight of each s-mer sequence is converted to hash values, they are sorted from large to small according to the hash values; if the sorting position and position offset of the last sorted s-mer sequence are the same, the current k-mer sequence is the first sampled k-mer sequence; If it does not exist, all s-mer sequences in the current k-mer sequence are converted to hash values. If the s-mer sequence with the smallest hash value has the same position and position offset in the current k-mer sequence, the current k-mer sequence is the first sampled k-mer sequence.

3. The method for comparing long genome sequences according to claim 2, characterized in that: The high-frequency s-mer sequence is an s-mer sequence whose occurrence frequency in the reference genome sequence is greater than a frequency threshold; the frequency threshold is based on the total frequency of occurrence of the s-mer sequence in the reference genome sequence.

4. The method for comparing long genome sequences according to claim 1, characterized in that: The operations of the chaining method in S3 are specifically as follows: Based on the candidate position set, a comparison position map is constructed with the horizontal axis being the query end position and the vertical axis being the reference end position; The last query end position in the query sequence, the anchor point in the comparison position graph, is used as the first chain starting point; In the first chaining region, the anchor point with the highest chaining score between the first chaining starting point and the neighborhood is used as the anchor point to be chained, and is connected to the first chaining starting point to obtain a first sequence position chain; Using the anchor point to be chained as the starting point of the first sequence position chain to obtain a second chain starting point, performing a connection operation based on the chain score in the second chain region to obtain a second sequence position chain; By analogy, after the connection of the last anchor point is completed, the aligned gene sequence position chain is obtained.

5. The method for comparing long genome sequences according to claim 4, characterized in that: The chaining score is obtained by the following formula: ， f(i) For the i The chain score of the anchor point is j The starting point of the chain i The chaining scores between anchor points are equal. f(j) For the j The chain score of the chain starting point is j The chaining score of the chaining starting point is equal to that of the previous chaining starting point. α(j,i) For the j The starting point of the chain i The number of matching bases between anchor points, β(j,i) For the vacancy penalty, k is the length of the k-mer sequence.

6. The method for comparing long genome sequences according to claim 5, characterized in that: The method for obtaining the number of matching bases is: Get the query end position difference and query end position difference between the j-th chain starting point and the i-th anchor point, and take the minimum value of the query end position difference, the query end position difference and the length of the k-mer sequence as the number of matching bases between the j-th chain starting point and the i-th anchor point.

7. The method for comparing long genome sequences according to claim 1, characterized in that: If there are multiple aligned gene sequence position chains, the aligned gene sequence position chains whose total position chain scores are greater than the total score threshold are sequentially subjected to the operations in S4 in descending order of the total position chain scores; the total position chain score is the chain score of the last chained anchor point in the aligned gene sequence position chain.

8. A genome long sequence alignment system, used to implement the genome long sequence alignment method according to claim 1, characterized in that: include: A reference genome index generation module is used to obtain a number of first sampling k-mer sequences by weighted syncmer sampling processing of the reference genome sequence; a number of first sampling k-mer sequences and corresponding positions in the reference genome sequence are converted into hash values ​​to obtain a reference genome index; The candidate position set generation module is used to query the sequence through weighted syncmer sampling processing to obtain several second sampling k-mer sequences; the position of the seq sequence in the query sequence in the second sampling k-mer sequence is used as the query start position; the sum of the query start position and the k-mer sequence length is used as the query end position; the position of the seq sequence in the reference genome index is used as the reference start position; the sum of the reference start position and the k-mer sequence length is used as the reference end position; all query end positions and corresponding reference end positions, as well as the k-mer sequence length, form a candidate position set; A module for generating a position chain of an aligned gene sequence, used for processing the candidate position set by a chain formation method to obtain a position chain of an aligned gene sequence; The genome long sequence comparison result generation module is used to obtain the genome long sequence comparison result based on the comparison gene sequence position chain with the query start position, query end position, reference start position and reference end position.

9. A genome long sequence comparison device, characterized in that: The method comprises a processor and a memory, wherein the processor implements the long genome sequence comparison method as described in any one of claims 1 to 7 when executing the computer program stored in the memory.

10. A computer-readable storage medium, characterized in that: Used to store a computer program, wherein when the computer program is executed by a processor, the method for aligning long genomic sequences as described in any one of claims 1 to 7 is implemented.

Images