A method for analyzing genomic sequencing data based on next-generation sequencing technology

Abstract

The present application relates to the technical field of data processing, and especially relates to a genome sequencing data analysis method based on second-generation sequencing technology, which obtains a base sequence to be compressed; sets a k value range of k-mer, obtains a k-mer substring of any k value in the base sequence for any k value in the k value range, forms a target set, and constructs a frequency distribution diagram of the target set; according to the data distribution characteristics in the frequency distribution diagram and the position distribution characteristics of each k-mer substring corresponding to a homozygous peak in the frequency distribution diagram in the base sequence, obtains a repetition evaluation value of any k value; obtains the repetition evaluation value of each k value in the k value range, obtains an optimal k value according to the repetition evaluation value of each k value in the k value range, takes the optimal k value as the length of a repeated string in a LZW dictionary compression algorithm, and compresses the base sequence, thereby improving the compression efficiency of the base sequence.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a method for analyzing genome sequencing data based on next-generation sequencing technology. Background Technology

[0002] Next-generation sequencing (NGS) technology, with its high throughput, speed, and low cost, has become the mainstream gene sequencing technology, finding wide application in life science research, medical diagnostics, drug development, and genetics. However, NGS generates extremely large amounts of data; sequencing a single genome can produce hundreds of gigabytes or even terabytes of data, placing immense pressure on data storage and transmission. Therefore, there is an urgent need for a compression algorithm with high compression efficiency suitable for NGS data to improve its storage and transmission capabilities. The LZW dictionary compression algorithm, as a mainstream compression algorithm, is often chosen as the compression algorithm for NGS technology.

[0003] The LZW dictionary compression algorithm uses a sliding window to search for duplicate strings in the data and constructs a dictionary table for these duplicate strings. It then represents the repeated strings with short codes, achieving lossless data compression. Second-generation sequencing data mainly consists of the bases A, T, C, and G, with a relatively random distribution. For example, the string TCAG may contain many duplicates, as may TCA and CAG within TCAG, and even TCAGT, which extends one position after TCAG, may contain even more duplicates. Therefore, when using the LZW dictionary compression algorithm to compress sequencing data and find duplicate strings, it is difficult to determine the string length, thus making it impossible to determine a suitable sliding window length, resulting in unstable compression efficiency.

[0004] Therefore, how to obtain the length of repeated strings in the LZW dictionary compression algorithm through data analysis of second-generation sequencing data, and improve the compression efficiency of second-generation sequencing data, has become an urgent problem to be solved. Summary of the Invention

[0005] In view of this, embodiments of the present invention provide a genome sequencing data analysis method based on second-generation sequencing technology to solve the problem of how to obtain the length of repeated strings in the LZW dictionary compression algorithm by analyzing second-generation sequencing data, thereby improving the compression efficiency of second-generation sequencing data.

[0006] This invention provides a method for analyzing genome sequencing data based on next-generation sequencing technology, the method comprising the following steps:

[0007] The genome was sequenced using next-generation sequencing technology to obtain the base sequence to be compressed;

[0008] Define a range of k values ​​for k-mer. For any k value within the range, obtain at least one k-mer substring of length equal to any k value in the base sequence to form a target set. Construct a frequency distribution map of the target set, where the horizontal axis of the frequency distribution map represents the frequency of occurrence of the k-mer substring, and the vertical axis represents the number of k-mer substrings corresponding to the frequency of occurrence.

[0009] Based on the data distribution characteristics in the frequency distribution map, obtain the frequency distribution characteristic value of any k value, obtain the homozygous peak in the frequency distribution map, obtain the spatial distribution characteristic value of any k value based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, and obtain the repetition evaluation value of any k value based on the frequency distribution characteristic value and spatial distribution characteristic value of any k value.

[0010] Obtain the repetition evaluation value for each value of k within the range of k values, obtain the optimal k value based on the repetition evaluation value for each value of k within the range of k values, and use the optimal k value as the length of the repeated string in the LZW dictionary compression algorithm to compress the base sequence.

[0011] Preferably, obtaining the frequency distribution feature value of any k value based on the data distribution characteristics in the frequency distribution map includes:

[0012] Obtain peak data points and trough data points in the frequency distribution map. For any peak data point, obtain the trough data point that is to the left of the peak data point and record it as the left adjacent trough point. Obtain the trough data point that is to the right of the peak data point and record it as the right adjacent trough point.

[0013] The data points between the left adjacent trough point and the right adjacent trough point, as well as the left adjacent trough point and the right adjacent trough point, form the local fluctuation range of any peak data point.

[0014] The sum of the ordinates of all data points within the local fluctuation interval is obtained to get the total number of the local fluctuation interval. The difference between the abscissas of the right adjacent trough point and the left adjacent trough point is obtained to get the frequency span value of the local fluctuation interval.

[0015] The total quantity value and frequency span value of the local fluctuation interval corresponding to each peak data point in the frequency distribution map are obtained respectively. Based on the difference in the total quantity value and frequency span value of the local fluctuation interval corresponding to each peak data point in the frequency distribution map, the frequency distribution feature value of any k value is obtained.

[0016] Preferably, the step of obtaining the frequency distribution feature value of any k value based on the difference in the total quantity value and the difference in the frequency span value of the local fluctuation interval corresponding to each peak data point in the frequency distribution map includes:

[0017] The peak data points in the frequency distribution map are divided into one homozygous peak and at least one heterozygous peak;

[0018] The total number of local fluctuation intervals corresponding to all heterozygous peaks is accumulated to obtain the total accumulated value. The ratio of the total number of local fluctuation intervals corresponding to the homozygous peak to the total accumulated value is obtained to obtain the quantitative characteristic ratio of any k value.

[0019] The frequency span values ​​of the local fluctuation intervals corresponding to all heterozygous peaks are accumulated to obtain the frequency span accumulation value. The ratio of the frequency span value of the local fluctuation interval corresponding to the homozygous peak to the frequency span accumulation value is obtained to obtain the frequency characteristic ratio of any k value.

[0020] The frequency distribution characteristic value of any k value is obtained by multiplying the quantity characteristic ratio and the frequency characteristic ratio.

[0021] Preferably, the step of obtaining the spatial distribution feature value of any k value based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence includes:

[0022] All k-mer substrings corresponding to the homozygous peak are recorded as high-repetition substrings. For any high-repetition substring, if any high-repetition substring has overlapping data with other high-repetition substrings, then any high-repetition substring is recorded as an overlapping substring.

[0023] Based on the number of overlapping substrings in all highly repeating substrings, the proportion of overlapping substrings in all highly repeating substrings is obtained, and the spacing distance between each two adjacent highly repeating substrings in the base sequence is obtained, and the spacing distance variance is obtained accordingly.

[0024] The spatial distribution characteristic value of any k value is obtained by multiplying the reciprocal of the proportion of overlapping substrings with the reciprocal of the variance of the interval distance.

[0025] Preferably, the step of obtaining the spacing distance between each pair of adjacent highly repetitive substrings in the base sequence includes:

[0026] For any two adjacent highly repetitive substrings in the base sequence, obtain the position number of the first base in the base sequence of the two adjacent highly repetitive substrings, calculate the absolute value of the difference between the position numbers of the first base in the base sequence of the two adjacent highly repetitive substrings, and obtain the interval distance of the two adjacent highly repetitive substrings in the base sequence.

[0027] Preferably, the step of obtaining the repeatability evaluation value of any k value based on the frequency distribution characteristic value and spatial distribution characteristic value of any k value includes:

[0028] The frequency distribution feature value is linearly normalized to obtain a frequency distribution feature normalized value. The repeatability evaluation value is linearly normalized to obtain a repeatability evaluation normalized value. The repeatability evaluation value for any value of k is obtained by adding the frequency distribution feature normalized value and the repeatability evaluation normalized value.

[0029] Preferably, the step of obtaining the optimal k value based on the repetition evaluation value of each k value within the range of k values ​​includes:

[0030] Among the repeatability evaluation values ​​for each k value within the range of k values, the maximum repeatability evaluation value is selected, and the k value corresponding to the maximum repeatability evaluation value is taken as the optimal k value.

[0031] The beneficial effects of the embodiments of the present invention compared with the prior art are as follows:

[0032] This invention uses next-generation sequencing technology to sequence the genome, obtaining a base sequence to be compressed; it sets a range of k values ​​for k-mers, and for any k value within this range, it extracts at least one k-mer substring of length equal to that k value from the base sequence, forming a target set; it then constructs a frequency distribution map of the target set, where the horizontal axis of the frequency distribution map represents the frequency of occurrence of the k-mer substring, and the vertical axis represents the number of k-mer substrings corresponding to each frequency; based on the data distribution characteristics in the frequency distribution map, it obtains the frequency distribution characteristic value for any k value, thus obtaining... The homozygous peaks in the frequency distribution map are selected. Based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, the spatial distribution characteristic value of any k value is obtained. Based on the frequency distribution characteristic value and spatial distribution characteristic value of any k value, the repetition evaluation value of any k value is obtained. The repetition evaluation value of each k value within the range of k values ​​is obtained. Based on the repetition evaluation value of each k value within the range of k values, the optimal k value is obtained. The optimal k value is used as the length of the repeating string in the LZW dictionary compression algorithm to compress the base sequence. Specifically, by analyzing the data distribution characteristics of any k value within the range of k values ​​in the frequency distribution map, and the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the frequency distribution map in the base sequence, the overall repetition of the k-mer substring under any k value is evaluated, obtaining the repetition evaluation value of any k value. The optimal k value is then obtained as the most suitable length of the repeating string in the LZW dictionary compression algorithm, improving the compression efficiency of the base sequence and thus enhancing the storage and transmission capabilities of second-generation sequencing data. Attached Figure Description

[0033] To more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0034] Figure 1 This is a flowchart of a genome sequencing data analysis method based on second-generation sequencing technology provided in Embodiment 1 of the present invention;

[0035] Figure 2 This is an example diagram of second-generation sequencing data provided in Embodiment 1 of the present invention. Detailed Implementation

[0036] Embodiments of this disclosure are described in detail below, with examples of these embodiments illustrated in the accompanying drawings. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain this disclosure, and should not be construed as limiting it.

[0037] It should be noted that the terms "first," "second," etc., used in this disclosure and the accompanying drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this disclosure described herein can be implemented in orders other than those illustrated or described herein. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this disclosure. Rather, they are merely examples of apparatuses and methods consistent with some aspects of this disclosure.

[0038] To illustrate the technical solution of the present invention, specific embodiments are described below.

[0039] See Figure 1 This is a flowchart of a genome sequencing data analysis method based on second-generation sequencing technology provided in Embodiment 1 of the present invention. Figure 1 As shown, the method may include:

[0040] Step S101: Use next-generation sequencing technology to sequence the genome to obtain the base sequence to be compressed.

[0041] Next-generation sequencing (NGS) technology is used to sequence the genome, obtaining raw data in FastQ format. This raw data is then parsed using the FastQC tool to obtain next-generation sequencing data, such as... Figure 2 As shown, the first row is the sequence attributes, the second row is the base sequence to be compressed, containing four types of base data: A, T, C, and G. The third row is the separator "+", and the fourth row is the quality score encoding.

[0042] Because next-generation sequencing (NGS) technology generates extremely large amounts of data, with single genome sequencing generating hundreds of GB or even TB of data, data storage and transmission face enormous data carrying capacity pressure. Therefore, there is an urgent need for a compression algorithm with strong compression efficiency suitable for NGS data to improve the storage and transmission capabilities of NGS data.

[0043] The LZW dictionary compression algorithm, as a mainstream compression algorithm, is often chosen for compression in next-generation sequencing technology. It searches for duplicate strings in the data by setting a sliding window and constructing a dictionary table for these duplicate strings, representing them with short codes to achieve lossless data compression. However, the base distribution in next-generation sequencing data is highly random. For example, the string TCAG may have many duplicates, as may TCA and CAG within TCAG, and even TCAGT, which extends one position after TCAG, may have even more duplicates. Therefore, when using the LZW dictionary compression algorithm to find duplicate strings in sequencing data, it is difficult to determine the string length, and consequently, the appropriate sliding window length, leading to unstable compression efficiency.

[0044] Therefore, in this embodiment, by analyzing the degree of repetition of strings of different lengths in the base sequence to be compressed, the appropriate length of the repeating strings in the LZW dictionary compression algorithm is obtained, thereby improving the compression efficiency of second-generation sequencing data.

[0045] Step S102: Set the range of k values ​​for k-mer. For any k value within the range of k values, obtain at least one k-mer substring of length equal to any k value in the base sequence to form a target set. Construct a frequency distribution map of the target set, wherein the horizontal axis of the frequency distribution map is the frequency of occurrence of the k-mer substring, and the vertical axis is the number of k-mer substrings corresponding to the frequency of occurrence.

[0046] To analyze the repetition of strings of different lengths in the base sequence to be compressed, it is necessary to first obtain strings of different lengths in the base sequence. In genome sequencing data analysis, k-mers are often used to segment base data, dividing the sequencing data into smaller, overlapping fragments. For example, when k=3, the k-mers of the sequence ATCGTA are ATC, TCG, CGT, and GTA. Therefore, in this embodiment, by setting the range of k values ​​for the k-mer, strings of different lengths in the base sequence are obtained.

[0047] Since the minimum length of a k-mer is 3, and the maximum value of k in a k-mer is normally around 21 in sequencing technology, this embodiment sets the range of k values ​​for k-mers to [3, 21]. This is not a restriction and can be set according to the specific implementation scenario. For any k value within the range, a sliding window of length k is established in the base sequence. Starting from the first base of the base sequence, it moves one base to the right each time. If the number of bases in the base sequence is N, then N-k+1 k-mer substrings of length k are obtained for any k value. Similarly, N-k+1 k-mer substrings of length k can be obtained for each k value within the range.

[0048] In a base sequence, the more frequently a k-mer substring appears, the greater the degree of repetition of the k-mer substring in the base sequence. Therefore, we can form a target set by taking N-k+1 k-mer substrings of length k for any value of k, and denote the k value corresponding to the target set as the current k value. By analyzing the frequency distribution characteristics of the k-mer substrings in the target set, we can analyze the degree of repetition of the k-mer substrings under the current k value.

[0049] To better analyze the frequency distribution characteristics of k-mer substrings within the target set, obtain the occurrence frequency of each k-mer substring in the target set, and the number of k-mer substrings corresponding to each occurrence frequency, a frequency distribution map of the target set is constructed. In the frequency distribution map of the target set, the horizontal axis represents the occurrence frequency of k-mer substrings, and the vertical axis represents the number of k-mer substrings corresponding to each occurrence frequency. Each point in the frequency distribution map represents the total number of k-mer substrings with the corresponding occurrence frequency. For example, (30, 200) means that the total number of k-mer substrings with a frequency of 30 is 200.

[0050] Thus, we obtain the frequency distribution map of the target set corresponding to the current value of k.

[0051] Step S103: Based on the data distribution characteristics in the frequency distribution map, obtain the frequency distribution characteristic value of any k value, obtain the homozygous peak in the frequency distribution map, obtain the spatial distribution characteristic value of any k value based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, and obtain the repetition evaluation value of any k value based on the frequency distribution characteristic value and the spatial distribution characteristic value of any k value.

[0052] Under the theoretical condition of no sequencing error, the frequency distribution map of the target set usually conforms to a Poisson distribution, that is, there is a main peak, and it may contain smaller peaks caused by other features in the genome. Here, the main peak is called the homozygous peak, and the other peaks are called heterozygous peaks, which are located before the main peak. Normally, the frequency distribution map also contains error peaks at the beginning and repetitive sequence peaks after the homozygous peak, but these will be removed as errors when the base sequence is collected. Therefore, the frequency distribution map in this embodiment only includes homozygous peaks and heterozygous peaks.

[0053] The AMPD peak lookup algorithm is used to obtain peak and trough data points in the frequency distribution map. The AMPD peak lookup algorithm is existing technology and will not be described in detail here. For any peak data point, the trough data point to the left of the peak data point is recorded as the left adjacent trough point, and the trough data point to the right of the peak data point is recorded as the right adjacent trough point. The data points between the left and right adjacent trough points, as well as the left and right adjacent trough points, form the local fluctuation interval of any peak data point. Each peak data point is the highest point of its local fluctuation interval, and the local fluctuation interval represents the clustering of a group of k-mer substrings with similar occurrence frequencies.

[0054] Because the homozygous peak is located towards the end of the frequency distribution map, the local fluctuation interval corresponding to the homozygous peak contains k-mer substrings with higher frequencies, while the earlier heterozygous peak contains k-mer substrings with lower frequencies. In the LZW dictionary compression algorithm, the higher the frequency of k-mer substrings, the higher the compression efficiency obtained by using the dictionary during the compression process. Therefore, when the proportion of the local fluctuation interval corresponding to the homozygous peak in the frequency distribution map is larger than that of the heterozygous peak, it indicates that more k-mer substrings are clustered in the homozygous peak at the current k value. This means that there are more k-mer substrings with higher frequencies among all k-mer substrings at the current k value, resulting in a higher frequency distribution of k-mer substrings in the target set and higher compression efficiency obtained by using the dictionary during the compression process.

[0055] Therefore, the frequency distribution feature value of the current k value can be obtained by analyzing the difference between the local fluctuation intervals corresponding to homozygous peaks and those corresponding to all heterozygous peaks in the frequency distribution map, as well as the difference in the number of k-mer substrings corresponding to homozygous peaks and those corresponding to all heterozygous peaks. Then, the repetition degree of the k-mer substrings under the current k value can be analyzed using this frequency distribution feature value. The method for obtaining the frequency distribution feature value of the current k value is as follows:

[0056] (1) For any peak data point, obtain the sum of the vertical coordinates of all data points in the local fluctuation interval of any peak data point to obtain the total number of the local fluctuation interval, obtain the difference between the horizontal coordinates of the right adjacent valley point and the left adjacent valley point to obtain the frequency span value of the local fluctuation interval.

[0057] In one embodiment, taking the i-th peak data point in the frequency distribution map as an example, the right adjacent trough point and the left adjacent trough point of the i-th peak data point are obtained to obtain the local fluctuation region of the i-th peak data point. The formula for calculating the total number of local fluctuation intervals and the frequency span value of the i-th peak data point is as follows:

[0058]

[0059] Where F is the total number of local fluctuation intervals for the i-th peak data point; Let be the ordinate of the j-th data point within the local fluctuation interval of the i-th peak data point; n is the number of data points within the local fluctuation interval of the i-th peak data point; P is the frequency span value of the local fluctuation interval of the i-th peak data point. Let x be the x-coordinate of the right adjacent trough point of the i-th peak data point; Let x be the x-coordinate of the left adjacent trough point of the i-th peak data point.

[0060] (2) Obtain the total quantity value and frequency span value of the local fluctuation interval corresponding to each peak data point in the frequency distribution map respectively. Based on the difference in the total quantity value and frequency span value of the local fluctuation interval corresponding to each peak data point in the frequency distribution map, obtain the frequency distribution feature value of the current k value.

[0061] Specifically, the peak data points in the frequency distribution map are divided into one homozygous peak and at least one heterozygous peak;

[0062] The total number of local fluctuation intervals corresponding to all heterozygous peaks is accumulated to obtain the total accumulated value. The ratio of the total number of local fluctuation intervals corresponding to the homozygous peak to the total accumulated value is obtained to obtain the quantitative characteristic ratio of the current value of k.

[0063] The frequency span values ​​of the local fluctuation intervals corresponding to all heterozygous peaks are accumulated to obtain the frequency span accumulation value. The ratio of the frequency span value of the local fluctuation interval corresponding to the homozygous peak to the frequency span accumulation value is obtained to obtain the frequency characteristic ratio of the current k value.

[0064] The frequency distribution characteristic value of the current k value is obtained by multiplying the quantity characteristic ratio and the frequency characteristic ratio.

[0065] In one embodiment, the formula for calculating the frequency distribution characteristic value of the current k value is:

[0066]

[0067] Where W is the frequency distribution characteristic value of the current value of k; This represents the total number of local fluctuation regions corresponding to homozygous peaks in the frequency distribution diagram. This represents the frequency span value of the local fluctuation region corresponding to the homozygous peak in the frequency distribution diagram. This represents the total number of local fluctuation regions corresponding to the z-th heterozygous peak in the frequency distribution diagram. is the frequency span value of the local fluctuation region corresponding to the z-th heterozygous peak in the frequency distribution map; m is the number of heterozygous peaks in the frequency distribution map.

[0068] It should be noted that, This represents the ratio of the number of values ​​for the current k value. The larger the value, the more k-mer substrings are corresponding to the homozygous peak; This represents the frequency characteristic ratio of the current value of k. The larger the value, the greater the distribution range of k-mer substring frequency in the homozygous peak compared to the heterozygous peak, indicating that more frequent k-mer substrings are clustered in the homozygous peak; when and When both are large, it means that there are more k-mer substrings with higher occurrence frequencies among all k-mer substrings under the current k value. The greater the repetition of k-mer substrings in the target set, the larger the frequency distribution feature value of the current k value, and the higher the compression efficiency obtained by using the dictionary during the compression process.

[0069] However, frequency distribution features are only analyzed based on the frequency of occurrence of k-mer substrings. The distribution of k-mer substrings in the base sequence may show some regional clustering, leading to overlapping of k-mer substrings. For example, when k=3, TCT and CTG in ATCTGC are both k-mer substrings within homozygous peaks, but these two substrings overlap in ATCTGC. When traversing strings using a dictionary, the step size of the sliding window is equal to the length of the sliding window, meaning the sliding windows are actually non-overlapping. When actually selecting strings for the dictionary using the sliding window, only one k-mer substring can be selected to generate the dictionary. Therefore, although overlapping k-mer substrings can provide a higher degree of repetition, it is still necessary to select k-mer substrings with a lower overlap rate with other k-mer substrings to be included in the compressed dictionary during actual data traversal.

[0070] Therefore, it is necessary to consider the distribution characteristics of highly repetitive k-mer substrings in the base sequence. The more uniform the distribution of highly repetitive k-mer substrings in the base sequence, the better the dictionary conversion effect of the corresponding base sequence. Therefore, it is necessary to obtain the spatial distribution characteristic value of the current k value based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, and combine the frequency distribution characteristics and positional distribution characteristics of the current k value to obtain the repetition evaluation value of the current k value. Then, the repetition degree of k-mer substrings under the current k value can be analyzed through the repetition evaluation value.

[0071] The method for obtaining the spatial distribution feature value of the current k value is as follows, based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence:

[0072] All k-mer substrings corresponding to the homozygous peak are recorded as high-repetition substrings. For any high-repetition substring, if any high-repetition substring has overlapping data with other high-repetition substrings, then any high-repetition substring is recorded as an overlapping substring.

[0073] For any two adjacent highly repeating substrings in the base sequence, obtain the position number of the first base in the base sequence of the two adjacent highly repeating substrings, calculate the absolute value of the difference between the position numbers of the first base in the base sequence of the two adjacent highly repeating substrings, and obtain the interval distance of the two adjacent highly repeating substrings in the base sequence.

[0074] Based on the number of overlapping substrings in all highly repeating substrings, the proportion of overlapping substrings in all highly repeating substrings is obtained, and the spacing distance between each two adjacent highly repeating substrings in the base sequence is obtained, and the spacing distance variance is obtained accordingly.

[0075] The spatial distribution feature value of the current k value is obtained by multiplying the reciprocal of the proportion of overlapping substrings with the reciprocal of the variance of the interval distance.

[0076] In one embodiment, the formula for calculating the spatial distribution characteristic value of the current k value is:

[0077]

[0078] Where E is the spatial distribution characteristic value of the current value of k; The variance of the interval distance; denoted as , where is the spacing between adjacent highly repetitive substrings in the homozygous peak within the base sequence; Q is the number of highly repetitive substrings in the homozygous peak; and q is the number of overlapping substrings in the homozygous peak.

[0079] It should be noted that, The variance of the interval distance. The smaller the value, the more uniform the distribution of highly repetitive substrings in the homozygous peak, and the larger the spatial distribution characteristic value of the current k value. It is the reciprocal of the proportion of overlapping substrings in homozygous peaks. The larger the value of k, the smaller the proportion of overlapping substrings in the homozygous peak, the more uniform the distribution of highly repetitive substrings in the homozygous peak, and the larger the spatial distribution characteristic value of the current k value.

[0080] Furthermore, combining the frequency distribution characteristics and location distribution characteristics of the current value of k, the method for obtaining the repetition evaluation value of the current value of k is as follows:

[0081] The frequency distribution feature value is linearly normalized to obtain a frequency distribution feature normalized value. The repeatability evaluation value is linearly normalized to obtain a repeatability evaluation normalized value. The repeatability evaluation value for the current k value is obtained by adding the frequency distribution feature normalized value and the repeatability evaluation normalized value.

[0082] In one implementation, the formula for calculating the repetition evaluation value of the current k value is:

[0083]

[0084] Where S is the repetition evaluation value of the current value of k; W is the frequency distribution characteristic value of the current value of k; E is the spatial distribution characteristic value of the current value of k; and norm() is the linear normalization function.

[0085] It should be noted that the larger W is, the more frequently occurring k-mer substrings exist among all k-mer substrings at the current k value. The greater the repetition of k-mer substrings in the target set, the higher the repetition evaluation value of the current k value, and the higher the compression efficiency obtained by using the dictionary during the compression process. The larger E is, the more uniformly the high-repetition substrings in the homozygous peaks in the frequency distribution diagram of the target set are distributed in the base sequence. The higher the repetition evaluation value of the current k value, and the higher the compression efficiency obtained by using the dictionary during the compression process.

[0086] Thus, the repetition evaluation value for the current value of k is obtained.

[0087] Step S104: Obtain the repetition evaluation value for each k value within the range of k values, obtain the optimal k value based on the repetition evaluation value for each k value within the range of k values, and use the optimal k value as the length of the repeated string in the LZW dictionary compression algorithm to compress the base sequence.

[0088] Based on the method for obtaining the repetition evaluation value of the current k value, the repetition evaluation value for each k value within the range of k values ​​is obtained. Since a larger repetition evaluation value indicates a greater degree of repetition of the k-mer substring at the corresponding k value, and a more uniform distribution of highly repetitive substrings in the base sequence at the corresponding k value, the maximum repetition evaluation value is selected from the repetition evaluation values ​​of each k value within the range of k values. The k value corresponding to the maximum repetition evaluation value is taken as the optimal k value, and the optimal k value is used as the length of the repeated string in the LZW dictionary compression algorithm to construct a dictionary. Since the dictionary construction method in the LZW dictionary compression algorithm is existing technology, it will be briefly described here. Description: Starting from the beginning of the base sequence, the string is truncated using the optimal k value as the sliding window length. Then, the sliding window moves to the next string with the optimal k value as the step size. When the truncated string is the same as a previously truncated string, the string is recorded in the dictionary. The first string recorded is recorded in the dictionary as 00000001 in eight-bit binary ASCII format, the second string is recorded as 00000010, the third string is recorded as 00000011, and so on, incrementing by one using binary counting until the sliding window ends and all duplicate data is recorded. This completes the dictionary construction and the compression of the base sequence.

[0089] In summary, this embodiment uses next-generation sequencing technology to sequence the genome, obtaining the base sequence to be compressed; it sets a range of k values ​​for k-mers, and for any k value within the range, it extracts at least one k-mer substring of length equal to that k value from the base sequence, forming a target set; it then constructs a frequency distribution map of the target set, where the horizontal axis of the frequency distribution map represents the frequency of occurrence of the k-mer substring, and the vertical axis represents the number of k-mer substrings corresponding to each frequency; based on the data distribution characteristics in the frequency distribution map, it obtains the frequency distribution characteristics for any k value. The algorithm obtains homozygous peaks from the frequency distribution map, and based on the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, obtains the spatial distribution characteristic value of any k value. Based on the frequency and spatial distribution characteristic values ​​of any k value, it obtains the repetition evaluation value of any k value. It then obtains the repetition evaluation value for each k value within the range of k values, and obtains the optimal k value based on these values. This optimal k value is used as the length of the repeating string in the LZW dictionary compression algorithm to compress the base sequence. By analyzing the data distribution characteristics of any k value within the range of k values ​​in the frequency distribution map, and the positional distribution characteristics of each k-mer substring corresponding to the homozygous peak in the frequency distribution map within the base sequence, the overall repetition of the k-mer substring under any k value is evaluated, resulting in a repetition evaluation value for any k value. This allows for the determination of the optimal k value as the most suitable length of the repeating string in the LZW dictionary compression algorithm, improving the compression efficiency of the base sequence and thus enhancing the storage and transmission capabilities of second-generation sequencing data.

[0090] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A method for analyzing genomic sequencing data based on a next-generation sequencing technology, characterized by, The method comprises: sequencing the genome using a second-generation sequencing technology to obtain a base sequence to be compressed; setting a k-mer k value range, for any k value in the k value range, obtaining at least one k-mer substring with the length of the any k value in the base sequence to form a target set, and constructing a frequency distribution graph of the target set, wherein the abscissa of the frequency distribution graph is the occurrence frequency of the k-mer substring, and the ordinate is the number of the k-mer substring corresponding to the occurrence frequency; obtaining a frequency distribution characteristic value of the any k value according to the data distribution characteristics in the frequency distribution graph, obtaining a homozygous peak in the frequency distribution graph, obtaining a spatial distribution characteristic value of the any k value according to the position distribution characteristics of each k-mer substring corresponding to the homozygous peak in the base sequence, and obtaining a repetition degree evaluation value of the any k value according to the frequency distribution characteristic value and the spatial distribution characteristic value of the any k value; obtaining the repetition degree evaluation value of each k value in the k value range respectively, obtaining an optimal k value according to the repetition degree evaluation value of each k value in the k value range, taking the optimal k value as the length of the repeated string in the LZW dictionary compression algorithm, and compressing the base sequence.

2. The method of claim 1, wherein the method is based on a next-generation sequencing technology. The method comprises: obtaining wave peak data points and wave valley data points in the frequency distribution graph, for any wave peak data point, obtaining a left adjacent wave valley data point of the wave peak data point, and obtaining a right adjacent wave valley data point of the wave peak data point; the data points between the left adjacent wave valley data point and the right adjacent wave valley data point and the left adjacent wave valley data point and the right adjacent wave valley data point form a local fluctuation interval of the any wave peak data point; obtaining the sum of the ordinates of all data points in the local fluctuation interval to obtain a total quantity value of the local fluctuation interval, and obtaining the difference between the abscissa of the right adjacent wave valley data point and the left adjacent wave valley data point to obtain a frequency span value of the local fluctuation interval; obtaining the total quantity value and the frequency span value of the local fluctuation interval corresponding to each wave peak data point in the frequency distribution graph respectively, and obtaining the frequency distribution characteristic value of the any k value according to the difference between the total quantity value and the frequency span value of the local fluctuation interval corresponding to each wave peak data point in the frequency distribution graph.

3. The method of claim 2, wherein the method is a method of analyzing genomic sequencing data based on a next-generation sequencing technology. The method comprises: dividing the wave peak data points in the frequency distribution graph into a homozygous peak and at least one heterozygous peak; accumulating the total quantity values of the local fluctuation intervals corresponding to all heterozygous peaks to obtain a total quantity accumulation value, obtaining a quantity characteristic ratio value of the any k value by obtaining the ratio of the total quantity value of the local fluctuation interval corresponding to the homozygous peak to the total quantity accumulation value; and Accumulate the frequency span values of the local fluctuation intervals corresponding to all hybrid peaks to obtain a frequency span accumulation value, obtain a ratio of the frequency span value of the local fluctuation interval corresponding to the homozygous peak to the frequency span accumulation value, and obtain a frequency feature ratio value of any k value; Obtain a frequency distribution feature value of any k value according to the product of the quantity feature ratio value and the frequency feature ratio value.

4. The method of claim 1, wherein the method is based on a next-generation sequencing technology. The space distribution feature value of any k value is obtained according to the position distribution feature of each k-mer substring corresponding to the homozygous peak in the base sequence, including: All k-mer substrings corresponding to the homozygous peak are recorded as high repetitive substrings, and for any high repetitive substring, if the any high repetitive substring overlaps with a high repetitive substring other than the any high repetitive substring, the any high repetitive substring is recorded as an overlapping substring; Obtain an overlapping substring proportion in all high repetitive substrings according to the number of overlapping substrings in all high repetitive substrings, obtain the interval distance of each two adjacent high repetitive substrings in the base sequence respectively, and obtain the interval distance variance correspondingly; Obtain the product of the reciprocal of the overlapping substring proportion and the reciprocal of the interval distance variance to obtain the space distribution feature value of any k value.

5. The method of claim 4, wherein the method is a method of analyzing genomic sequencing data based on a next-generation sequencing technology. The interval distance of each two adjacent high repetitive substrings in the base sequence is obtained respectively, including: For any two adjacent high repetitive substrings in the base sequence, the position number of the first base in the any two adjacent high repetitive substrings in the base sequence is obtained respectively, the absolute value of the difference of the position number of the first base in the any two adjacent high repetitive substrings in the base sequence is calculated, and the interval distance of the any two adjacent high repetitive substrings in the base sequence is obtained.

6. The method of claim 1, wherein the method is based on a next-generation sequencing technology. The repeat degree evaluation value of any k value is obtained according to the frequency distribution feature value and the space distribution feature value of any k value, including: The frequency distribution feature value is linearly normalized to obtain a frequency distribution feature normalized value, the repeat degree evaluation value is linearly normalized to obtain a repeat degree evaluation normalized value, and the repeat degree evaluation value of any k value is obtained according to the addition result between the frequency distribution feature normalized value and the repeat degree evaluation normalized value.

7. The method of claim 1, wherein the method is based on a next-generation sequencing technology. The best k value is obtained according to the repeat degree evaluation value of each k value in the k value range, including: The maximum repeat degree evaluation value is selected from the repeat degree evaluation values of each k value in the k value range, and the k value corresponding to the maximum repeat degree evaluation value is taken as the best k value.

Images