A long-read umi sequence correction and consensus sequence generation method

Abstract

The embodiment of the application discloses a long-read UMI sequence correction and consensus sequence generation method, comprising: obtaining at least one batch of to-be-processed reads; for each batch: from each read in the current batch, a first UMI to be corrected is respectively intercepted; each first UMI is deduplicated to obtain each second UMI after deduplication; the similarity of each two second UMIs is calculated according to the base difference of each two second UMIs and the average quality of each base position; the similarity threshold value used for distinguishing whether the current batch is similar is determined, and the similarity threshold value is corrected according to the global average quality of the second UMI of the current batch; a graph is constructed with each second UMI as a node and a node relationship with a similarity higher than the corrected similarity threshold value as an edge, and pruning is performed according to the abundance difference of the second UMIs at both ends of each edge, each closed subgraph in the pruned graph is divided into each UMI group, and a consensus sequence and a consensus quality of each UMI group.

Description

Technical Field

[0001] The embodiments of the present invention relate to the interdisciplinary field of bioinformatics and genomics, and in particular to a method for long read UMI sequence correction and consensus sequence generation. Background Technology

[0002] In recent years, next-generation sequencing technologies, represented by single-cell sequencing and long-read sequencing (such as PacBio SMRT sequencing and Oxford Nanopore sequencing), have continued to develop. These technologies can resolve complete transcript isoforms at the single-cell level with unprecedented resolution, greatly promoting the development of fields such as oncology, immunology, and developmental biology. To accurately quantify individual RNA molecules, UMI (Unique Molecular Identifier) ​​technology has been widely used. However, long reads also present a challenge in single-cell sequencing: the error rate of raw reads is higher than that of traditional short-read sequencing.

[0003] These sequencing errors can occur not only in the gene sequence itself but also in UMI sequences, causing the same original RNA molecule to be incorrectly labeled with multiple similar UMIs. This leads to misclassification as multiple independent molecules in subsequent analyses, severely overestimating the complexity of gene expression. To address this issue, various UMI correction algorithms have been proposed in the prior art, with the directional algorithm used by UMI-tools serving as the industry benchmark. This algorithm constructs a UMI network graph and, based on the core assumption that "sequencing errors always point from high-abundance UMIs to low-abundance UMIs," merges similar low-abundance UMIs into high-abundance UMIs, thereby achieving correction.

[0004] However, existing technologies, such as UMI-tools, suffer from the following technical flaws: their core alignment models are based on fixed-length Hamming distance or fixed-length edit distance, which rests on the fundamental assumption that UMI length remains constant. On traditional short-read sequencing platforms (such as Illumina), sequencing errors are primarily base substitutions, and UMI length typically remains unchanged, making this assumption largely valid. However, on long-read sequencing platforms (such as PacBio SMRT sequencing and Oxford Nanopore sequencing), insertion and deletion errors (collectively known as Indels) become the dominant error types, causing changes in the actual length of the UMI sequence. When UMIs become longer or shorter due to Indel errors, the fixed-length distance model completely fails: UMI sequences originating from the same original molecule and supposedly identified as identical cannot be correctly aligned and merged due to length differences, leading to severe false negatives—true homologous UMIs are incorrectly identified as different molecules and thus segmented into different UMI groups. The "missed merging" phenomenon directly leads to a serious overestimation of molecular counts, which greatly reduces the accuracy of gene expression quantification and causes an intolerable systemic bias in long-read single-cell sequencing applications.

[0005] This deficiency reveals the need for a shift in UMI correction technology: from fixed-length sequence alignment to variable-length sequence alignment. In order to truly solve the UMI correction problem in the era of long-read sequencing, a completely new technical solution is urgently needed. It must be able to: (1) process variable-length UMI sequences and achieve indel-aware intelligent alignment; (2) make full use of sequencing data information and dynamically adapt to the data quality of different platforms and batches; (3) improve the fidelity of single-molecule sequences while completing accurate molecule counting, and lay a high-quality data foundation for downstream analysis. Summary of the Invention

[0006] This invention provides a method for long read length UMI sequence correction and consensus sequence generation to solve the above-mentioned technical problems.

[0007] In a first aspect, embodiments of the present invention provide a method for long-read UMI sequence correction and consensus sequence generation, including:

[0008] Obtain at least one batch of read segments to be processed;

[0009] Perform the following operations for each batch:

[0010] S210. Extract the first UMI to be corrected from each read segment in the current batch; remove duplicates from each first UMI to obtain each second UMI after deduplication; and calculate the similarity between each pair of second UMIs based on the base difference between each pair of second UMIs and the average mass of each base position.

[0011] S220. Based on the distribution of similarity, determine the similarity threshold used to distinguish whether the current batch is similar, and correct the similarity threshold based on the global average quality of the second UMI of the current batch.

[0012] S230. Using each second UMI as a node, construct a graph with the node relationship with similarity higher than the corrected similarity threshold as an edge, and prune according to the abundance difference of the second UMI at both ends of each edge, and divide each closed subgraph in the pruned graph into each UMI group.

[0013] S240. For each UMI group: fill the reads containing each second UMI within the group with a uniform length and align them into an alignment matrix with each read as a row and each base position in the read as a column; calculate the quality weight, position weight, homopolymer weight, and coverage weight of each row in the same column of the alignment matrix, and determine the consensus bases and consensus quality of each column according to each weight, and arrange the consensus bases of each column into the consensus gene sequence of the current UMI group.

[0014] In a second aspect, embodiments of the present invention provide an electronic device, the electronic device comprising:

[0015] One or more processors;

[0016] Memory, used to store one or more programs;

[0017] When the one or more programs are executed by the one or more processors, the one or more processors implement the long read length UMI sequence correction and consensus sequence generation method described in any embodiment.

[0018] Thirdly, embodiments of the present invention also provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the long read length UMI sequence correction and consensus sequence generation method described in any embodiment.

[0019] In summary, this invention provides a method for long-read UMI sequence correction and consensus sequence generation, aiming to overcome the shortcomings of existing technologies, such as inability to handle variable-length UMI sequences, insufficient information utilization, lack of adaptability, and failure to improve sequence quality when processing high-error-rate sequencing data. The method of this embodiment can handle variable-length UMI sequences and, by introducing base quality-weighted distance, adaptive threshold calculation, and downstream consensus sequence generation, achieves ultra-high precision correction of UMI sequences and deep error correction of single-molecule sequences. Specifically, the method of this embodiment can achieve the following beneficial effects:

[0020] 1. This method offers variable-length sequence alignment capability and higher correction accuracy. It overcomes the limitations of existing technologies based on fixed-length distance alignment by implementing indel-aware variable-length sequence alignment, thus solving the UMI correction problem in the era of long-read sequencing. Instead of using simple edit distance, this method introduces a base quality (Q-score) to weight the calculated UMI distance. This approach simultaneously considers base substitutions, insertions, or deletions, assigning a weight based on its Q-score to each differential base. A mismatch of a high-quality base will incur a much higher "distance penalty" than a low-quality base mismatch, enabling the algorithm to accurately distinguish between genuine UMI differences and high-probability sequencing errors.

[0021] 2. Adaptive calculation of clustering thresholds. Addressing the quality differences across datasets, this embodiment abandons a fixed merging threshold. Before clustering the UMIs of each cell (or gene), this method first rapidly analyzes the abundance distribution, sequence similarity distribution, and average base quality of the UMIs within that subset. Based on these characteristics, a "weighted distance" threshold best suited to the current data quality is dynamically calculated, thereby achieving "personalized" processing for each data subset and maximizing the accuracy of the correction.

[0022] 3. Achieving high-fidelity consensus sequence generation. For each UMI cluster formed after correction and merging, the method in this embodiment further utilizes the complete sequence of all reads within the cluster and the corresponding base quality information. By calling multiple sequence alignment to arrange these reads, and based on the depth, consistency, and average quality of each base site, a "high-fidelity consensus sequence" with a theoretically much lower error rate than any single original read is calculated and generated. Attached Figure Description

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

[0024] Figure 1 This is a flowchart of a long-read UMI sequence correction and consensus sequence generation method provided in an embodiment of the present invention;

[0025] Figure 2 This is a flowchart of another long-read UMI sequence correction and consensus sequence generation method provided by an embodiment of the present invention;

[0026] Figure 3This is a UMI calibration accuracy-simulation cross-platform comparison chart provided by an embodiment of the present invention;

[0027] Figure 4 This is a consensus sequence fidelity-simulated cross-platform data comparison chart provided by an embodiment of the present invention;

[0028] Figure 5 This is a comparison chart of the InDel error rate of a homopolymer region provided in an embodiment of the present invention;

[0029] Figure 6 This is a comparison chart showing an enhanced detection capability for low-abundance molecules provided by an embodiment of the present invention;

[0030] Figure 7 This is a comparison chart showing how an embodiment of the present invention reduces the detection threshold for rare variants;

[0031] Figure 8 This is a comparison chart showing the trend of quantitative deviation with Indel error rate provided by an embodiment of the present invention;

[0032] Figure 9 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0033] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0034] In the description of this invention, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first," "second," and "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0035] In the description of this invention, it should also be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0036] Figure 1 This is a flowchart illustrating a method for long-read UMI sequence correction and consensus sequence generation according to an embodiment of the present invention. The method is executed by an electronic device. Figure 1 As shown, the method specifically includes:

[0037] S10. Obtain at least one batch of read segments to be processed.

[0038] This step first obtains at least one batch of reads to be processed, which serves as the data source for the entire method. Each batch includes multiple raw reads, each containing a UMI and the gene sequence from the original molecule.

[0039] Since the data quality varies from batch to batch, the specific parameters involved in subsequent operations on each batch will also be different, thus demonstrating the adaptive capability of the method in this embodiment to different batches of data.

[0040] S20. Perform operations S210-S240 for each batch to cluster reads from the same original molecule (i.e., those that should carry the same UMI) into one class, and generate a consensus sequence and consensus quality for each class of reads. Specifically:

[0041] S210. Extract the UMI to be corrected from each read segment in the current batch; remove duplicates from each UMI to obtain distinct UMIs after deduplication; and calculate the similarity between each pair of deduplicated UMIs based on the base difference between each pair of deduplicated UMIs and the average mass of each base position.

[0042] In this embodiment, it is assumed that the UMI length and start position of each read segment are the same and known. This step first reads the default UMI from each read segment based on the known UMI length and start position. For example, if the default UMI start position is the first position of the read segment and the default UMI length is 10 bits, then the first 10 bits of the sequence of each read segment are read as the UMI of each read segment. Since the UMI itself in the long read sequence may have base substitutions, additions, or deletions, the UMI sequence extracted here may not be the real UMI, but a UMI to be corrected. Subsequently, a series of operations will be performed to correct the real UMIs that should be the same to be consistent (i.e., clustered into one class).

[0043] Meanwhile, the extracted UMIs may be duplicated, meaning that the same UMI may be extracted from multiple read segments. Therefore, this embodiment first deduplicates the extracted UMIs to transform the set of extracted UMIs into a set of UMIs without duplicate elements. For ease of distinction and description, this embodiment refers to each UMI that may be duplicated directly extracted from the read segment as the first UMI, and each UMI that is distinct after deduplication as the second UMI.

[0044] Furthermore, this embodiment can also generate the following three types of metadata for each second UMI:

[0045] 1. UMI abundance Defined as "the first UMI in the current batch is exactly equal to the kth second UMI". The total number of reads for the UMI candidate molecular tag is the number of times it has been observed.

[0046] 2. Average mass at each base position in UMI Defined as sharing the same second UMI All reads in the first The average sequencing quality score (i.e., Phred quality score, or simply quality) at each base position is calculated using the following formula:

[0047] (1)

[0048] in, Indicates the index of the second UMI. , Indicates the number of the second UMI; This is a base position index within UMI. This is the default fixed length for UMI; This indicates that traversing the first UMI equals All reads, Indicates that the first UMI equals The first of all reads The first of the reads The mass of each base position.

[0049] 3. Reads Reference List: Retain all reads references for each second UMI, so that the complete sequence and quality information can be called in the subsequent consensus sequence generation stage.

[0050] Through the deduplication and recombination process described above, the current batch of reads dataset can be reorganized into a lookup table with the second UMI as the node. Each node carries three types of information: abundance, average quality at each base position, and read citations. This recombination operation does not discard any original information about the reads; it only changes the way the data is organized.

[0051] Furthermore, after deduplication of UMIs, this step calculates the similarity between every two deduplicated UMIs, providing a basis for subsequent determination of whether UMIs should be merged. In one specific implementation, the following operations can be performed for every two second UMIs:

[0052] Step 1: Determine the similarity between the two current second UMI sequences based on the number of positions with identical bases. The more positions with identical bases, the greater the similarity. Similarity is used to compare whether two second UMIs are "similar" enough to be merged; this step calculates a basic similarity. Optionally, the two current second UMIs can be denoted as... Then their basic similarity can be calculated using the following formula. :

[0053] (2)

[0054] in, express and The Hamming distance between them, that is and The number of different bases between them. (Through...) calculate and The number of identical bases between them, then divided by the uniform UMI length. Normalize the results to If the two second UMIs are exactly the same, ,but If the two second UMIs are completely different, ,but .

[0055] It should be noted that the bases in this embodiment refer to... Any of the four base types. A base at a given position can be... In any of the following, the mass of the base at that position is the mass of that base position. Regarding the number of different bases and the number of identical bases between two second UMIs, it can be understood as follows: Iterate through each base position of the two second UMIs. If the bases at the same position are identical (e.g., both are A, or both are C), then increment the number of identical bases by 1; if the bases at the same position are different (e.g., one is A, and the other is C), then increment the number of different bases by 1.

[0056] Step 2: Based on the quality of the positions where the two second UMIs have different bases, the similarity is corrected to obtain the final similarity between the two second UMIs. Optionally, this step first calculates a quality penalty value: This embodiment quantifies the confidence of "dissimilarity" by introducing base quality. If the two second UMIs have different bases at a specific position, a penalty is imposed based on the sequencing quality of these different bases.

[0057] Optionally, the formula for calculating the quality penalty value is:

[0058] (3)

[0059] in, express and Quality penalty value between, and They represent The Middle Bases at each base position; and They represent The Middle Average mass of each base position; Indicates an indicator function, when The Middle When the bases at each base position are different ,when The Middle When the bases at each base position are the same ; This represents the scaling factor, used to control the overall intensity of the quality penalty. (By...) Iterate through all base positions in UMI, then use the indicator function I. The system filters for differential sites, calculating a penalty only when two second UMIs differ in the bases at that position; otherwise, no penalty is applied. When base differences occur, the system uses... Calculate the average Phred mass fraction of the two mismatched bases, and multiply the average mass value by a scaling factor. This process generates a penalty value for that location; the higher the average quality, the larger the penalty value. Finally, the penalty values ​​calculated for all mismatch locations are summed to obtain the total quality penalty. .

[0060] Then, calculate the current two second UMI sequences according to the following formula. and Final similarity :

[0061] (4)

[0062] in, The purpose of this is to ensure that the similarity does not become negative. When the quality penalty is too large and exceeds the basic similarity, the two second UMIs will be considered completely dissimilar.

[0063] Steps one and two are performed for each pair of second UMIs, ultimately yielding the final similarity between all pairs of second UMIs. Unless otherwise specified, the similarity calculated in subsequent embodiments refers to the final similarity.

[0064] S220. Based on the distribution of similarity, determine the similarity threshold used to distinguish whether the current batch is similar, and correct the similarity threshold based on the global average quality of the second UMI of the current batch.

[0065] This embodiment determines an adaptive similarity threshold for the current batch based on the similarity between the second UMIs. That is, to determine whether two second UMIs are "similar" enough to be merged, the similarity threshold is dynamically calculated based on the data characteristics of the current batch. Optionally, the similarity threshold for whether two second UMI sequences should be clustered together can be denoted as... This threshold depends on the data quality and error patterns of the current sequencing batch.

[0066] Optionally, this embodiment proposes a hybrid decision-making method to determine the similarity threshold. The specific principle is as follows: typically, all second UMIs in the same batch... The distribution exhibits a bimodal pattern: the left peak consists of numerous truly unrelated second UMI pairs, which typically have low similarity; the right peak consists of a small number of "child UMIs" (both referring to second UMIs) generated by sequencing errors and highly similar to the "parent UMI," exhibiting high similarity. Threshold It is located in a valley between two mountain peaks. Therefore, this embodiment uses the following three methods to generate the threshold. :

[0067] Method 1: Generate a similarity density curve based on the similarity distribution; and determine the trough position between the two peaks in the density curve as a similarity threshold. Optionally, this can be applied to all... By applying the kernel density estimation method, a smooth distribution density function curve is generated. In the density curve A derivative-based peak detection algorithm is used to identify local maxima on the density curve, and the two highest peaks are defined as the left main peak. And the right main peak The horizontal axis (similarity value) is defined as follows: and If, after searching for the dominant peak, the data distribution is found to be not bimodal (e.g., unimodal, multimodal, or no peak), it will directly return None, indicating that subsequent decisions will not consider the kernel density estimation strategy; then, trough location will be performed to determine a similarity threshold. The calculation formula is:

[0068] (5)

[0069] in, This represents a similarity threshold. It is calculated by extracting the density curve portion between two main peaks, inverting the curve, and then using a peak detection algorithm to determine the coordinates of the "peak" (corresponding to the trough of the original curve) of the inverted curve. This trough position is the similarity threshold obtained using the KDE method.

[0070] Furthermore, peak identification can be achieved by calculating the zero point of the first derivative and the negative value of the second derivative of the density function. A point is identified as a peak if it meets the following conditions: (1) the first derivative is zero or close to zero; (2) the second derivative is negative (a characteristic of concave functions); and (3) the function value is higher than the local maximum value of the adjacent region. The same mathematical principle is used for trough location, which is based on the inverted density curve (i.e., Peak detection is applied, and the peak value detected corresponds to the trough position of the original curve.

[0071] Method 2: Using the percentile method, determine another similarity threshold based on the expected percentile of the similarity threshold in the distribution. Optionally, the expected percentile can be the 65th percentile, which can be used to represent all... Sort from lowest to highest, and take the fixed 65th percentile. The calculation formula is as follows:

[0072] (6)

[0073] in, This represents another similarity threshold.

[0074] Method 3: Set a safety lower bound. Preset an absolute minimum similarity threshold. Exemplary 0.62.

[0075] The three methods correspond to three similarity thresholds. For ease of distinction and description, this embodiment uses the threshold obtained using method one. This is called the first similarity threshold, and the threshold obtained using method two will be used. This is called the second similarity threshold. Therefore, this embodiment first determines a basic similarity threshold based on these three similarity thresholds. .

[0076] In one specific implementation, the larger value between the first similarity threshold and the second similarity threshold is first taken as the base similarity threshold. That is:

[0077] (7)

[0078] If the KDE method fails, then Considered invalid. Direct value .

[0079] Then, the basic similarity threshold and the safety lower limit are taken. The larger value in the range is used as the final basic similarity threshold. This threshold is used to initially distinguish whether the second UMI is similar. That is:

[0080] (8)

[0081] After obtaining the final basic similarity threshold, a global quality adjustment factor is calculated to adjust the final threshold based on the global average quality of the current batch of data to ensure that it better matches the characteristics of the current batch of data.

[0082] In one specific implementation, firstly, the following operation can be performed for each second UMI sequence: the average quality of each base position in the current second UMI sequence across each read segment is averaged to obtain the average quality of the current second UMI sequence. Optionally, the average quality of all base positions in the k-th second UMI sequence calculated using formula (1) is then used... Taking the average, we obtain the average quality of the k-th second UMI sequence. .

[0083] Then, based on the gene abundance of each second UMI sequence, a weighted average of the average quality of each second UMI sequence is calculated to obtain the global average quality of each second UMI sequence. Optionally, the abundance of each second UMI maintained in S210 can be used. The global average quality is obtained by weighting the average quality of all second-level UMIs in the current batch. :

[0084] (9)

[0085] Finally, the similarity threshold is adjusted based on the global average quality. A higher global average quality results in a higher adjusted similarity threshold (i.e., a more stringent condition). Optionally, the global average quality can be compared with the reference quality. By comparison, the global quality adjustment factor is obtained. :

[0086] (10)

[0087] in, The default value is 30. This represents the adjustment factor, which defaults to 0.01. The final threshold correction is then:

[0088] (11)

[0089] Unless otherwise specified, the similarity thresholds used below refer to the corrected similarity thresholds described here. .

[0090] S230. Using each second UMI as a node, construct a graph with the node relationship with similarity higher than the corrected similarity threshold as an edge, and prune according to the abundance difference of the second UMI at both ends of each edge, and divide each closed subgraph in the pruned graph into a UMI group.

[0091] This step performs adaptive UMI correction and clustering. It uses a graph theory-based approach combined with a directional pruning strategy to cluster UMIs. By adaptively determining the clustering threshold, it distinguishes true UMI sequences from sequencing errors. It also prunes based on the abundance of UMIs, grouping reads originating from the same initial molecule together.

[0092] Specifically, graph-based clustering and targeted pruning aim to group together UMI sequences associated due to sequencing errors, while avoiding the merging of independent, genuine UMI errors. This is achieved by networking UMI relationships and incorporating an adaptive similarity threshold. Unreliable connection edges are filtered out using UMI abundance.

[0093] In one specific embodiment, the process may include the following steps:

[0094] Step 1: Initialize the graph. Specifically, initialize the graph framework G. ,in It is a set of nodes. It is a set of edges. This embodiment transforms each unique UMI molecular entity into a node in a graph. If the set of all unique UMI sequences (i.e., the second UMI) in the current batch is considered... Expressed as a formula: For Each of the second UMIs Create a corresponding node in the graph. All nodes constitute a node set. Each of the nodes Corresponding to a second UMI sequence Then, an abundance attribute is assigned to each node. .

[0095] Step 2: Construct the initial candidate relationship graph. After initializing the graph, the calculated adaptive similarity threshold can be used. Connect the nodes to form the initial structure of the graph. ,in, This represents the initial edge structure. Specifically, for any two distinct nodes... and Between them, if their corresponding second UMI similarity Greater than the preset If the condition is met, then an edge is created. The corresponding formula is:

[0096] (12)

[0097] At this time It contains all possible connections, but it is mixed with real "parent-child" errors and fake "similar neighbor" relationships.

[0098] Step 3: Targeted pruning of the initial candidate relationship graph. This step is one of the core inventive points of this embodiment. By reviewing each candidate edge and combining it with abundance information to determine its rationality, the graph is optimized. In this embodiment, one edge... High reliability requires the following condition to be met: the abundance of one node must be significantly higher than that of another node; for example, the abundance of the second UMI at one end of an edge must be a multiple of the abundance of the second UMI at the other end. Optionally, this condition can also be expressed as follows:

[0099] (12)

[0100] in, and These are nodes and The abundance. It should be noted here that abundance, from the perspective of the second UMI, is denoted as... From the perspective of graph nodes, the same abundance is denoted as The two satisfy the identity relation This refers to the second UMI in the dataset being equal to Total number of reads. It is a multiplier factor, with a default value of 2.2, which requires that the abundance of one UMI must be more than 2.2 times that of another. It is a constant, with a default value of 1, used to provide buffering when dealing with low-abundance UMI.

[0101] Based on this condition, each edge in the graph is traversed and checked; any edge that does not meet the condition is removed from the graph. Therefore, the final set of edges after pruning is... for:

[0102] (13)

[0103] Therefore, the final result is shown in the figure. This step corresponds to the similarity threshold adjustment in S220 based on the overall data quality of the current batch. The essence of the similarity threshold adjustment in S220 is to relax the screening conditions, extract more second UMI candidates into the graph, and avoid insufficient candidate information due to inappropriate thresholds. Here, the strict screening based on edge credibility tightens the screening conditions and ensures the screening quality on the basis of sufficient candidate information.

[0104] Step 4: Extract UMI groups based on the pruned graph. After the previous pruning step, this embodiment extracts the UMI groups from the final graph. Find all connected components (i.e., subgraphs in the graph where there is a path between any two nodes and the subgraph is not connected to any other part of the graph by edges) in the final graph. All UMI group sets This can be expressed as:

[0105] (14)

[0106] Where Q represents the number of UMI groups. For each group including component groups (q=1,2,…,Q), select the node with the highest abundance in the group as the main UMI of the group, and collect the reads of all the second UMIs in the group for the next stage of consensus sequence generation.

[0107] S240. Perform the following operations for each UMI group to identify the consensus gene sequence and consensus quality of each UMI group:

[0108] Step 1: For each UMI group: pad the reads containing the second UMI within the group to a uniform length and align them into an alignment matrix with each read as a row and each base position in the read as a column.

[0109] To maintain high efficiency when processing a large number of UMI groups, this embodiment employs a heuristic progressive alignment strategy, converting a UMI group into an alignment matrix. Simultaneously, homopolymer regions are identified in the reference sequences selected from the UMI groups for subsequent comprehensive weight calculations.

[0110] In one specific implementation, the set of all reads within the current UMI group is first defined as... Where N is the total number of reads in the group. Then, the longest sequence (read) from the current UMI group is selected as the reference sequence. ,in, Indicates reading segment The length of this sequence, where n represents the segment index (n = 1, 2, ..., N), determines the width of the final alignment matrix. .

[0111] Simultaneously, homopolymer regions in the reference sequence are identified. Optionally, a subsequence is defined as a homopolymer if:

[0112] Condition 1: Uniformity must be satisfied; all bases in the identified region must be identical. That is, bases must exist. This makes a continuous sequence of positions in the reference sequence place (among which) ) bases ,in, and These represent the start and end positions of this continuous sequence of positions, respectively.

[0113] Condition 2, Length Condition: Length of the identified homopolymer region The threshold must be reached or exceeded. ,Right now: , The default value is 3.

[0114] Condition 3, Maximality Condition: The identified region cannot extend to either end; that is, it is the longest continuous identical base region, ensuring that a large homopolymer region is not mistakenly identified as multiple small overlapping regions.

[0115] Identified homopolymer regions The positions of each base are stored for use in subsequent comprehensive weight calculations: .

[0116] Then, define the sequence alignment transformation function. This function is used to convert any sequence Convert to a sequence aligned with the reference sequence, with a length of [length missing]. new sequence of The rule for this change is to left-align the sequence and fill it with empty spaces. For example, for any sequence... The transformed sequence Defined as: It is a line of length The sequence, its first The bases at each base position Determined by the following piecewise function:

[0117] (15)

[0118] That is, location index In sequence If the length is within the range, then the base at position m in the new sequence is equal to the base at position m in the original sequence. The position of the base, if indexed Exceeding the sequence The length of the new sequence is then... The position is then filled with the empty character "-".

[0119] For sets Each sequence in All apply transformation functions Ultimately, we can obtain a new sequence with all sequences of length 1. set By arranging all the new sequences sequentially, with each new sequence as a row, we can obtain a dimension... The comparison matrix .

[0120] Step 2: Calculate the quality weight, position weight, homopolymer weight, and coverage weight of each row in the same column of the alignment matrix, and determine the consensus bases and consensus quality of each column based on each weight. The consensus bases of each column are arranged to form the consensus gene sequence of the current UMI group.

[0121] This step generates a consensus sequence based on quality weighting of the alignment matrix. Specifically, this step combines sequence features to calculate a comprehensive weight score and outputs a consensus sequence through voting, thereby generating a consensus sequence that best represents the original sequence from all associated reads. Simultaneously, to ensure the reliability of the generated consensus sequence, its quality is also output, providing information for downstream applications of the formula sequence.

[0122] To achieve the above objective, in one specific embodiment, the comparison matrix can be targeted. For each column in the table, perform the following calculations:

[0123] S2-1. Initialization preparation based on the alignment matrix. Let the current column be the m-th column. First, count the number of non-"-" bases in this column, as the coverage of that column. Then, iterate through each read (i.e., each row in the matrix) in that column. ), to obtain the Phred mass fraction at its base positions. If the corresponding position is "-" (i.e., a gap), or If the value is lower than the preset quality value (default is 10), the base position will be considered to have a low contribution and will be skipped in subsequent calculations.

[0124] S2-2. For each confirmed valid base position in the current column, calculate the overall weight of each base position sequentially. The overall weight includes quality weight, position weight, homopolymer weight, and coverage weight. The calculation methods for each weight are as follows:

[0125] Regarding quality weights: The quality weights of each row in the current column can be calculated using the following formula:

[0126] (16)

[0127] in, This represents the mass of the base position in the nth row of the mth column. This represents the mass weight of the base position in the nth row of the mth column. This represents the maximum quality value among all rows in the m-th column, used to standardize the Phred score. ; This represents the nonlinear mapping parameter, used to nonlinearly map the Phred quality score to the weight range of [0.1,1]. Optionally, k is fixed at 10. The processing significance is that bases below Phred20 are quickly suppressed to near the lower weight limit of 0.1, while bases above Phred20 are quickly raised to the upper weight limit of 1.0, ensuring that reliable and unreliable bases can be significantly distinguished near Phred 20.

[0128] In the quality weighting, a non-linear mapping based on the sigmoid function is used to achieve higher weights for bases with higher quality scores.

[0129] Regarding position weights: The position weights of each row in the current column can be calculated using the following formula:

[0130] (17)

[0131] in, This represents the positional weight of the base in the nth row of the mth column. This indicates the length of the nth line before it is padded; This represents the attenuation coefficient, used to control the rate at which the weight decreases as it deviates from the center.

[0132] Position weights are used to reduce the influence of bases located at the ends of the sequencing read (because the error rate of sequences at the ends of long reads is usually higher), and to give higher weights to bases in the central region.

[0133] Regarding homopolymer weights: The homopolymer weights for each row in the current column can be calculated using the following formula:

[0134] (18)

[0135] in, This represents the homopolymer weight at the base position in the nth row of the mth column. This indicates the homopolymer region of the reference sequence; h represents the penalty factor, which is a positive number less than 1; when a base is located in a homopolymer region, its weight is multiplied by this factor to reduce its influence, and the default h is 0.6.

[0136] Regarding coverage weights: The coverage weights of each row in the current column can be calculated using the following formula:

[0137] (19)

[0138] in, This represents the coverage weight of each row in the m-th column; represents the base coverage of column m, which is the total number of valid (non-gap, high-quality) bases in that column.

[0139] Coverage weights are used to give a slight weight bonus to locations with higher coverage.

[0140] After calculating the weights for each of the above items, multiplying the weights for each base position together yields the overall weight for each base position:

[0141] (20)

[0142] in, This represents the overall weight of the base position in the nth row and mth column of the alignment matrix.

[0143] S2-3. Perform weighted majority voting based on the comprehensive weight of each base position. Optionally, firstly, for each base, sum the comprehensive weights of all rows containing the current base in the same column to obtain the quality score of the current base in that column. Optionally, the bases can be... The quality score in the current column can be expressed as: ,but:

[0144] (twenty one)

[0145] in, This is an indicator function. When the nth read segment contains the base in the mth column... When it is b, When the base in the m-th column of the nth read segment is not b, .

[0146] Then, the base with the highest overall weight in the current column is selected as the consensus base for that column. Optionally, the base with the highest quality score in the current column can also be selected. It can be represented as follows:

[0147] (twenty two)

[0148] Next, based on the overall weight of the consensus bases in each row of the current column, a weighted average of the consensus bases in each row is calculated to obtain the consensus quality of the current column. Optionally, the specific formula is as follows:

[0149] (twenty three)

[0150] in, This represents the consensus quality of the m-th column.

[0151] After performing operations S2-1 to S2-3 on all columns, each column corresponds to a consensus base and a consensus quality. Arranging the consensus bases of all columns sequentially yields the complete consensus sequence of the current UMI group; arranging the consensus quality of all columns sequentially reflects the quality of the consensus sequence. Consensus quality provides a confidence metric for downstream applications, helping to filter out low-quality consensus sequences or mask low-quality bases.

[0152] Similarly, after performing steps one and two of S240 on all UMI groups, the consensus sequence and consensus quality of each UMI group can be obtained. After performing S10-S20 on all batches of read segments, the UMI sequence correction and consensus sequence generation for each batch can be adaptively completed based on the data quality of each batch.

[0153] The entire process described above can also be combined Figure 2 Please refer to the flowchart shown for comprehension.

[0154] In summary, this embodiment provides a method for long-read UMI sequence correction and consensus sequence generation, aiming to overcome the shortcomings of existing technologies, such as inability to handle variable-length UMI sequences, insufficient information utilization, lack of adaptability, and failure to improve sequence quality when processing high-error-rate sequencing data. The method of this embodiment can handle variable-length UMI sequences and, by introducing base quality-weighted distance, adaptive threshold calculation, and downstream consensus sequence generation, achieves ultra-high precision correction of UMI sequences and deep error correction of single-molecule sequences. Specifically, the method of this embodiment can achieve the following beneficial effects:

[0155] 1. Introducing a base quality (Q-score) weighted calculation for UMI distance. This embodiment no longer uses simple edit distance, but instead calculates a "weighted distance" between two UMIs using a comprehensive model. This approach considers base substitutions, insertions, or deletions simultaneously, assigning a weight based on its Q-score to each differential base. A mismatch of a high-quality base will incur a much higher "distance penalty" than a low-quality base mismatch, enabling the algorithm to accurately distinguish between genuine UMI differences and high-probability sequencing errors.

[0156] 2. Adaptive calculation of clustering thresholds. Addressing the quality differences across datasets, this embodiment abandons a fixed merging threshold. Before clustering the UMIs of each cell (or gene), this method first rapidly analyzes the abundance distribution, sequence similarity distribution, and average base quality of the UMIs within that subset. Based on these characteristics, a "weighted distance" threshold best suited to the current data quality is dynamically calculated, thereby achieving "personalized" processing for each data subset and maximizing the accuracy of the correction.

[0157] 3. Achieving high-fidelity consensus sequence generation. For each UMI cluster formed after correction and merging, the method in this embodiment further utilizes the complete sequence of all reads within the cluster and the corresponding base quality information. By calling multiple sequence alignment to arrange these reads, and based on the depth, consistency, and average quality of each base site, a "high-fidelity consensus sequence" with a theoretically much lower error rate than any single original read is calculated and generated.

[0158] Compared with the prior art, the method of this embodiment has the following advantages:

[0159] 1. This method offers variable-length sequence alignment capability and higher correction accuracy. It overcomes the limitations of existing technologies based on fixed-length alignment by achieving indel-aware variable-length sequence alignment, thus solving the UMI correction problem in the era of long-read sequencing. When processing standard Illumina next-generation sequencing data, the UMI correction accuracy of this method is comparable to existing technologies (such as UMI-tools), both reaching over 98%. Especially when processing long-read data with high indel error rates, tests based on the simulated ONT platform show that the UMI correction accuracy of this method reaches 94.5%, while the baseline method, due to its fixed-length alignment model's inability to handle variable-length UMI, experiences a sharp drop in accuracy to 82.0%, a decisive difference of 12.5 percentage points. This breakthrough enables the molecular quantification accuracy of long-read single-cell sequencing to reach a practical level for the first time.

[0160] 2. High-fidelity sequence error correction and genetic disease feature detection capabilities. This embodiment's method, by introducing a multidimensional weighted consensus model, transforms UMI correction from a simple "molecular counting tool" into a "sequence error correction and variant detection tool," thus improving sequence fidelity. On Illumina data, the consensus sequence identity generated by this embodiment's method is almost identical to the baseline method (both exceeding 99.8%). However, when simulating ONT sequencing data, this embodiment's method reduces the original error rate from 8% to below 1.2%, generating a consensus sequence identity as high as 98.8%. Particularly noteworthy is the error correction capability in homopolymer regions: homopolymer regions (such as consecutive AAA or TTT repeats) are a "major problem area" in traditional sequencing technologies and are also hotspots for pathogenic mutations in many genetic diseases. This embodiment's method significantly reduces the InDel error rate in homopolymer regions from over 20% in the baseline method to 4%, achieving a 5-fold performance improvement. This high-fidelity error correction capability is of great significance for the detection of clinical genetic / rare diseases and the characterization of tumor somatic mutations, and can provide an accurate molecular detection data foundation.

[0161] Furthermore, to verify the technical effectiveness of the above method, this embodiment also designed a simulated data generation system based on the characteristics of real sequencing platforms. This system can accurately reproduce the error characteristics of different sequencing platforms, including substitution errors, insertion / deletion errors (Indel), high error rates in homopolymer regions, and error distribution patterns related to base mass fractions.

[0162] The core parameter settings of the system include:

[0163] 1. Error rate model:

[0164] PCR amplification error rate: 1×10-5 (0.001%), consistent with the typical performance of high-fidelity polymerase;

[0165] Sequencing error rate: 3% (0.03), based on actual test data from the Oxford Nanopore platform;

[0166] Context relevance: Error rate improved by 1.5 times in homopolymer regions, 1.2 times in GC-rich regions, and decreased to 0.9 times in AT-rich regions;

[0167] Error type distribution: Supports pure replacement mode (Illumina compatible) or mixed mode (including Indel).

[0168] 2. Indel error characteristics:

[0169] Indel percentage of total errors: adjustable from 0% (pure replacement) to 80% (high indel scenario), default is 80% to simulate long read length platforms;

[0170] Insertion / deletion ratio: 1:1 (50 / 50 allocation), consistent with the measured distribution of Nanopore and PacBio;

[0171] Length distribution: mainly single-base indels (~90%), with a small number of multi-base indels;

[0172] Location preference: The incidence of indels in homopolymer regions is significantly increased.

[0173] 3. Homopolymer characteristics:

[0174] Length distribution: 2-6 bases, following an exponential decay distribution with probabilities of [0.4, 0.3, 0.2, 0.08, 0.02];

[0175] Detection threshold: The minimum length is set to 2bp to ensure that all potentially high-risk areas are captured;

[0176] Frequency of occurrence: Approximately 15% of UMI sequences contain at least one homopolymer region;

[0177] Base composition: A / T homopolymer accounts for approximately 60%, and G / C homopolymer accounts for approximately 40%.

[0178] 4. Base mass fraction distribution:

[0179] Correct bases: Q30-Q41 (Phred mass fraction), corresponding to an error probability of < 0.1%;

[0180] Incorrect bases: Q5-Q15, corresponding to an error probability range of 3%-32%;

[0181] Quality-error correlation: realistically simulating the quality control logic of a sequencer by assigning low quality fractions to incorrect bases;

[0182] The quality distribution conforms to the statistical characteristics of real sequencing data.

[0183] To fairly assess the advantages of this invention over existing technologies, this embodiment designs a dual-mode data generation strategy:

[0184] Experimental data: Includes a complete error model (replacement + indel + homopolymer + quality weighting), with an indel ratio of up to 80%, truly reflecting the complex error patterns of long-read sequencing.

[0185] Control group data: The --no-indels parameter generates equal-length UMI sequences containing only substitution errors, which is compatible with the processing premise of UMI-tools.

[0186] Baseline method: UMI-tools directional algorithm, using a fixed Hamming distance threshold of 1, which is the standard configuration and best practice of this tool.

[0187] The key significance of this comparative design is that the UMI-tools algorithm, based on fixed Hamming distance, cannot handle variable-length sequences (the UMI length changes after indels), thus its performance drops significantly on experimental data. This invention, by introducing multiple sequence alignment (MSA) and a quality-weighted mechanism, can effectively handle variable-length sequences and complex error patterns, thereby maintaining stable high performance even in high-indel scenarios.

[0188] The final verification results are as follows: Compared with existing technologies (such as UMI-tools), the method in this embodiment has achieved improvements in UMI deduplication accuracy, consensus sequence fidelity, and downstream analysis applications, specifically reflected in:

[0189] 1. Improved UMI deduplication accuracy and adaptive error correction. This embodiment overcomes the limitations of existing technologies (such as UMI-tools) that rely on a fixed Hamming distance threshold (e.g., distance = 1) for UMI grouping. Existing technologies are prone to over-correction (incorrectly merging true UMIs) or under-correction (failing to merge sequencing errors) when dealing with libraries of varying sequencing error rates and complexities. This embodiment avoids the bias caused by a fixed threshold. Ultimately, the method in this embodiment achieves more accurate UMI grouping than UMI-tools under different library quality and sequencing depth conditions, improving the accuracy of single-molecule counting, especially when processing data with complex error patterns or high-depth sequencing data. Figure 3As can be seen, in a simulation scenario where the Indel error rate is as high as 70%, the UMI correction accuracy of the method in this embodiment can still be maintained at 94.5%, while the accuracy of the benchmark method drops to 82.0% due to its model limitations, showing a decisive performance gap between the two.

[0190] 2. High-fidelity consensus sequence generation and error suppression for complex sequences. After obtaining accurate UMI grouping, existing technologies such as UMI-tools typically use simple majority voting to generate consensus sequences. This method cannot effectively handle systematic biases and errors during sequencing. This embodiment overcomes this limitation by proposing a multi-dimensional comprehensive weighting model to construct high-fidelity consensus sequences. This comprehensive weighting model intelligently suppresses sequencing noise from multiple physical and statistical dimensions. The final result is that... Figure 4 Even in areas where traditional methods are most problematic, such as low coverage, read length ends, or the presence of homopolymers, this invention can generate consensus sequences with higher fidelity. Figure 5 In processing PacBio or ONT platform data with high Indel error rates, the consensus sequences generated by the method in this embodiment have an InDel error rate in the homopolymer region that is 5-6 times lower than the benchmark method. This provides a data foundation superior to existing technologies for subsequent analysis such as variant detection and transcript splicing, and effectively reduces the false positive rate.

[0191] 3. Enhanced downstream analytical capabilities and sensitivity for detecting rare events. Due to the more accurate UMI quantification and higher-fidelity consensus sequences of the method in this embodiment, the depth and breadth of downstream bioinformatics analysis are directly improved. Specifically:

[0192] In terms of quantitative expression: it allows for more accurate quantification of gene / transcript expression, especially improving the detection sensitivity and quantification accuracy for low-expression genes. Figure 6 Simulation data show that the method in this embodiment consistently and stably provides lower quantification bias and higher detection recall than the baseline method. This improvement is particularly critical when gene expression abundance is extremely low, enabling more reliable identification and quantification of marker genes in rare cell subpopulations.

[0193] In variant detection: High-fidelity consensus sequences effectively filter out a large number of random sequencing errors, significantly improving the signal-to-noise ratio of somatic mutations or RNA editing events at the single-cell level. This makes it possible to detect low-frequency (subclonal) variants within tumors and discover key mutations in rare cell subpopulations, which is often difficult to achieve with traditional tools such as UMI-tools due to excessive background noise. Figure 7Simulation experiments demonstrate that the method in this embodiment can reliably detect real biological variants with frequencies as low as 30-40% in sequencing reads (detection rate >94%), while the traditional majority voting method has a detection rate of less than 30% under the same conditions. This invention significantly lowers the detection threshold for rare variants while controlling the false positive error rate to below 0.6%, providing crucial technical support for precision medicine research.

[0194] 4. Quantitative Bias Control and Molecular Count Stability. One of the core values ​​of UMI technology is achieving absolute quantification of RNA molecules, that is, accurately inferring the original molecule number by counting the number of different UMI families. However, existing technologies (such as UMI-tools) suffer from serious systematic quantitative bias when processing sequencing data with high indel error rates, a problem particularly prominent in third-generation sequencing platforms (PacBio, ONT). The root cause is that UMI-tools' clustering strategy based on a fixed Hamming distance can only handle base substitution errors and cannot correctly identify sequence length changes caused by indel errors. When the indel error rate in sequencing data increases, multiple UMI sequence copies originating from the same original molecule will have length differences due to insertion or deletion errors, resulting in the inability to calculate the Hamming distance or an abnormally large value. This causes sequences that should cluster into one UMI family to be incorrectly split into multiple "false families," resulting in a systematic overestimation of molecular counts, or "false overestimation."

[0195] The method in this embodiment fundamentally solves this quantitative bias problem by introducing an Indel-based perceptual distance metric and a quality score weighting strategy based on multiple sequence alignment. Figure 8 The simulation clearly demonstrates the comparison of the quantitative accuracy of this method as the Indel error rate increases from 0% to 10%.

[0196] Regarding the quantitative stability of Qscore-UMI: the observed / expected ratio generated by the method in this embodiment remains consistently within the range of 1.0 ± 0.05, indicating that the molecule counting results remain accurate regardless of changes in the Indel error rate, and there is no systematic bias. This is due to the fact that multiple sequence alignment can correctly handle variable-length sequence alignment, incorporating the sequence length changes caused by Indel errors into the similarity calculation. Combined with the mass fraction weighting mechanism, it can accurately identify homologous UMI sequences and avoid spurious segmentation.

[0197] In the UMI-tools method, systematic overestimation occurs: In contrast, the baseline method's quantification ratio increases exponentially with the Indel error rate. When the Indel error rate reaches 10% (close to the typical level of ONT R9 chemistry), the observed / expected ratio is as high as approximately 2.3, meaning that the molecular count is overestimated by 130%. This artificial overestimation directly leads to false positives in differential expression analysis, potentially causing misjudgments of rare cell types or low-abundance markers in clinical practice.

[0198] It should be noted that the information and data involved in this application are all authorized by the users or fully authorized by all parties, and the collection, use and processing of the relevant data must comply with the relevant laws, regulations and standards of the relevant countries and regions, and corresponding operation portals are provided for users to choose to authorize or refuse.

[0199] Figure 9 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention, such as... Figure 9 As shown, the device includes a processor 60, a memory 61, an input device 62, and an output device 63; the number of processors 60 in the device can be one or more. Figure 9 Taking a processor 60 as an example; the processor 60, memory 61, input device 62, and output device 63 in the device can be connected via a bus or other means. Figure 9 Taking the example of a connection between China and Israel via a bus.

[0200] The memory 61, as a computer-readable storage medium, can be used to store software programs, computer-executable programs, and modules, such as the program instructions / modules corresponding to the long read length UMI sequence correction and consensus sequence generation method in this embodiment of the invention. The processor 60 executes various functional applications and data processing of the device by running the software programs, instructions, and modules stored in the memory 61, thereby implementing the aforementioned long read length UMI sequence correction and consensus sequence generation method.

[0201] The memory 61 may primarily include a program storage area and a data storage area. The program storage area may store the operating system and at least one application program required for a given function; the data storage area may store data created based on terminal usage. Furthermore, the memory 61 may include high-speed random access memory and non-volatile memory, such as at least one disk storage device, flash memory, or other non-volatile solid-state storage device. In some instances, the memory 61 may further include memory remotely located relative to the processor 60, which can be connected to the device via a network. Examples of such networks include, but are not limited to, the Internet, intranets, local area networks, mobile communication networks, and combinations thereof.

[0202] Input device 62 can be used to receive input digital or character information, and to generate key signal inputs related to user settings and function control of the device. Output device 63 may include display devices such as a display screen.

[0203] This invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the long read length UMI sequence correction and consensus sequence generation method of any embodiment.

[0204] The computer storage medium of this invention can be any combination of one or more computer-readable media. A computer-readable medium can be a computer-readable signal medium or a computer-readable storage medium. A computer-readable storage medium can be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of computer-readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable computer disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage device, magnetic storage device, or any suitable combination thereof. In this document, a computer-readable storage medium can be any tangible medium that contains or stores a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0205] Computer-readable signal media may include data signals propagated in baseband or as part of a carrier wave, carrying computer-readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. Computer-readable signal media may also be any computer-readable medium other than computer-readable storage media, capable of sending, propagating, or transmitting programs for use by or in connection with an instruction execution system, apparatus, or device.

[0206] Program code contained on a computer-readable medium may be transmitted using any suitable medium, including but not limited to wireless, wire, optical fiber, RF, etc., or any suitable combination thereof.

[0207] Computer program code for performing the operations of this invention can be written in one or more programming languages ​​or a combination thereof. Programming languages ​​include object-oriented programming languages—such as Java, Smalltalk, and C++—as well as conventional procedural programming languages—such as C or similar programming languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0208] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; 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 or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A method for long-read UMI sequence correction and consensus sequence generation, characterized in that, include: Obtain at least one batch of read segments to be processed; Perform the following operations for each batch: S210. Extract the first UMI to be corrected from each read segment in the current batch; remove duplicates from each first UMI to obtain each second UMI after deduplication; and calculate the similarity between each pair of second UMIs based on the base difference between each pair of second UMIs and the average mass of each base position. S220. Based on the distribution of similarity, determine the similarity threshold used to distinguish whether the current batch is similar, and correct the similarity threshold based on the global average quality of the second UMI of the current batch. S230. Using each second UMI as a node, construct a graph with the node relationship with similarity higher than the corrected similarity threshold as an edge, and prune according to the abundance difference of the second UMI at both ends of each edge, and divide each closed subgraph in the pruned graph into each UMI group. S240. For each UMI group: fill the reads containing each second UMI within the group with a uniform length and align them into an alignment matrix with each read as a row and each base position in the read as a column; calculate the quality weight, position weight, homopolymer weight, and coverage weight of each row in the same column of the alignment matrix, and determine the consensus bases and consensus quality of each column according to each weight, and arrange the consensus bases of each column into the consensus gene sequence of the current UMI group.

2. The method according to claim 1, characterized in that, The calculation of the similarity between two second UMIs based on the base difference and the average quality of each base position includes: For every two second UMIs: The similarity between the two current second UMIs is determined by the number of positions with the same bases. The more positions with the same bases, the greater the similarity. The similarity is corrected based on the quality of the two second UMIs with different base positions, to obtain the final similarity between the two second UMIs.

3. The method according to claim 2, characterized in that, The step of correcting the similarity based on the quality of the two second UMIs with different base positions to obtain the final similarity between the two second UMIs includes: Calculate the current two UMI sequences using the following formula. and Final similarity : in, express and Similarity; express and Quality penalty value between; Indicates the uniform length of a UMI sequence; and They represent The Middle Bases at each base position; and They represent The Middle Average mass of each base position; Indicates the indicator function, when the first... When the bases at each base position are different When the first When the bases at each base position are the same ; This represents the scaling factor, used to control the overall impact strength of the quality penalty.

4. The method according to claim 1, characterized in that, The step of determining the similarity threshold for distinguishing whether a batch is similar based on the distribution of similarity includes: Based on the distribution of similarity, a density curve of similarity is generated; and the position of the trough between the two peaks in the density curve is determined as the first similarity threshold. A second similarity threshold is determined based on the expected percentile of the similarity threshold in the distribution; The larger of the first similarity threshold and the second similarity threshold is taken as the intermediate threshold; then the larger of the intermediate threshold and the preset threshold lower limit is taken as the similarity threshold used to distinguish whether they are similar.

5. The method according to claim 1, characterized in that, The step of correcting the similarity threshold based on the global average quality of the second UMI in the current batch includes: For each second UMI sequence: the average quality of the current second UMI sequence is obtained by averaging the average quality of each base position in each read segment; Based on the abundance of each second UMI sequence, the average quality of each second UMI sequence is weighted and averaged to obtain the global average quality of the second UMI sequences in the current batch. The similarity threshold is corrected based on the global average quality, wherein the higher the global average quality, the higher the corrected similarity threshold.

6. The method according to claim 1, characterized in that, The pruning based on the abundance difference of the second UMI at both ends of each side includes: Prune edges that do not meet the following condition: The abundance of the second UMI at one end of the current edge is a multiple of the abundance of the second UMI at the other end.

7. The method according to claim 1, characterized in that, The calculation of the quality weight, position weight, homopolymer weight, and coverage weight for each row in the same column of the alignment matrix includes: For each column of the alignment matrix, perform at least one of the following operations: Calculate the quality weight of each row in the current column using the following formula: in, This represents the mass of the base position in the nth row of the mth column. This represents the mass weight of the base position in the nth row of the mth column. This represents the maximum mass among all rows in the m-th column. Represents nonlinear mapping parameters; Calculate the position weight of each row in the current column using the following formula: ， in, This represents the positional weight of the base in the nth row of the mth column. This indicates the length of the nth line segment before it is filled. This represents the attenuation coefficient, used to control the rate at which the weight decreases as it deviates from the center; Calculate the homopolymer weights for each row in the current column using the following formula: in, This represents the homopolymer weight at the base position in the nth row of the mth column. The longest read segment in the alignment matrix represents the homopolymer region; h represents the penalty factor, where h is a positive number less than 1. Calculate the coverage weight of each row in the current column using the following formula: ， in, This represents the coverage weight of each row in the m-th column. represents the base coverage of column m.

8. The method according to claim 1, characterized in that, The determination of consensus bases and consensus quality for each column based on their respective weights includes: Multiply the weights of each row in the same column together to get the combined weight of each row; For each base, the combined weights of all rows containing the current base in the same column are summed to obtain the mass fraction of the current base in that column. The base with the highest overall weight in each column is selected as the consensus base for each column. Based on the comprehensive weight of the consensus bases in each column in each row, the consensus quality of each column is obtained by weighted averaging of the consensus bases in each column in each row.

9. An electronic device, characterized in that, include: One or more processors; Memory, used to store one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors implement the long read length UMI sequence correction and consensus sequence generation method according to any one of claims 1-8.

10. A computer-readable storage medium, characterized in that, It stores a computer program that, when executed by a processor, implements the long read length UMI sequence correction and consensus sequence generation method as described in any one of claims 1-8.

Images