Genome sequence compression method and system

Abstract

The application relates to the technical field of bioinformatics data processing, in particular to a genomic sequence compression method and system. The method comprises the following steps: acquiring genomic sequencing data; determining difference sites by aligning the sequencing data with a reference genome; distinguishing the difference sites as sequencing errors or real variations through a time series difference neural network model; performing differential compression coding according to the recognition result; and constructing a variation detection friendly format to support fast variation query through a multi-level index structure and a variation metadata table. The method provided by the application can accurately distinguish sequencing errors and real variations through a time series difference neural network model, and important biological variation information can be protected while improving the compression efficiency by adopting a differential compression coding strategy.

Description

Technical Field

[0001] This application relates to the field of bioinformatics data processing technology, and in particular to a method and system for genomic sequence compression. Background Technology

[0002] With the development of high-throughput sequencing technology, the demand for large-scale population genomics research in precision medicine and genetics research is increasing. It requires efficient data compression and support for rapid variation frequency statistics, while preserving complete variation information to analyze genetic diversity.

[0003] Traditional genomic data compression techniques mainly include lossless compression (such as CRAM and BAM) and lossy compression. CRAM and BAM formats utilize sequence alignment and quality scores to compress data using standard algorithms, supporting variant detection and frequency analysis. Lossy compression further reduces data volume through simple quality score filtering, making it suitable for storage-constrained scenarios. However, in large-scale population genomics studies, existing technologies still have some shortcomings in meeting the requirements of balancing high compression rates, preservation of rare variants, and rapid variant lookup. Summary of the Invention

[0004] To address or at least partially address the aforementioned technical problems, this application provides a genome sequence compression method and system. This system can distinguish between sequencing errors and real variations using a temporal differential neural network model. It employs a differential compression coding strategy to improve compression efficiency while protecting important biological variation information. Furthermore, it constructs a variation detection-friendly format for rapid variation querying, providing an efficient and accurate data compression and analysis solution for large-scale population genomics research.

[0005] In a first aspect, this application provides a method for genome sequence compression, the method comprising the following steps: Acquire genome sequencing data, which includes base sequences and corresponding quality fractions; The sequencing data is compared with a reference genome to determine the differential sites and read positions; The differentially expressed sites are identified to distinguish between sequencing errors and true variants; Differential compression coding is performed based on the identification results: lossless compression coding is used for true variants, simplified coding is used for sequencing errors, and a zero-storage strategy is used for indifferent sites, thereby generating compressed data files.

[0006] Optionally, the identification of differential sites includes the following steps: Based on sequencing data and alignment results, the sequence feature matrix of each differentially expressed site is extracted; By analyzing the sequencer's operating status information, the chemical state vector corresponding to each differential site is obtained; The sequence feature matrix and chemical state vector are input into a temporal difference neural network model, which includes a temporal feature extraction module and a differential attention classification module. Among them, the temporal feature extraction module uses a long short-term memory network to model the temporal dependencies of sequences; Among them, the differential attention classification module focuses on sequence difference features and sequencing error features through a global attention network and a noisy attention network, respectively, and outputs the error probability and mutation probability corresponding to each difference site; Based on preset error thresholds and preset mutation thresholds, the error probability and mutation probability of each differential site are compared with the corresponding thresholds to determine whether the differential site is a sequencing error or a true mutation.

[0007] Optionally, the sequence feature matrix includes: Unithermal encoding of bases, first error rate, local GC content, homopolymer length, position encoding, base type indicator, and reference sequence matching flag; The first error rate is a sequencing reliability indicator calculated based on the quality score.

[0008] Optionally, the chemical state vector includes: Second error rate, fluorescence crosstalk coefficient, synchronization error, temperature normalization, reagent replenishment cycle normalization, and cycle normalization; The second error rate is the base error rate calculated based on the sequencing reaction cycle.

[0009] Optionally, the structure of the time-series feature extraction module includes: Input layer: The sequence feature matrix and the chemical state matrix are concatenated column-wise to form the input matrix; First Long Short-Term Memory (LSTM) network layer: Receives the input matrix and processes temporal information through forget gate, input gate, and output gate mechanisms; The second long short-term memory network layer receives the output of the first long short-term memory network layer and further extracts deep temporal features; Output layer: The hidden state output of the second long short-term memory network layer is used as a temporal feature vector for the differential attention classification module.

[0010] Optionally, the structure of the differential attention classification module includes: Global Attention Network: Includes query matrix generation layer, key matrix generation layer, and value matrix generation layer. It receives temporal feature vectors, calculates attention weights through the softmax function, and outputs global attention features. Noise Attention Network: Includes query matrix generation layer, key matrix generation layer, and value matrix generation layer, specifically for identifying sequencing error patterns, receiving temporal feature vectors and outputting noise attention features; Weight mapping layer: Receives global attention features and noisy attention features, and maps them into scalar weights through fully connected layers respectively; Differential score calculation layer: Calculates the differential score based on the scalar weights and time-series weight factors; Classification output layer: Receives the difference scores and converts them into the error probability and mutation probability of each differential site through a fully connected layer and a softmax function.

[0011] Optionally, the method further includes constructing a variant detection-friendly format: Compressed data that can be quickly located at genomic positions by constructing a multi-level index structure; A variant metadata table is constructed, and the compressed data is reorganized into a structured index that supports fast querying. The metadata table contains sample identifier, genome location, variant type, reference base, observed base, confidence level, and compressed data offset information. A multidimensional index is built based on the aforementioned variant metadata table, including a location index, a confidence index, and a variant type index, to quickly filter the genomic locations located by the above multi-level index structure; When performing a mutation query, the location of the compressed data block in the target region is first located through the multi-level index structure. Then, the mutation metadata table is used to pre-filter based on conditions such as confidence level and mutation type. Only the compressed data block corresponding to the filtering results is accessed and parsed. Provides a format conversion interface to convert query results into a standard variant call format and output them.

[0012] Optionally, the differential compression encoding includes: For locations identified as true variants, a complete variant record is constructed and lossless compression is performed using Huffman coding. The variant record includes chromosome number, genomic location, variant type, reference base, observed base, confidence level, and quality score. For locations identified as sequencing errors, a graded compression strategy is adopted based on the error probability. For positions determined to be indifferent, base information and mass fraction are not stored; instead, reconstruction is performed using a reference sequence.

[0013] In a second aspect, this application provides a genome sequence compression system, the system comprising a processor and a memory, the memory storing at least one instruction, at least one program, code set or instruction set, the at least one instruction, the at least one program, the code set or instruction set being loaded and executed by the processor to implement the genome sequence compression method as described in any of the first aspects.

[0014] The technical solution provided in this application has the following advantages compared with the prior art: One of its beneficial effects and its working principle is as follows: Existing lossless compression methods retain all original data, including a certain proportion of sequencing errors, resulting in low compression efficiency. While lossy compression methods improve compression ratios through simple quality threshold filtering, this filtering strategy may inadvertently delete genuine biological variants. This is because existing methods cannot accurately distinguish between sequencing errors and genuine variants, and can only employ either a conservative full-preservation strategy or an aggressive simple filtering strategy.

[0015] This application identifies sequencing errors and true variants, uses 1-2 bit simplified encoding for high-confidence sequencing errors and Huffman lossless encoding for true variants, achieving a lower compression rate while ensuring a high retention rate of true variants, thereby reducing storage costs.

[0016] Its second beneficial effect and its working principle are as follows: Existing error and variant identification methods lack accuracy in complex sequence environments, particularly exhibiting a high false positive rate for rare variants, such as homopolymers and low-complexity sequence regions. This problem arises because existing error and variant identification methods ignore the temporal flow characteristics and physicochemical properties of the sequencing process, relying solely on static sequence information and quality scores for judgment. They fail to fully utilize the differences in error patterns at different reaction cycles of the sequencer.

[0017] This application proposes a temporal differential neural network model, which for the first time introduces the temporal flow characteristics and chemical kinetics of the sequencing process into variant identification. First, at the chemical kinetics modeling level, by analyzing the physicochemical processes of the sequencer, a chemical state parameter model is established, including the second error rate, fluorescence crosstalk coefficient, and synchronization error, quantifying the error rate variation patterns and chemical reagent decay modes across different sequencing cycles. Second, at the temporal feature extraction level, an LSTM network is used to model the temporal dependencies of sequences, capturing long-window sequence context information and quality evolution patterns, identifying temporal features that traditional methods cannot detect. Third, at the differential attention mechanism level, by designing two parallel modules—global attention and noise attention—to focus on the global differential features of the sequence and the specific patterns of sequencing errors, respectively, accurate separation and comprehensive judgment of complex features are achieved. This method utilizes the physicochemical characteristics and temporal evolution laws of the sequencing process to improve the identification accuracy in rare variants and complex sequence regions.

[0018] Therefore, through time-series modeling and chemical kinetic analysis, this application is able to more accurately distinguish between biological real variations and technical sequencing errors, reducing the false positive rate and ensuring that important low-frequency variation information is fully protected.

[0019] Its third beneficial effect and its working principle are as follows: Traditional genome compression formats (such as BAM and CRAM) require complete or partial decompression before variant detection, which is time-consuming and consumes a large amount of temporary storage space, resulting in low efficiency. Furthermore, the new compression format proposed in this application differs significantly from traditional formats (such as BAM and CRAM), and existing variant detection tools cannot directly parse it, requiring decompression and reconstruction into the traditional format, further exacerbating the efficiency problem.

[0020] This application implements a collaborative query optimization process by designing a variant detection-friendly format: First, a B+ tree index structure is used for rapid location; then, an innovative variant metadata table is used to pre-screen data within the compressed domain based on multiple dimensions such as confidence level and variant type; only the compressed data blocks corresponding to the selected high-quality variants are accessed; finally, a VCF format adapter is used to convert the query results into an industry-standard format. This hierarchical collaborative query mechanism avoids the resource waste caused by accessing all data, shortens the time for variant frequency queries, and provides an efficient solution for large-scale population genomics research. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the genome sequence compression method provided in this application embodiment. Detailed Implementation

[0022] The technical solutions in this application will now be described with reference to the accompanying drawings.

[0023] Many specific details are set forth in the following description to provide a thorough understanding of this application, but this application may also be implemented in other ways different from those described herein; obviously, the embodiments in the specification are only some embodiments of this application, and not all embodiments. It should be noted that, unless otherwise specified, the embodiments of this application and the features in the embodiments can be combined with each other.

[0024] Firstly, such as Figure 1 As shown, this application provides a genome sequence compression method. The following describes the implementation process of this method in detail, using the genome compression process of processing whole-genome sequencing data from 10,000 individuals for rapid variant frequency statistics as an example.

[0025] S101: Data Acquisition and Feature Extraction.

[0026] We analyzed whole-genome data from 10,000 individuals from different regions to rapidly count variant frequencies and identify patterns of genetic diversity in the population. The original sequencing data was enormous and required compression to reduce storage requirements while still supporting rapid variant lookup. The study used the Illumina NovaSeq platform, with a sequencing depth of 30× per sample and a read length of 150 base pairs (bp).

[0027] The system reads sequencing data from the FASTQ file of the first sample, SAMPLE_001. The FASTQ file is a standard text format for storing biological sequences and their quality information, acquired through the Illumina NovaSeq platform.

[0028] Each read segment in the FASTQ file includes: a base sequence and a quality fraction string.

[0029] Examples of base sequence and mass fraction strings are shown below: ATCGATCGATCGAAAAATCGATCGATCG… (150 bases); IIIIIIIIHHHHHFFFFEEDDDCCCCBB… (150 quality characters); Each character in the quality score string corresponds to a sequencing quality of one base. The characters represent the Phred quality score using ASCII encoding, and the conversion formula is: For example, the ASCII code for the character "I" is 73, and the corresponding Phred quality score is: 73-33=40.

[0030] To identify potential variants, the system uses BWA-MEM (Burrows-Wheeler Aligner's MEM algorithm) to align reads with the human reference genome GRCh38, generating the final alignment position of the reads on the reference sequence and the corresponding reference sequence information, determining the differential sites between the reads and the reference sequence, and generating the alignment information mentioned in S103.

[0031] S1011: Obtain the sequence feature matrix.

[0032] The construction of the sequence feature matrix consists of the following components, forming a 150×10 feature matrix: One-hot encoding (4-dimensional): One-hot encoding is used for the 150bp base sequence to facilitate neural network processing. A is encoded as [1,0,0,0], T as [0,1,0,0], C as [0,0,1,0], and G as [0,0,0,1]. Each base occupies 4 dimensions, forming a 150×4 submatrix. First error rate (1-dimensional): The Phred quality score Q is converted into a first error rate P to quantify the sequencing reliability of each base: Where P is the first error rate and Q is the Phred quality score.

[0033] For example, the quality character "F" at position 13 corresponds to a Phred quality score of 37, and the calculated first error rate P = 0.0002. This forms a 150×1 probability vector. Local GC content (1-dimensional): Regions with high GC content are more prone to sequencing errors. For each position, the ratio of G to C bases within a 5bp window (11bp total) before and after it is calculated. For example, the window sequence at position 13 is “TCGAAAAATCG”, with a GC content of 3 / 11 ≈ 0.27. A 150×1 GC content vector is then generated. Homopolymer length (1-dimensional): Homopolymer regions (such as AAAAA) are prone to polymerase slippage, leading to errors. Identify the homopolymer length at each location; for example, location 13 is within the "AAAAA" pentamer, with a homopolymer length of 5. Form a 150×1 length vector. Position encoding (1-dimensional): The relative positions of the bases in the read segment are normalized to the [0,1] interval. The relative position of position 13 is 13 / 150 ≈ 0.087. This forms a 150×1 position vector. Base type indication (1-dimensional): The base types are encoded as numerical values: A=0.25, T=0.50, C=0.75, G=1.00. Position 13 is an A base, encoded as 0.25. This forms a 150×1 type vector. Reference sequence matching flag (1-dimensional): The position is marked as 0 if it matches the reference sequence, and 1 if it doesn't, directly indicating the location of potential variants or errors. Reference sequence matching markers can assist time-series difference neural network models in initially screening differentially expressed sites that require further analysis. Position 13 is an insertion site, marked as 1. This forms a 150×1 marker vector. By combining the above seven feature components, a complete 150×10 sequence feature matrix is ​​constructed, which serves as the input to the temporal difference neural network model in S102.

[0034] S1012: Obtain the chemical state matrix.

[0035] By analyzing the sequencer log files and the temporal patterns of sequencing errors, a 6-dimensional chemical state vector is generated for each base position, forming a 150×6 chemical state matrix, which provides key input for subsequent temporal feature extraction.

[0036] First, extract the timing metadata from the sequencer log file. Taking the read from sample SAMPLE_001 as an example, this read comes from flow cell batch 20231215, channel 1, coordinates (200, 300). Each base position corresponds to a sequencing reaction cycle. For example, position 13 corresponds to the 13th sequencing reaction cycle, the reaction temperature is 65.2°C, and the reagent replenishment status is the 13th cycle (reagent is replenished every 25 cycles).

[0037] Based on the time-dependent nature of sequencing errors, the system calculates the following chemical state parameters and generates a 6-dimensional vector: Second error rate (1-dimensional): Based on the decay pattern of chemical reagents over time, the error rate changes with the sequencing cycle (increasing over time): Where c is the number of sequencing reaction cycles, This is the second error rate for that period; For example, for position 13 (c=13), since c<50, the second error rate is 0.001; This ultimately forms a second error rate vector of 150×1; Fluorescence crosstalk coefficient (1-dimensional): In the Illumina four-color fluorescence system, A bases use red fluorescence and C bases use green fluorescence. Signals from adjacent channels can interfere with each other. In the Illumina four-color fluorescence system, four different fluorescent dyes are used to label four bases: A (red fluorescence, excitation / emission wavelength approximately 660 / 680 nm), C (green fluorescence, approximately 520 / 540 nm), G (orange fluorescence, approximately 580 / 600 nm), and T (yellow fluorescence, approximately 540 / 560 nm). Due to spectral overlap, signal crosstalk occurs between adjacent channels, the severity of which depends on the base combination and the sequencing cycle. Fluorescence crosstalk coefficient Calculate using the following formula: Where c represents the number of sequencing cycles, Indicates the preceding base. Indicates the current base, This represents the time-related crosstalk factor, reflecting the effect of signal attenuation during the sequencing cycle. This represents the base correlation crosstalk factor, which reflects the effect of spectral overlap between base combinations. Among them, time-related crosstalk factor The definition is as follows: Among them, base-related crosstalk factor The definition is as follows: For example, in a sequencing reaction cycle of c=13, bases =“A”, base In the case of "A", since c < 30, =0.01, and =1.0, therefore the fluorescence crosstalk coefficient for: The final result is a 150×1 fluorescence crosstalk coefficient vector; Synchronization error (1-dimensional): During sequencing, some DNA molecules may complete their reactions earlier or later, leading to asynchronous signaling. Where c represents the number of sequencing cycles, Let c be the synchronization error probability of the c-th cycle. =0.0001 is the synchronization error coefficient. For c=13: This ultimately results in a 150×1 synchronization error vector; Temperature normalized value (1-dimensional): Normalize the reaction temperature by dividing it by the maximum value of 100°C; For example, at location 13, the temperature is 65.2°C, and the normalized value is 65.2 / 100 = 0.652.

[0038] The final result is a 150×1 temperature normalized vector; Normalized value of reagent replenishment cycle (1-dimensional): The normalized value of the reagent replenishment cycle (in this example, replenishment is once every 25 cycles).

[0039] For example, position 13 (13th cycle), reagent age = 13 / 25 = 0.52; The final result is a 150×1 reagent replenishment cycle normalized vector; Periodic normalized value (1-dimensional): Normalize by dividing the sequencing cycle by the read length of 150; For example, for position 13 (cycle=13), the normalized value is approximately 13 / 150 ≈ 0.087; This ultimately results in a 150×1 periodic normalized vector; Based on the above six chemical features, a six-dimensional chemical state vector is generated for each base position, forming a 150×6 chemical state matrix, which serves as the input to the temporal difference neural network model in S102.

[0040] S102: Temporal difference neural network model processing.

[0041] The 150×10 sequence feature matrix and 150×6 chemical state matrix obtained from S101 are input into the temporal differential neural network model. This model accurately determines whether each differential position is a sequencing error or a real variant through two core modules: temporal feature extraction and differential attention classification.

[0042] The temporal difference neural network model is an end-to-end supervised learning model, in which LSTM temporal feature extraction and differential attention classification are the two core modules. During training, only the model's final output (error probability P) is evaluated. 错误 and the probability of mutation P 变异 The intermediate results do not need to be labeled. The model learns the complete feature mapping relationship between sequencing errors and real variants through end-to-end training.

[0043] S1021: Extract temporal features.

[0044] Input the sequence feature matrix and chemical state matrix of a read segment into the time series feature extraction module to output a time series feature vector.

[0045] The temporal feature extraction module uses a two-layer long short-term memory (LSTM) network to model the temporal dependencies of sequences.

[0046] Input data preparation: The 150×10 sequence feature matrix and the 150×6 chemical state matrix are concatenated column-wise to form a 150×16 input matrix. Each row represents the complete feature of a base position, including sequence features and chemical state information.

[0047] LSTM network configuration: A two-layer unidirectional LSTM structure is used, with an input dimension of 16, a hidden layer dimension of 256, and a dropout rate of 0.2. LSTM captures long-range dependencies in sequences through a gating mechanism, the core of which includes a forget gate (determining how much historical information to retain), an input gate (determining how much new information to receive), and an output gate (determining how much information to output). The specific structure includes: Input layer: The 150×10 sequence feature matrix and the 150×6 chemical state matrix are concatenated column by column to form a 150×16 input matrix; The first LSTM layer has an input dimension of 16 and a hidden layer dimension of 256. It processes temporal information through forget gate, input gate, and output gate, and outputs the first hidden state. The second LSTM layer receives the hidden state of the first layer and further extracts deep temporal features, with an output dimension of 256. Output layer: The final hidden state of the second LSTM layer is used as the temporal feature vector h. 13 Output.

[0048] Temporal difference neural network models can use LSTM networks to identify the temporal features of bases: For example, for position 13 in the example, the LSTM network identifies the following key temporal features by analyzing its context sequence and chemical state: Sequence characteristics: Positions 10-14 form AAAAA pentamers, which is a typical characteristic of polymerase slippage; Quality characteristics: The quality score in this region remained stable (Phred values ​​were all >35), indicating good sequencing quality; Chemical characteristics: The early cycle (cycle=13) is characterized by a low error rate (0.001) and low crosstalk (0.01); The LSTM outputs a 256-dimensional temporal feature vector h. 13 This vector comprehensively encodes the sequence context, quality distribution, and chemical evolution information at position 13. Simultaneously, the hidden state of the LSTM propagates along the sequence direction, allowing each position to utilize the accumulated information from previous positions.

[0049] S1022: Differential attention classification.

[0050] The temporal feature vector obtained in S1021 is input into the differential attention classification module, which outputs the error probability and mutation probability to achieve accurate differentiation between sequencing errors and real mutations, and generates detailed classification description information.

[0051] The core idea of ​​the differential attention mechanism is to use two parallel multi-head attention network modules to focus on different feature patterns: a global attention network module that focuses on all sequence difference features and a noise attention network module that specifically identifies sequencing error features. The global attention network module is calculated as follows: in, , , These represent the query matrix, key matrix, and value matrix of the global attention network module, respectively. =32 indicates the dimension of the key vector, used to scale the dot product to stabilize the gradient. The softmax function normalizes the attention scores to a probability distribution, ensuring that the output is a weighted vector of values.

[0052] The noise attention network module is calculated as follows: in, , , These represent the query matrix, key matrix, and value matrix of the noisy attention network module, respectively. =32 indicates the dimension of the key vector, used to scale the dot product to stabilize the gradient. The softmax function normalizes the attention scores to a probability distribution, ensuring that the output is a weighted vector of values.

[0053] Output vector of the global attention network module and the output vector of the noise attention network module Mapped to scalar weights via fully connected layers, the scalar weights include global attention weights. and noise attention weights ; Among them, global attention weight This reflects the salience of global features; Noise attention weight This reflects the contribution of noise characteristics; For example, at a specific sequencing location, the model outputs... =0.85, indicating that the global features are significant; =0.30, indicating that the noise characteristics of the early cycle are relatively weak.

[0054] Based on global attention weights and noisy attention weights, combined with temporal weight factors, a difference score (DiffScore) is calculated to synthesize the contributions of both: Where λ = 0.5 is a learnable equilibrium parameter. =1.0 is the time-series weighting factor calculated based on chemical state. Substitute the value: Based on the difference score, the final classification probability is calculated through a fully connected layer and a softmax function.

[0055] A high difference value (0.70) indicates that the position is more likely to be a true variant. For example, for the base at position 13 in the example, the temporal difference neural network model outputs: P 错误 =0.25, P 变异 =0.85. Where P 错误 and P 变异 The sequencing error probability and the true mutation probability are independently represented, and they are not additive. Using a temporal difference neural network model, the system outputs a classification probability for each difference position and determines that the position is a true mutation based on a preset mutation threshold of 0.85.

[0056] S1023: Post-processing analysis and descriptive information generation.

[0057] Based on the classification output of the temporal difference neural network model, the system generates detailed descriptive information through an independent post-processing module, providing complete data input for S103 differential compression coding.

[0058] By analyzing the read and reference sequence data, basic positional information can be extracted, including position number, genomic coordinates, observed base, reference base, and quality score. For example, for base at position 13, its basic positional information is: position number 13, genomic coordinates chr1:1234567 (chr1 indicates that the base is located on the first chromosome, and 1234567 indicates that the read corresponds to the 1234567th base in the reference sequence; genomic coordinates are used to locate the read's position in the reference sequence so that subsequent tasks can reconstruct the genome sequence based on the reference sequence), observed base GA, reference base G (from the reference sequence), and quality score 37.

[0059] Based on the classification probability output of the temporal difference neural network model, the error probability and mutation probability at a specific position can be obtained. For example, for the base at position 13, P... 错误 =0.25, P 变异 =0.85.

[0060] According to the error probability and mutation probability, the classification result of the specific position is determined by combining the preset mutation threshold and preset error threshold: If P 错误 is greater than or equal to the preset error threshold of 0.85, it is determined as an error; If P 变异 is greater than or equal to the preset mutation threshold of 0.85, it is determined as a mutation; If P 错误 is less than the preset error threshold of 0.85 and P 变异 is less than the preset mutation threshold of 0.85, it is determined as an indifferent base.

[0061] According to the position classification result, the system performs judgment on the specific error type or mutation type for positions determined as errors or mutations. The judgment of the specific error type or mutation type is performed by automatically comparing the observed base sequence with the reference sequence through a sequence alignment algorithm, and the specific steps are: Extract the sequence window of the specific position and 5bp before and after it, then use a local alignment algorithm (such as Smith-Waterman) to identify difference patterns, and finally perform automatic classification according to the alignment result; For sequencing errors, there are four types set in total: single base substitution, insertion, deletion, and unknown; For real mutations, there are four types set in total: single base substitution, insertion, deletion, and others.

[0062] In summary, the final data includes chromosome number, genomic position, mutation type, error type, reference base, observed base, confidence level, and quality score. This ensures that the S103 differential compression coding step has complete, accurate and traceable data input.

[0063] S103: Differential compression coding.

[0064] According to the classification result, the system executes corresponding compression strategies for different types of sites: For P 变异 positions where P is greater than or equal to the preset mutation threshold of 0.85 are determined as real mutations, a complete mutation record is constructed for them and lossless compression is performed. The mutation record includes: chromosome number (chr1), genomic position, mutation type, reference base, observed base (GA), confidence level, and quality score. All data contained in the mutation record can be obtained from step S102.

[0065] For example, for position 13 in this example, the P output by the temporal differential neural network model 变异 is 0.85, this position is determined as a real mutation, and the mutation record is constructed as follows: <chr1, 1234579, insertion, G, GA, 0.85, 37> The genomic position "1234579" is obtained by adding the specific base position 13-1 in the read segment to the start position 1234567 of the read segment.

[0066] Considering the significant frequency differences in the content of mutation records, such as the fact that single base substitutions typically account for 85% of all mutation types, Huffman coding is used to losslessly compress the mutation locations. The core idea of ​​Huffman coding is to encode nodes based on their frequency of occurrence. The node with the highest frequency of occurrence is encoded with 1 bit (0), the node with the second highest frequency of occurrence is encoded with 2 bits (10), and so on, to achieve the compression function.

[0067] The specific encoding process includes: Chromosome numbers are encoded using a fixed 5-bit code, such as chr1 being encoded as 00001; Incremental coding is used to encode genomic locations, that is, coding is done with an increment relative to the previous variant location, in order to save storage space; The variant type, reference base, observed base, confidence level, and quality score are encoded. The variant type and reference base are encoded using 2 bits (corresponding to 4 variant types and 4 reference bases), the observed base is encoded using 4 bits (there may be "insertion" variant types, the corresponding observed base is "GA" which requires 4 bits of encoding), the confidence level is encoded using 7 bits, and the quality score is encoded using 4 bits.

[0068] The complete encoding of the mutation record is 31 bits in size.

[0069] For P 错误 Positions with an error value greater than or equal to a preset error threshold of 0.85 are identified as sequencing errors, and graded compression is performed based on the error probability. The specific process includes: For P 错误 Positions with an error value greater than or equal to a preset high error threshold of 0.9 are identified as high-confidence errors, and the most simplified encoding is applied to these positions. This most simplified encoding involves directly marking the erroneous base at that position as "N" (unknown base), without preserving the original base and quality information; only a 1-bit flag is required. For P 错误 Positions with an error value greater than or equal to a preset error threshold of 0.85 and less than a preset high error threshold of 0.9 are identified as medium-confidence errors, and the minimum necessary information for these positions is retained for possible subsequent analysis. This retained minimum necessary information includes the error type and the observed base at that position, both encoded using 2 bits, requiring a total of 4 bits for encoding.

[0070] For P 错误 Less than 0.85 and P变异 Positions with a mass fraction less than 0.85 are considered indistinguishable positions, completely identical to the reference sequence, and a zero-storage strategy is applied to these positions. This zero-storage strategy includes not storing the base information and mass fraction of indistinguishable positions. Based on their indistinguishability from the reference sequence, they can be quickly reconstructed from the reference sequence without requiring any storage.

[0071] For each read, alignment information also needs to be stored to locate the start position of the read in the reference sequence. This alignment information is obtained in step S101. The core function of the alignment information is to support genome sequence reconstruction. Because this method employs a zero-storage strategy for indifferent positions, if subsequent research tasks require the base information at these positions, reconstruction using the reference sequence is necessary. The specific reconstruction process includes: Based on the alignment information, the starting position of the read in the reference sequence is determined. Then, combined with the stored variation and error information, the bases in the corresponding region of the reference sequence are replaced or corrected to the actual sequencing results, thereby completely restoring the original genome sequence.

[0072] S104: Construct a variant detection-friendly format.

[0073] Existing fast mutation lookup algorithms are designed for traditional BAM / CRAM and VCF data compression formats and cannot directly handle the novel compressed data format generated in steps S101-S103 of this embodiment. To support fast mutation lookup, a specialized mutation detection-friendly format and corresponding index structure need to be designed.

[0074] To support rapid mutation queries, the system constructs a multi-level index structure and a standardized metadata format based on the differentiated compressed data of S103.

[0075] B+ tree index construction process: Root node: forms 24 second-level branches, corresponding to the 24 chromosomes of the human genome (chromosomes 1-22, X and Y sex chromosomes). The second layer: Each chromosome is divided into 50 megabase (Mb) intervals, forming the third layer of branches; The third layer: Each 50Mb interval is further subdivided into 5 10Mb sub-intervals, forming leaf nodes; Leaf nodes: Each leaf node manages 10 1000bp data blocks and records the offset address and length of all mutated compressed data within that interval.

[0076] Each node in the index contains: a start position, an end position, and pointers to child nodes (internal nodes) or data block addresses (leaf nodes).

[0077] For example, for the insertion mutation example at position 13, the index location path is: root node, chr1 branch, 1-50Mb interval, 1-10Mb sub-interval, containing the data block at positions 1234567, and the specific compressed data block is quickly located through the 4-level index.

[0078] Construct the mutation metadata table: Traditional indexing methods (such as B+ trees) can provide location lookup functionality, but they cannot use the confidence information output by the S102 neural network for filtering, nor can they directly process the variant classification results after S103 differential compression.

[0079] The core function of the mutated metadata table is to reorganize the differentially compressed data from step S103 into a structured index that supports fast queries, avoiding the need to parse all compressed data blocks for each query. The steps for building metadata records include: The compressed data from step S103 is parsed into standardized metadata records. These metadata records include sample ID, chromosome, location, variant type, reference base, observed base, confidence level (variation probability), quality score, and compressed data offset (the position of the variant within the compressed file). A multidimensional index structure is then established based on these metadata records. The position index maps chr1:1234567 to record 1, chr1:1234653 to record 2, and so on; The confidence index maps conditions with a confidence level greater than or equal to 0.95 to record 1, conditions with a confidence level greater than or equal to 0.9 and less than 0.95 to record 2, and so on. The variant type index maps the insert type to record 1, the replacement type to record 2, and so on.

[0080] By constructing a variant metadata table, the query process can be simplified to three steps: first, records that meet the conditions are directly filtered in the metadata table; second, only the corresponding compressed data blocks are accessed for the filtered results to reduce disk I / O operations; and finally, the query results are returned directly, avoiding a large amount of decompression calculation. Simultaneously, the system achieves direct and efficient querying of the new S103 compression format, fully utilizing the information output by the S102 temporal difference neural network model for filtering—a core advantage that existing indexing algorithms based on the BAM / CRAM format cannot achieve.

[0081] VCF format adapter implementation: The primary function of the VCF format adapter is to ensure seamless integration of the new compression format generated by S103 with the existing ecosystem of bioinformatics tools. Since the vast majority of genome analysis tools (such as GATK, SAMtools, and Picard) are developed based on the standard VCF format, directly using the new format would lead to compatibility issues.

[0082] The conversion process includes: Extract variant records from compressed data; Field mapping and conversion are performed according to VCF standard format requirements; Generate standard VCF record output.

[0083] For example, the insertion variation at position 13 is transformed into: chr1 1234579 . G GA 85 PASS AF=0.0001 DP=30 In this table, chr1 represents the first chromosome, 1234579 (1234567+13-1) represents the position, "." indicates no dbSNP number, G is the reference base, GA is the post-mutation sequence (observed bases from the mutation record), 85 is the quality score (based on a confidence level of 0.85), PASS indicates passing the filter (based on a confidence level greater than or equal to 0.85), AF=0.0001 represents the allele frequency in the population (calculated based on statistics from 10,000 samples), and DP=30 represents the sequencing depth at this position (from sequencing parameter settings).

[0084] With the VCF adapter, users can directly use existing genome analysis tools for downstream analysis, including variant annotation, population genetics analysis, and variant effect prediction, achieving full compatibility between the new compression format and existing toolchains.

[0085] S105: Perform variation frequency analysis.

[0086] After performing differential compression on 10,000 samples, the researchers conducted a variation frequency analysis.

[0087] First, a B+ tree index constructed using S104 is used for rapid location. The multi-level index structure of the B+ tree supports efficient location from the chromosome level to specific genomic locations, quickly finding the location of compressed data blocks in the target region by following the path from the root node to the leaf node.

[0088] Secondly, pre-screening is performed using the mutation metadata table in S104. The metadata table integrates the confidence output of the S102 temporal difference neural network model and the mutation classification results of S103 differential compression, supporting rapid screening based on multiple dimensions such as confidence threshold and mutation type, thus avoiding the need to parse each compressed data block one by one.

[0089] The system then reads the differentially compressed data blocks from S103 as needed. Based on the pre-screening results of the metadata table, only compressed data blocks containing the target mutations are accessed to reduce disk I / O operations. The system parses the mutation records in the compressed data blocks and uses S103's Huffman coding to achieve efficient decompression.

[0090] Next, the system performs statistical calculations on the mutation frequency. Based on the confidence information output by the S102 neural network, quality filtering is performed, counting only high-quality mutations that meet preset threshold requirements to ensure the reliability of the frequency analysis results.

[0091] Finally, the system outputs standardized results through the S104 VCF format adapter. This converts the internally compressed query results into the standard VCF format, ensuring compatibility with existing genome analysis toolchains and supporting subsequent downstream applications such as variant annotation and population genetics analysis.

[0092] Through the above process, this method achieves efficient variant frequency analysis. Compared to traditional methods that require complete decompression of the BAM file and re-detection of variants, this method utilizes pre-processed compressed data and an index structure to improve query performance while maintaining the accuracy of the analysis results. The entire system provides a complete solution for large-scale population genomics research, from data compression to rapid analysis.

[0093] Through a complete implementation process, the system successfully achieved efficient compression of large-scale genome sequencing data. The compressed data consists of real variant data (using Huffman coding), sequencing error locations (using 1-2 bit simplified coding), indexes, and metadata. Furthermore, the system improved the retention rate of real variants, especially rare variants, and the accuracy of sequencing error identification to a certain extent.

[0094] Secondly, embodiments of this application provide a genome sequence compression system, including a processor and a memory, wherein the memory stores at least one instruction, at least one program, code set, or instruction set, and the at least one instruction, the at least one program, the code set, or instruction set is loaded and executed by the processor to implement the genome sequence compression method as described in the above embodiments.

[0095] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element. Moreover, in the description of the embodiments of this application, unless otherwise stated, " / " signifies "or," for example, A / B can mean A or B; "and / or" in this document is merely a description of the relationship between related objects, indicating that three relationships can exist, for example, A and / or B can represent: A alone, A and B simultaneously, and B alone. Furthermore, in the description of the embodiments of this application, "multiple" refers to two or more.

[0096] The above description is merely a specific embodiment of this application, enabling those skilled in the art to understand or implement this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments described herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for compressing genome sequences, characterized in that, The method includes the following steps: Acquire genome sequencing data, which includes base sequences and corresponding quality fractions; The sequencing data is compared with a reference genome to determine the differential sites and read positions; The differentially expressed sites are identified to distinguish between sequencing errors and true variants; Differential compression coding is performed based on the identification results: lossless compression coding is used for true variants, simplified coding is used for sequencing errors, and zero-storage strategy is used for indifferent sites, thereby generating compressed data files; The identification of differential sites includes the following steps: Based on sequencing data and alignment results, the sequence feature matrix of each differentially expressed site is extracted; By analyzing the sequencer's operating status information, the chemical state vector corresponding to each differential site is obtained; The sequence feature matrix and chemical state vector are input into a temporal difference neural network model, which includes a temporal feature extraction module and a difference attention classification module. The chemical state vector includes: Second error rate, fluorescence crosstalk coefficient, synchronization error, temperature normalization, reagent replenishment cycle normalization, and cycle normalization; The second error rate is the base error rate calculated based on the sequencing reaction cycle. Among them, the temporal feature extraction module uses a long short-term memory network to model the temporal dependencies of sequences; Among them, the differential attention classification module focuses on sequence difference features and sequencing error features through a global attention network and a noisy attention network, respectively, and outputs the error probability and mutation probability corresponding to each difference site; Based on preset error thresholds and preset mutation thresholds, the error probability and mutation probability of each differential site are compared with the corresponding thresholds to determine whether the differential site is a sequencing error or a true mutation.

2. The genome sequence compression method according to claim 1, characterized in that, The sequence feature matrix includes: Unithermal encoding of bases, first error rate, local GC content, homopolymer length, position encoding, base type indicator, and reference sequence matching flag; The first error rate is a sequencing reliability index calculated based on the quality score.

3. The genome sequence compression method according to claim 1, characterized in that, The structure of the temporal feature extraction module includes: Input layer: The sequence feature matrix and the chemical state matrix are concatenated column-wise to form the input matrix; First Long Short-Term Memory (LSTM) network layer: Receives the input matrix and processes temporal information through forget gate, input gate, and output gate mechanisms; The second long short-term memory network layer receives the output of the first long short-term memory network layer and further extracts deep temporal features; Output layer: The hidden state output of the second long short-term memory network layer is used as a temporal feature vector for the differential attention classification module.

4. The genome sequence compression method according to claim 3, characterized in that, The structure of the differential attention classification module includes: Global Attention Network: Includes query matrix generation layer, key matrix generation layer, and value matrix generation layer. It receives temporal feature vectors, calculates attention weights through the softmax function, and outputs global attention features. Noise Attention Network: Includes query matrix generation layer, key matrix generation layer, and value matrix generation layer, specifically for identifying sequencing error patterns, receiving temporal feature vectors and outputting noise attention features; Weight mapping layer: Receives global attention features and noisy attention features, and maps them into scalar weights through fully connected layers respectively; Differential score calculation layer: Calculates the differential score based on the scalar weights and time-series weight factors; Classification output layer: Receives the difference scores and converts them into the error probability and mutation probability of each differential site through a fully connected layer and a softmax function.

5. The genome sequence compression method according to claim 1, characterized in that, The method also includes constructing a variant detection-friendly format: Compressed data that can be quickly located at genomic positions by constructing a multi-level index structure; A variant metadata table is constructed, which reorganizes the compressed data into a structured index that supports fast querying. The metadata table contains sample identifier, genome location, variant type, reference base, observed base, confidence level, and compressed data offset information. A multidimensional index is built based on the aforementioned variant metadata table, including a location index, a confidence index, and a variant type index, to quickly filter the genomic locations located by the above multi-level index structure; When performing a mutation query, the location of the compressed data block in the target region is first located through the multi-level index structure. Then, the mutation metadata table is used to pre-filter based on confidence level and mutation type conditions. Only the compressed data blocks corresponding to the filtered results are accessed and parsed. Provides a format conversion interface to convert query results into a standard variant call format and output them.

6. The genome sequence compression method according to claim 1, characterized in that, The differential compression encoding includes: For locations identified as true variants, a complete variant record is constructed and lossless compression is performed using Huffman coding. The variant record includes chromosome number, genomic location, variant type, reference base, observed base, confidence level, and quality score. For locations identified as sequencing errors, a graded compression strategy is adopted based on the error probability. For positions determined to be indifferent, base information and mass fraction are not stored; instead, reconstruction is performed using a reference sequence.

7. A genome sequence compression system, characterized in that, The system includes a processor and a memory, the memory storing at least one instruction, at least one program, a code set, or an instruction set, the at least one instruction, the at least one program, the code set, or the instruction set being loaded and executed by the processor to implement the genome sequence compression method as described in any one of claims 1-6.

Images