A single-cell sequencing data feature sequence detection method suitable for multiple library structures

Abstract

The application provides a single-cell sequencing data feature sequence detection method suitable for multiple library structures. The method is adapted to multiple different requirements by defining a library structure and an output structure, and uses prefix tree-based fuzzy matching to improve processing speed.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and relates to methods for processing single-cell sequencing data, and more specifically to a method for detecting characteristic sequences in single-cell sequencing data that is adaptable to various library structures. Background Technology

[0002] Single-cell sequencing libraries typically include several specialized regions, such as barcode regions, capture sequence regions, and UMI regions. These regions do not contain insert fragments and are only used to distinguish cells, label UMIs, or are introduced during library construction; they have no special function. These regions need to be located and processed differently before actual analysis to adapt to the input requirements of different analytical workflows.

[0003] Currently, many software programs can process single-cell library transcriptome data, such as STARsolo, Cellranger, and UMItools. These programs generally rely on one premise and two methods to identify specific sequences. The premise is that the input library structure must conform to the software's presets; for example, Cellranger can only process 10x library structures, and STARsolo requires the R1 structure to be a cell barcode + unique molecular identifier (UMI). The identification methods are generally two: 1) Sequence identification based on specified positions, for example, identifying the first to 16 bases of R1 as barcodes and the 17th to 25th bases as UMIs; 2) Sequence identification based on a whitelist, meaning that if a sequence in a specified region matches a sequence in the whitelist, it is considered a successful identification. From the above analysis, it can be seen that these software programs support only a limited range of library structures, can only identify a few pre-set structures, have few custom options, and cannot customize the output format. Furthermore, the aforementioned software typically uses hash tables and Hamming distance for fuzzy matching of the whitelist, which can be time-consuming when high fault tolerance is required or the whitelist is large. Therefore, researchers need to develop targeted sequence processing algorithms based on their own library designs. Summary of the Invention

[0004] To address the issues of low flexibility and high processing time in existing technologies, this invention provides a method for detecting feature sequences in single-cell sequencing data that adapts to various library structures. By customizing the library structure and output structure, it can accommodate a variety of different needs, and uses fuzzy matching based on prefix trees to improve processing speed.

[0005] The technical solution adopted in this invention is: a method for detecting feature sequences in single-cell sequencing data that is adaptable to various library structures, comprising the following steps:

[0006] Prepare a library structure file to define the library structure attributes of the feature sequences in the reads to be filtered; the library structure attributes include at least the whitelisted sequences;

[0007] Prepare an output structure file to define the output pattern attributes of the feature sequences in the filtered reads; read each library structure attribute from the library structure file as a filtering condition attribute to form the filtering conditions; among them, the whitelist sequences in the library structure attributes are read through a prefix tree constructed for fuzzy sequence matching.

[0008] Read each output mode attribute from the output structure file and use it as the output structure;

[0009] Acquire single-cell sequencing data and read each read sequentially. For each read, query and match each feature sequence according to the filtering conditions, obtain the query results, and record the query sequence, matching sequence, matching start position, and matching end position as the filtering results of the corresponding query conditions. If the read passes the filtering of all feature sequences, it is considered to conform to the library structure and is regarded as a filtered read. Among them, for feature sequences configured with whitelist sequences, the prefix tree is used for fault-tolerant matching in the target region of the read, and the number of mismatches allowed does not exceed the mismatch tolerance.

[0010] For the filtered reads, the specified sequence feature content is output to the corresponding file according to the format defined in the output structure file.

[0011] Preferably, the library structure attributes include at least one or more of the following: feature name, position label, start position, end position, whitelist sequence, length, and mismatch tolerance.

[0012] Preferably, the output pattern attributes include at least one or more of the following: output position label, output order, feature name, retrieval method, output type, output sequence start position, and output sequence end position.

[0013] Preferably, the construction method of the prefix tree is determined based on whether the start position and the end position of the feature sequence are fixed: when the start position is not fixed but the end position is fixed, a reverse prefix tree is constructed; otherwise, a forward prefix tree is constructed.

[0014] Preferably, the method for constructing the prefix tree includes: initializing a node; for each sequence in the whitelist sequence, adding each base of the sequence as a child node in sequence; if the current base is the last base of the sequence, setting the sequence termination attribute of the current node to true and setting the sequence character attribute of the node to the full length of the sequence, otherwise setting the sequence termination attribute of the previous node to false and setting the current node as a new child node.

[0015] Preferably, the method for querying and matching each feature sequence sequentially according to the filtering conditions includes: calculating the actual start position and actual end position of the feature sequence in the corresponding reads based on the position label, start position, end position, and length defined in the library structure file, and obtaining the query sequence accordingly; if the feature sequence is configured with a whitelist sequence, then the prefix tree is used on the query sequence for fault-tolerant matching, and the number of mismatches allowed does not exceed the mismatch tolerance, and the query result is obtained; if the feature sequence is not configured with a whitelist sequence, then the actual start position and actual end position of the target region are directly used as the query result.

[0016] Preferably, the method for fault-tolerant matching using the prefix tree includes: if there is no child node in the current node that matches the base to be matched, then the search is performed by querying other child nodes of the current node and incrementing the mismatch count until a match is successful or the accumulated mismatch count exceeds the mismatch tolerance.

[0017] Preferably, the method further includes deduplicating the filtered reads: grouping the filtered reads according to their corresponding barcode and UMI; calculating a consensus sequence based on the target sequence for multiple reads in the same group; and outputting the consensus sequence as the representative sequence of the group when the support of the consensus sequence is greater than or equal to a preset threshold.

[0018] Preferably, the support is the ratio of the number of reads supporting the consensus sequence to the total number of reads in the group, and the preset threshold is 0.9.

[0019] Preferably, the start or end position of the feature sequence is defined by any of the following methods (1) to (3):

[0020] (1) Defined based on a fixed location;

[0021] (2) Define it relative to its positional attributes relative to the dependent structure and in combination with the offset;

[0022] (3) is defined as an uncertain position.

[0023] Preferably, the whitelist sequence is provided in any of the following (a) to (d):

[0024] (a) Directly specify one or more fixed sequences;

[0025] (b) Provided by referencing an external file containing a list of possible valid sequences;

[0026] (c) No specific whitelist sequence;

[0027] (d) Obtain the sequence from the header information of read using regular expressions.

[0028] Preferably, the library structure attribute also includes an N-policy, which defines whether a sequencing result of N is considered a mismatch.

[0029] Preferably, the method is applicable to single-cell transcriptome sequencing data and / or single-cell genome sequencing data.

[0030] The present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the method described herein.

[0031] The beneficial effects of this invention are:

[0032] This invention provides a highly flexible, fast-processing method for feature sequence detection in single-cell sequencing data that can simultaneously adapt to genomic libraries. This method achieves adaptation to various single-cell sequencing library structures through two customizable configuration files (library structure file and output structure file). During whitelist matching, the method uses a prefix tree-based fuzzy matching method to adapt to larger whitelists and higher mismatch capacity, thus enabling rapid multi-threaded processing of large amounts of data. Furthermore, for sequencing data such as genomics, this method can also perform deduplication based on barcode+UMI+seq to adapt to genomic libraries. Attached Figure Description

[0033] Figure 1 This is a simplified flowchart of a single-cell sequencing data feature sequence detection method adapted to various library structures, provided by an embodiment of the present invention.

[0034] Figure 2 This is a simplified flowchart of reading the library structure in a single-cell sequencing data feature sequence detection method adapted to multiple library structures provided in an embodiment of the present invention.

[0035] Figure 3 This is a simplified flowchart of the process for reading the output structure in a single-cell sequencing data feature sequence detection method adapted to multiple library structures, provided in an embodiment of the present invention.

[0036] Figure 4 This is a simplified flowchart of reading whitelist sequences using a prefix tree in a single-cell sequencing data feature sequence detection method adapted to multiple library structures provided in an embodiment of the present invention.

[0037] Figure 5This is a simplified flowchart of read filtering in a single-cell sequencing data feature sequence detection method adapted to multiple library structures provided in an embodiment of the present invention.

[0038] Figure 6 This is a simplified flowchart of read deduplication in a single-cell sequencing data feature sequence detection method adapted to multiple library structures provided in an embodiment of the present invention.

[0039] Figure 7 This is a simplified flowchart of the reads output in a single-cell sequencing data feature sequence detection method adapted to multiple library structures provided in an embodiment of the present invention. Detailed Implementation

[0040] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0041] This invention provides a method for detecting feature sequences in single-cell sequencing data that is compatible with various library structures. Its workflow can be summarized into six main stages: preparing the library structure, preparing the output structure, reading the library structure, reading the output structure, filtering reads, and outputting. In a preferred embodiment, the method also includes a read deduplication step.

[0042] refer to Figure 1This method first requires two configuration files: a library structure file defining the position of the feature sequences to be detected (such as barcodes and UMIs) in the reads, whitelisted sequences, and other library structure attributes; and an output structure file defining the output format. The library structure file is then read as filter conditions to obtain the filter criteria, and the output structure file is read as the output structure to obtain the output structure. During the library structure reading process, a prefix tree search model is constructed for each whitelisted feature sequence to support fast, error-tolerant matching. The end sequencing results of each read from the single-cell sequencing data are read sequentially. After checking the filter conditions in order, queries and matches are performed sequentially within the target regions of each feature sequence according to the order defined in the library structure. For whitelisted feature sequences, a fuzzy search is performed using the constructed prefix tree; for features without a whitelist, the sequence at the specified position is directly extracted. Reads that successfully match all features are considered filtered and retained. Based on the output structure configuration, the specified feature sequences are output to the corresponding files in a predetermined order and format, thus generating files adapted to different downstream analysis needs. In a preferred embodiment, the method further includes a deduplication step for reads of genomic data, which groups reads by a combination of barcode and UMI. For multiple reads in the same group, a consensus sequence is generated by comparing the similarity of their target sequences. Only when the support of the consensus sequence exceeds a preset threshold (e.g., 90%) is it considered a unique molecule and output.

[0043] In one specific embodiment, specific steps for preparing a library structure file are provided. The library structure file is used to define the library structure attributes of the feature sequences in the reads to be filtered. The library structure attributes include feature name, position label, start position, end position, whitelist sequence, length, mismatch tolerance, and ambiguous bases. The library structure file also includes a sequence number for each feature sequence. In one specific embodiment, the library structure file is a tab-delimited text file. As shown in Table 1, a library structure file sequentially includes sequence number (ID), feature name (name), position label (pos), start position (start), end position (end), whitelist sequence (seq), length (length), mismatch tolerance (mismatch), and ambiguous bases (Npolicy). When matching the whitelist, Npolicy defines whether a sequencing result of N is considered a mismatch. In one or more embodiments, if the Npolicy item is not specified, the default value is hard (strict matching); otherwise, it is soft (relaxed matching).

[0044] Table 1. Library Structure Files

[0045]

[0046] Specifically, the `pos` field is separated by a colon (`:`). The first part is the type, and the second part is the value. The type part accepts any one of C, D, or R. C indicates that the feature sequence is located in a fixed sequencing file, and the value represents the tag of the sequencing file; D indicates that the position of the feature sequence is the same as the position of the dependent structure, and the value represents the name of the dependent structure; R indicates that the position of the structure is opposite to that of the dependent structure, and the value represents the name of the dependent structure. Only conversion between R1 / R2 and R1 / R2 is supported.

[0047] Specifically, the `start` field is separated by a colon (`:`). The first part is the type, and the second part is the value. The type part accepts any one of C, D, or N. C represents that the starting position of the feature sequence is fixed, and the value represents the position; D represents that the starting position of the feature sequence needs to refer to the dependency structure, and the value is separated by a comma (`,`) and consists of: the name of the dependency structure, the attribute of the dependency structure, and the offset; N represents that the starting position of the feature sequence is uncertain, and the value is fixed as "N".

[0048] Specifically, the `end` field is separated by a colon (`:`). The first part is the type, and the second part is the value. The type part accepts any one of C, D, or N. C represents that the termination position of the feature sequence is fixed, and the value represents the position; D represents that the termination position of the feature sequence needs to refer to the dependency structure, and the value is separated by commas, which are: the name of the dependency structure, the attribute of the dependency structure, and the offset; N represents that the termination position of the feature sequence is uncertain, and the value is fixed as "N".

[0049] Specifically, the seq items are separated by ":", with the first part being the type and the second part being the value. The type part accepts any one of C / F / N / H. C represents a fixed whitelist for the structure, with values ​​separated by "," representing a few fixed sequences. F represents a fixed whitelist for the structure, but stored in a file, with the value representing the filename, which is a text file separated by "," (as shown in Table 2). N represents a structure without a whitelist, with the value being a fixed "N". H represents obtaining the sequence directly at the head of the read operation using a regular expression, with the value being the regular expression.

[0050] Table 2. Whitelist

[0051]

[0052] Specifically, the length field is separated by a colon (:), with the first part being the type and the second part being the value. The type part accepts either C or N. C represents a fixed length structure with the value being the length; N represents a variable length structure with the value being "N".

[0053] In one specific embodiment, detailed steps for preparing the output structure file are provided. The output structure file is used to define the output pattern attributes of the feature sequences in the filtered reads. The output pattern attributes include output position labels, output order, feature name, retrieval method, output type, output sequence start position, and output sequence end position. The output structure file also includes a sequence number for each feature sequence. In one specific embodiment, the output structure file is a tab-delimited text file. As shown in Table 3, an output structure file sequentially includes a sequence number (ID), position label (pos), output order (order), feature name (name), retrieval method (seq), output type (type), output sequence start position (start), and output sequence end position (end). The retrieval method (seq) defines whether the output sequence is a detected sequence (T) or a searched sequence (Q).

[0054] Table 3. Output Structure Files

[0055]

[0056] like Figure 2 The illustrated embodiment provides specific steps for reading the library structure. Each library structure attribute in the library structure file is read sequentially as a filter condition attribute, forming the filter conditions. That is, name, pos, start, end, length, Npolicy, and mismatch are read sequentially, and finally the whitelist sequence (seq) is read. After reading, the library structure attributes of each feature sequence are converted into filter conditions for reads. Specifically, during the reading of the whitelist sequence, each possible sequence is read sequentially as a prefix tree (tre) for subsequent fuzzy sequence matching. The construction method of the prefix tree depends on whether the starting position of the feature sequence is fixed. When the starting position is fixed, a forward prefix tree is constructed; when the starting position is uncertain, a backward prefix tree is constructed. That is, if the type of start is N, the read sequence is a backward prefix tree; otherwise, it is a forward prefix tree.

[0057] .

[0058] like Figure 4The illustrated embodiment provides specific steps for reading a whitelist sequence using a prefix tree. The method for constructing the prefix tree includes: initializing a node; for each sequence in the whitelist, adding each base as a child node to the prefix tree sequentially; if the current base is the last base of the sequence, setting the sequence termination attribute of the current node to true and the sequence character attribute of the node to the full sequence length; otherwise, setting the sequence termination attribute of the previous node to false and setting the current node as a new child node. More specifically, it includes the following steps: First, initialize the prefix tree by creating an empty root node. Determine whether a forward prefix tree or a backward prefix tree needs to be constructed. If it is a backward prefix tree, reverse the sequence to ensure a correct match with the reversed sequence (processedseq) when searching forward from the end of the read. Then, read each sequence in the whitelist sequentially, adding each base (including the fuzzy base N) to the tree starting from its first base. During this process, the current node is initialized as the root node. Then, node creation, traversal, and sequence termination marking are performed according to the following logic: Starting from the first base of the sequence, check if the current node has a child node corresponding to that base. If not, create a new child node representing the current base. If it already exists, set that child node as the current node. Then, move to the next base of the sequence and repeat this process, i.e., check if the new current node has a child node corresponding to the next base, and so on, which is equivalent to extending downwards along a path from the root node. When adding each base, it is necessary to determine whether it is the last base of the sequence. If the current base is the end of the sequence, mark the sequence termination attribute of the current node as TRUE and set the sequence character attribute of the node to the full sequence length. If the current base is not the end of the sequence, mark the sequence termination attribute of the current node as FALSE and update the current node to the newly created child node so that the next base in the sequence can be added. Through the above steps, all sequences in the whitelist are integrated into the prefix tree, and each sequence corresponds to a path from the root node to a certain terminal node in the tree, thus adapting to fast prefix matching. Furthermore, by allowing traversal of different child nodes during the matching process (i.e., allowing base mismatches), fault-tolerant matching is supported, thereby adapting to larger whitelists and higher mismatch capacity, and enabling multi-threaded fast processing of large amounts of data.

[0059] like Figure 3 The illustrated embodiment provides specific steps for reading the output structure. Each output pattern attribute (pos, order, name, seq, type, start, and end) is read sequentially as the output structure, and all output structures are sorted in descending order by order.

[0060] like Figure 5The illustrated embodiment provides specific steps for read filtering. These steps can be summarized as follows: sequentially read all ends of the sequencing results for each read (including R1, R2, or any custom tags); sequentially check each filtering condition; determine the query sequence region, perform sequence matching and querying; record the query results and determine the filtering criteria.

[0061] The method for determining the query sequence region includes the following steps: Based on the definition of the `pos` item in the filter conditions, determine which sequencing file (e.g., R1 or R2) the sequence template will be obtained from (if the type of `pos` is R (dependency), then look for the dependent filter structure: if the dependent `pos` is R1, then actually search the R2 file; and vice versa). If the type of the `start` (start position) or `end` (end position) of the filter conditions is D (dependency), then look for the filtering result of the sequence feature it depends on (e.g., the end position of the dependent UMI feature), and add a defined offset to this value to calculate the expected position of the current feature sequence. Based on the start position, end position, and feature length, calculate the expected start position (`startq`) and expected end position (`endq`) of the query sequence to obtain the expected query interval. The calculation formula is as follows:

[0062]

[0063] Then, the template of the query sequence is obtained according to pos, and the length L of the template sequence is calculated. The expected query interval is compared with the actual length of the template, and the actual start position (startqf) and actual end position (endqf) of the query sequence are calculated. The calculation formula is as follows:

[0064]

[0065] If startqf equals -1, the query is considered unsuccessful and the feature matching is deemed to have failed. If startqf is not equal to -1, the query sequence seqq is obtained from the template sequence according to startqf and endqf, and sequence matching and query are performed accordingly.

[0066] The sequence matching and query method includes the following steps. If the filter condition is not configured with a whitelist (i.e., the prefix tree tre is empty), the actual search interval startqf and endqf are directly used as the matching result. If a whitelist is configured (tre is not empty), the range of the starting position s is set to 1 to L, and the prefix tree is used to search the query sequence in seqq sequentially from left to right: the position of the temporary query sequence seqqt is from s to endqf, the current node is set as the root node of the prefix tree, the current base is set as the first base of seqqt, the current mismatch base count mis is set to 0, and the prefix tree traversal of seqqt begins. The traversal logic includes: if the sequence termination attribute of the current node is TRUE, the match is successful, and seqt is the sequence character attribute of the current node; if the mismatch count mis has exceeded the maximum fault tolerance max_mis, the search is terminated; if the sequence termination attribute of the current node is FALSE, the search is attempted to continue to match the next base. The methods for continuing the search include: querying the child nodes of the current node; if the current base exists in a child node, then set the current node as a child node that meets the conditions, keep mis unchanged, and continue the search; if the child nodes of the current node contain N, and the filtering conditions tolerate N, then set the current node as child node N, keep mis unchanged, and continue the search; querying all child nodes of the current node, searching all the queried child nodes, setting the current node as a child node, and incrementing mis by 1; once a match is found at a certain position s, stop the search and do not continue to query other positions s.

[0067] After obtaining the query results, calculate the starting position startm and the ending position endm of the match, using the following formula:

[0068]

[0069] Record at least one of the query sequence seqt, the matching sequence seqqt, the start position of the match, and the end position of the match as the filtering result of the corresponding query condition; if the read successfully passes the matching of all filtering conditions, it is determined that the read conforms to the document structure and is marked as the filtered read; record the matching results of all features of the read for subsequent output or deduplication steps.

[0070] like Figure 6The illustrated embodiment provides a read deduplication step. This step can be summarized as follows: grouping the filtered reads according to their corresponding barcode and UMI; for multiple reads within the same group, calculating a consensus sequence based on their target sequence; when the support of the consensus sequence is greater than or equal to a preset threshold, outputting the consensus sequence as the representative sequence of that group. More specifically, it includes the following steps: For all filtered reads, using their corresponding barcode+UMI string as a unique identifier (ID), and using this ID as the key. For each read to be processed, perform the following operations: Check if the ID exists in the database; if the ID does not exist, create the following four subkeys for this ID: best_seq (to store the current best sequence, i.e., the consensus sequence), count (representing the total number of reads in this ID group), best_count (representing the number of supporters for the current best sequence), and check_list (a list). The subkey best_seq records the current read, the subkey count is set to 1, the subkey best_count is set to 1, and the current read is added to check_list; if the ID already exists, increment the count value of that ID by 1, and add the current read to the check_list of that ID. Compare the current read with the best_seq of this ID to see if they are sufficiently similar (by default, less than 5% base inconsistency is considered sufficient similarity between the two sequences). If they are sufficiently similar, increment the best_count of the ID by 1. After storing all reads in the database, the database is traversed, and the final output is determined based on the consistency of the sequences within the group: if the best_count / count of an ID is >= 0.9, then the sequence of that ID is considered to be best_seq and output as the unique representative sequence of the group; if the best_count / count of an ID is < 0.9, then check_list is traversed, and the possible sequences for each ID are output.

[0071] like Figure 7The illustrated embodiment provides a reads output step. For filtered reads, the specified sequence feature content is output to the corresponding file according to the format defined in the output structure file. Specifically, this includes the following steps: Output each output structure in order of order. Find the corresponding filter item in the filtering results according to type. If seq is T, the template sequence is target seq (matching sequence); if seq is Q, the template sequence is query seq (query sequence). If name is whole_seq, the output start position starto is 1, and the output end position endo is the length of the template sequence; if name is not whole, the output start position starto is the start attribute of the structure, and the output end position endo is the end attribute of the structure. Calculate starto and endo; if starto is N, change starto to 1; if endo is N, change endo to the length of the template sequence. Write the sequence fragment from starto to endo of the template sequence to the corresponding output file specified by the pos field. Iterate through all filtered reads or deduplicated reads (ID), repeating the above steps until all data has been processed and output according to the output structure file.

[0072] In one specific embodiment, a device for detecting feature sequences of single-cell sequencing data adapted to multiple library structures is also provided, including a processor and a memory; the processor and the memory are connected via a communication bus; wherein, the processor is used to call and execute a program stored in the memory; the memory is used to store the program, which is at least used to execute the single-cell sequencing data feature sequence detection method adapted to multiple library structures.

[0073] In one specific embodiment, bacdrop ​​article data (Ma P, AmemiyaHM, He LL, Gandhi SJ, Nicol R, Bhattacharyya RP, Smillie CS, Hung DT. Bacterial droplet-based single-cell RNA-seq reveals antibiotic-associated heterogeneous cellular states. Cell. 2023 Feb 16;186(4):877-891.e14. doi:10.1016 / j.cell.2023.01.002. Epub 2023 Jan 27. PMID: 36708705; PMCID:PMC10014032.) was constructed. The library structure file shown in Table 4 and the output structure file shown in Table 5 were constructed. The SRR15174673 sample was analyzed, and the filtering results shown in Table 6 were obtained.

[0074] Table 4. Bacdrop ​​Library Structure Files

[0075]

[0076] Table 5. Bacdrop ​​Output Structure Files

[0077]

[0078] Table 6. Bacdrop ​​Filtering Results

[0079]

[0080] Table 6 describes the overall results of the filtering analysis, where: processed_reads is the number of input reads, passed_reads is the number of output reads, barcode1 is the number of reads that failed to be filtered due to the barcode1 structure, and barcode2 is the number of reads that failed to be filtered due to the barcode2 structure. The overall pass rate for filtered reads is above 97%, consistent with the results described in the article.

[0081] In one specific embodiment, smrandom-seq article data (Xu Z, Wang Y, Sheng K, Rosenthal R, Liu N, Hua X, Zhang T, Chen J, Song M, Lv Y, Zhang S, Huang Y, Wang Z, Cao T, Shen Y, Jiang Y, Yu Y, Chen Y, Guo G, Yin P, WeitzDA, Wang Y. Droplet-based high-throughput single microbe RNA sequencing by smRandom-seq. Nat Commun. 2023 Aug 23;14(1):5130. doi: 10.1038 / s41467-023-40137-9. PMID: 37612289; PMCID: PMC10447461.) was constructed. The library structure file shown in Table 7 and the output structure file shown in Table 8 were constructed. The CRR794623 sample was analyzed, and the filtering results shown in Table 9 were obtained.

[0082] Table 7. smrandom-seq library structure file

[0083]

[0084] Table 8. smrandom-seq output structure file

[0085]

[0086] Table 9. smrandom-seq filtering results

[0087]

[0088] Table 9 describes the overall results of the filtering analysis, where: processed_reads is the number of input reads, passed_reads is the number of output reads, barcode1 is the number of reads that failed to be filtered due to the barcode1 structure, and barcode2 is the number of reads that failed to be filtered due to the barcode2 structure. The overall pass rate for filtered reads is over 90%, consistent with the results described in the article.

[0089] In one specific embodiment, a genome data is processed according to the above method to obtain the filtering results shown in Table 10. After deduplication of each barcode, the deduplication results shown in Table 11 are obtained.

[0090] Table 10. Results of Filtering Genomic Data

[0091]

[0092] Table 10 describes the overall results of the filtering analysis, where: processed_reads is the number of input reads, passed_reads is the number of output reads, bc1 is the number of reads that failed to be filtered due to the barcode1 structure, bc2 is the number of reads that failed to be filtered due to the barcode2 structure, and w1 is a fixed sequence between bc1 and bc2.

[0093] Table 11. Deduplication Results of Genomic Data

[0094]

[0095] Table 11 describes the overall results after deduplication, where: SAG_id is the barcode, read_count is the number of reads, UMI_count is the number of UMIs, and deduplicate_count is the number of reads after deduplication.

[0096] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope of the present invention.

Claims

1. A method for detecting characteristic sequences in single-cell sequencing data that is adaptable to various library structures, characterized in that, Includes the following steps: Prepare a library structure file to define the library structure attributes of the feature sequences in the reads to be filtered; the library structure attributes include whitelisted sequences and mismatch tolerance. Prepare the output structure file to define the output pattern attributes of the feature sequences in the filtered reads; Each library structure attribute in the library structure file is read as a filter condition attribute to form a filter condition; wherein, the whitelist sequence in the library structure attribute is read through a constructed prefix tree for fuzzy sequence matching; the construction method of the prefix tree is determined according to whether the starting position of the feature sequence is fixed. When the starting position is fixed, a forward prefix tree is constructed, and when the starting position is not fixed, a reverse prefix tree is constructed. Read each output mode attribute from the output structure file and use it as the output structure; Acquire single-cell sequencing data and read each read sequentially. For each read, query and match each feature sequence according to the filtering conditions, obtain the query results, and record the query sequence, matching sequence, matching start position, and matching end position as the filtering results of the corresponding query conditions. If the read passes the filtering of all feature sequences, it is considered to conform to the library structure and is regarded as a filtered read. Among them, for feature sequences configured with whitelist sequences, the prefix tree is used for fault-tolerant matching in the target region of the read, and the number of mismatches allowed does not exceed the mismatch tolerance. For the filtered reads, the specified sequence feature content is output to the corresponding file according to the format defined in the output structure file.

2. The method as described in claim 1, characterized in that, The library structure attributes also include at least one or more of the following: feature name, position label, start position, end position, and length; The output pattern attributes include at least one or more of the following: output position label, output order, feature name, retrieval method, output type, output sequence start position, and output sequence end position.

3. The method as described in claim 1, characterized in that, The method for constructing the prefix tree includes: initializing a node; for each sequence in the whitelist sequence, adding each base of the sequence as a child node in sequence; if the current base is the last base of the sequence, then setting the sequence termination attribute of the current node to true and setting the sequence character attribute of the node to the full length of the sequence, otherwise setting the sequence termination attribute of the current node to false and setting the current node as a new child node.

4. The method as described in claim 1 or 2, characterized in that, The method for querying and matching each feature sequence in sequence according to the filtering conditions includes: calculating the actual start position and actual end position of the feature sequence in the corresponding reads based on the position label, start position, end position, and length defined in the library structure file, and obtaining the query sequence accordingly; if the feature sequence is configured with a whitelist sequence, then the prefix tree is used on the query sequence for fault-tolerant matching, and the number of mismatches allowed does not exceed the mismatch tolerance, and the query result is obtained; if the feature sequence is not configured with a whitelist sequence, then the actual start position and actual end position of the target region are directly used as the query result.

5. The method as described in claim 1, characterized in that, The method for fault-tolerant matching using the prefix tree includes: if there is no child node in the current node whose child node matches the base to be matched, then the search is performed by querying other child nodes of the current node and incrementing the mismatch count until a match is successful or the accumulated mismatch count exceeds the mismatch tolerance.

6. The method as described in claim 1, characterized in that, The method further includes deduplicating the filtered reads: grouping the filtered reads according to their corresponding barcode and UMI; calculating a consensus sequence based on the target sequence for multiple reads in the same group; and outputting the consensus sequence as the representative sequence of the group when the support of the consensus sequence is greater than or equal to a preset threshold.

7. The method as described in claim 6, characterized in that, The support level is the proportion of the number of reads supporting the consensus sequence to the total number of reads in the group, and the preset threshold is 0.

9.

8. The method as described in claim 1 or 2, characterized in that, Includes at least one of the following (a) to (b), (a) The start or end position of the feature sequence is defined by any of the following methods (1) to (3): (1) Based on the definition of a fixed position, (2) Based on its positional attributes relative to the dependent structure and combined with the offset, a relative definition is made. (3) Defined as an uncertain position; (ii) The whitelist sequence is provided in any of the following ways (a) to (d): (a) Directly specify one or more fixed sequences, (b) Provided by referencing an external file containing a list of possible valid sequences. (c) No specific whitelist sequence, (d) Obtain the sequence from the header information of read using regular expressions; (iii) The library structure attributes also include fuzzy bases; (iv) The method described herein is applicable to single-cell transcriptome sequencing data and / or single-cell genome sequencing data.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1 to 8.

Images