High-throughput sequencing library host rapid identification method, system, device and storage medium

Abstract

The application discloses a high-throughput sequencing library host rapid identification method, system, device and storage medium, belonging to the field of biological information sequencing data analysis. The nucleotide homologous alignment of contig sequences is obtained in the 6th standard output format, and three types of mitochondrial marker sequences, COX1, Cytb and D-loop, are matched. The three types of filtering, alignment length, similarity and coverage, are performed. The species classification pedigree is recalled from the classification database, and invalid sequences are removed. The number of marker gene species, the number of matched sequences and the average similarity are counted according to the species, the hit marker gene species are sorted from high to low according to the number, and the species with a hit marker gene species number greater than or equal to 2 are determined as high-confidence hosts. The application completes the construction of a pre-database, generates a classification database and a marker keyword library in a resident memory, supports custom traceable nodes and adjustable filtering thresholds, and realizes full-automatic offline batch operation, so that the efficiency is greatly improved, and the application is suitable for wild animal monitoring, pathogen tracing and other host rapid identification.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics sequencing data analysis technology, and in particular to a method, system, device, and storage medium for rapid host identification of high-throughput sequencing libraries. Background Technology

[0002] Existing host identification and analysis methods suffer from several practical problems. First, manual analysis is inefficient for large-scale sample volumes. With the decreasing cost of high-throughput sequencing, research projects often involve hundreds or even thousands of samples, with each library's BLAST alignment containing tens to hundreds of thousands of sequence hits. Manual screening and processing require significant manpower and time, severely hindering project progress. Second, manual analysis results suffer from poor consistency. Different researchers have varying interpretations of screening criteria, making cross-sectional comparisons between different batches and from different personnel difficult, hindering data integration. Third, existing online identification tools have limited applicability, mostly focusing on single-sequence species identification. They cannot directly interface with batch alignment results from metagenomic sequencing and do not support custom screening parameters and taxonomic ranks, failing to meet the diverse analytical needs of pathogen surveillance scenarios. Fourth, taxonomic information retrieval is inefficient. Currently, obtaining species classification information largely relies on online calls to the NCBI Taxonomy database interface, which is significantly affected by network speed during batch analysis and is prone to interruption due to interface limitations. Fifth, host identification has a high professional threshold, requiring researchers to master bioinformatics operations, taxonomy knowledge, and sequence alignment skills, which is not conducive to its promotion and application in grassroots monitoring networks. Summary of the Invention

[0003] To overcome the above shortcomings, the present invention aims to provide a method, system, device and storage medium for rapid identification of hosts in high-throughput sequencing libraries, so as to solve the defects of existing sequencing library host identification processes, such as cumbersome process, high dependence on manual operation, high professional threshold, low efficiency of batch analysis, fixed identification order, inability to screen target species groups in a targeted manner and inconsistent identification result standards.

[0004] To achieve the above objectives, the first aspect of the present invention provides a method for rapid host identification of high-throughput sequencing libraries, including an automatic pre-database construction step and an automated analysis step consisting of S1-S6; the automatic pre-database construction step is completed before the execution of S2, and the program generates a classification database containing four types of memory-mapped structures and resident in memory by filtering animal classification nodes with taxonomic domain identifier = 5, and simultaneously constructs a tag keyword library that includes all annotation variants of COX1, Cytb, and D-loop; The automated analysis steps include: S1, Obtain the alignment file with the field structure of the Class 6 standard output format, which is the output of the contiguous group sequence after nucleotide homology alignment; S2, read the comparison file, match the keywords corresponding to the mitochondrial DNA barcode markers of three types of animals: COX1, Cytb, and D-loop, and screen to obtain target gene matching sequences containing only the target marker; The alignment file is read, and for the target sequence title field of each alignment record, case-insensitive substring matching is performed. Keywords corresponding to the three types of animal mitochondrial DNA barcode markers—COX1, Cytb, and D-loop—are matched sequentially. A sequence matching any one of these marker keywords in any position within the title is considered a target gene match sequence, is retained, and its corresponding gene type is marked. Records that do not match any marker keyword are directly discarded and do not proceed to the next step. Finally, a set of target gene match sequences containing only the target markers is obtained. S3, perform a three-layer linkage filtering on the target gene matching sequence in sequence, which includes alignment length, sequence similarity and alignment coverage, to obtain a valid candidate sequence that simultaneously meets the three-layer threshold constraints; The target gene matching sequence is subjected to a three-layer linkage filtering process: alignment length, sequence similarity, and alignment coverage. Each layer of filtering is independent, and sequences that do not meet the threshold of the current layer are directly eliminated and do not proceed to the next layer of filtering. Only sequences that pass all three layers of threshold constraints are considered valid candidate sequences. The filtering rules for each layer are as follows: Alignment length filtering: The alignment length value in the alignment record is directly extracted and compared with the preset minimum alignment length threshold. If the value is greater than or equal to the threshold, it passes; otherwise, it is eliminated. Sequence similarity filtering: The sequence similarity percentage value in the alignment record is extracted and compared with the preset minimum similarity threshold. If the value is greater than or equal to the threshold, it passes; otherwise, it is eliminated. Alignment coverage filtering: The query-side coverage of each sequence is calculated according to the formula: alignment coverage = (alignment length ÷ total length of contiguous sequences) × 100%, and compared with the preset minimum coverage threshold. If the value is greater than or equal to the threshold, it passes; otherwise, it is eliminated. S4. Extract the species classification number of the valid candidate sequence, verify the classification number to which the species classification number belongs to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back upwards step by step. Remove non-animal sequences and sequences with incomplete classification systems to obtain a complete classification sequence. First, the sequence is checked against a pre-stored set of animal classification numbers using an algorithm with time complexity to determine if it belongs to an animal group. If it does not, it is directly identified as a non-animal sequence and removed. For sequences belonging to an animal group, a hierarchical backtracking process is initiated. The specific process is as follows: the species classification number of the current sequence is set as the starting node, and the lineage string is initialized to empty. Through the parent node mapping relationship in the classification database, the classification number of the parent node above the current classification number is queried, and the taxonomic rank and species name of the node are obtained simultaneously. These are then concatenated to the beginning of the lineage string in the format "species name (taxonomic rank);". The current classification number is updated to the parent node classification number. This step is repeated. When the current classification number is equal to the preset root classification number, the backtracking stops. The root node information is added to the head of the lineage to generate a complete hierarchical classification lineage. If, during the backtracking process, a parent node does not exist or the root classification number cannot be reached, it is identified as a sequence with an incomplete lineage and is removed. Finally, a complete parsed sequence of all lineages is obtained. S5, group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity in each group; species with ≥2 marker gene types are identified as high-confidence hosts; when ranking candidate species, they are first arranged from highest to lowest according to the number of marker gene types hit, with those having more types ranked higher; when different species have the same number of marker gene types hit, these species are ranked from highest to lowest according to their comprehensive score; when all candidate species hit only one marker gene, all candidate species are directly ranked from highest to lowest according to their comprehensive score; comprehensive score = number of matching sequences × average similarity; The complete sequence analysis of the phylogenetic tree was grouped by species using species classification numbers as unique identifiers. For each species group, three core indicators were calculated: the number of marker gene types hit (counting all marker gene types in all sequences of that species after deduplication, ranging from 1 to 3), the number of matching sequences (the total number of all valid sequences in that species), and the average similarity (the sum of similarity values ​​of all valid sequences in that species ÷ the number of matching sequences, rounded to three decimal places). Species with ≥2 hit marker gene types were classified as high-confidence hosts. Candidate species were ranked according to the following two-tiered rules: the first priority was the number of hit marker gene types, with higher numbers ranking higher; the second priority was the comprehensive score. When different species had the same number of hit marker gene types, they were ranked from highest to lowest based on the comprehensive score, which was calculated as: comprehensive score = number of matching sequences × average similarity. A higher comprehensive score indicated higher host confidence and a higher ranking. When all candidate species hit only one marker gene, all candidate species were directly ranked from highest to lowest based on their comprehensive score.

[0005] S6 outputs multiple sets of candidate host information and generates a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

[0006] The second aspect of the present invention provides a rapid host identification system for high-throughput sequencing libraries, including a database construction module, a data input module, a gene screening module, a threshold filtering module, a classification and tracing module, a host determination module, and a result output and statistics module; The database construction module is used to perform the automatic construction of the pre-database, filter animal classification nodes with classification domain identifier = 5, generate a classification database containing four types of memory mapping and resident in memory, and simultaneously build a tag keyword library that includes all annotation variants of COX1, Cytb, and D-loop. The data input module is used to execute the acquisition of alignment files with the nucleotide homology alignment output of contig sequences and the field structure of the sixth type of standard output format. The gene screening module is used to read the alignment file, match the keywords corresponding to the mitochondrial DNA barcode markers of three types of animals (COX1, Cytb, and D-loop), and screen to obtain target gene matching sequences containing only the target marker. The threshold filtering module is used to perform a three-layer linkage filtering on the target gene matching sequence in sequence, which considers alignment length, sequence similarity, and alignment coverage, to obtain an effective candidate sequence that simultaneously meets the three-layer threshold constraints. The classification and tracing module is used to extract the species classification number of the valid candidate sequence, verify the species classification number to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back upwards step by step, eliminating non-animal sequences and sequences with incomplete system, and obtaining a complete system analysis sequence. The host determination module is used to group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity in each group; species with ≥2 marker gene types are identified as high-confidence hosts; when ranking candidate species, they are first arranged from highest to lowest by the number of marker gene types, with those having more types ranked higher; when different species have the same number of marker gene types, these species are ranked from highest to lowest by their comprehensive score; when all candidate species have only one marker gene, all candidate species are directly ranked from highest to lowest by their comprehensive score; comprehensive score = number of matching sequences × average similarity; The result output and statistics module is used to output multiple sets of candidate host information and generate a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

[0007] A third aspect of the present invention provides an electronic device for rapid identification of hosts in high-throughput sequencing libraries, wherein the system stores a computer program, and the program executes the aforementioned identification method when it runs.

[0008] A fourth aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the program, when executed by a processor, implements the aforementioned identification method.

[0009] The beneficial effects of this application, based on the above technical solutions, are as follows: This invention addresses the shortcomings of existing host identification processes, such as high reliance on manual labor, low batch processing efficiency, and susceptibility to human error. It directly reads nucleotide homology alignment results from contiguous sequences, automating the entire process of marker screening, sequence filtering, species lineage tracing, and host determination without human intervention, filling a gap in automated analysis tools. This invention simultaneously identifies three types of mitochondrial barcode markers: COX1, Cytb, and D-loop. Combined with multi-layer adjustable filtering thresholds, it achieves multi-gene cross-validation, effectively avoiding sequence omissions and species misidentification, and reducing false positives. Furthermore, this invention supports customizable classification tracing termination nodes via root classification numbers, enabling precise species identification or rapid screening of higher-order groups as needed, selectively removing irrelevant sequences. Coupled with a lightweight, resident memory classification database that only loads animal classification nodes, it reduces memory usage, improves lineage tracing speed, and is suitable for large-scale sample analysis. The overall tool has a low operational threshold, requiring no complex bioinformatics background, and uses a unified judgment standard throughout, eliminating human subjectivity and ensuring standardized and reproducible identification results. Traditional manual analysis of a single sequencing library host takes 1-2 hours. The fully automated process of this invention can reduce single sample identification to a few minutes, greatly improving the efficiency of analyzing and processing large batches of monitoring samples. Attached Figure Description

[0010] Figure 1 This is a flowchart of a method for rapid host identification of high-throughput sequencing libraries provided in an embodiment of the present invention; Figure 2 This is a flowchart illustrating the preloading process for the classification database and tagging keywords in an embodiment of the present invention. Figure 3 This is a flowchart of mitochondrial marker gene screening and three-level threshold filtering provided in the embodiments of the present invention; Figure 4 This is a flowchart of the automatic tracing process for classification lineages provided in an embodiment of the present invention; Figure 5 This is a flowchart of the host species determination and statistical ranking output provided in the embodiments of the present invention; Figure 6 This is a block diagram of a high-throughput sequencing library host rapid identification system provided in an embodiment of the present invention. Detailed Implementation

[0011] To fully illustrate the implementation process of this invention, the technical solution of this invention will be described below with reference to the corresponding implementation process in the accompanying drawings. It should be understood that the following implementation process is only used to explain this invention and is not intended to limit the scope of protection of this invention; simple modifications and parameter adjustments made based on the overall concept of the method of this invention all fall within the scope of protection claimed by this invention.

[0012] Some technical terms in this invention are defined as follows: contig: A contig sequence is a gapless, continuous nucleic acid sequence obtained by splicing together short reads from overlapping sequencing. It is a core sequence unit used for species alignment in metagenomic / metavinomic sequencing.

[0013] de novo assembly: This refers to a bioinformatics analysis process that assembles long contiguous sequences from raw sequencing reads without a reference genome, relying solely on the overlap relationships between these reads.

[0014] BLASTn: A nucleic acid sequence homology alignment tool used to align a nucleic acid sequence to be tested with a nucleic acid reference database and output homology information such as sequence similarity and matching position. In this invention, it is used for species homology retrieval between contiguous group sequences and nucleotide reference databases.

[0015] E-value: The expected value of the alignment, a parameter of the BLAST alignment, represents the probability of randomly matching equally similar sequences; the smaller the value, the higher the credibility of the sequence homology matching. This indicator is retained in the standardized result table of this invention to assist in screening credible matches.

[0016] O(1): Algorithm time complexity indicator, representing that dictionary / set key-value lookup has a fixed time consumption and does not increase the retrieval time with the increase of the number of database entries. The animal species classification number set of this invention adopts this query method to achieve millisecond-level classification verification.

[0017] Command line parameters: Adjustable control parameters passed in from outside the program during runtime; --min-length is used to customize the minimum threshold for alignment length, and --min-ident is used to customize the minimum threshold for sequence similarity. Users can modify the filtering criteria as needed to adapt to different samples.

[0018] Comparison Coverage: This invention uniformly adopts the overlap group sequence query sequence coverage standard, and the calculation method is as follows: The comparison length is calculated as: (length of alignment ÷ length of complete contiguous sequence) × 100%. There is another method in the industry for calculating the coverage of the Subject reference sequence, but this invention does not use this standard throughout the process, thereby avoiding false positive identification of species caused by matching short conserved reference fragments.

[0019] The complete operating process of this invention will be described in detail below with reference to actual test samples.

[0020] 1. Test Sample The test sample was a mixed droppings sample from wild birds. After high-throughput sequencing and de novo assembly, a total of 9209 contig sequences were obtained, with sequence lengths ranging from 201 to 19783 base pairs.

[0021] 2. Preparation of taxonomic database Open the official FTP download page for the NCBI Taxonomy database and download the taxdmp.tar.gz categorical data compressed package. After decompression, you will get nodes.dmp (node ​​data file) and names.dmp (name data file). In this example, the target group is specified as Aves (root taxonomic number = 8782). The tool automatically parses the nodes.dmp file and preloads only the taxonomic nodes under Aves to build an in-memory dictionary, including taxid - parent node mapping dictionary, taxid - taxonomic rank dictionary, and taxid - species name dictionary.

[0022] This tool supports custom traceability termination nodes, allowing the root classification number to be modified to Anseriformes 8838 to rebuild the classification database; it loads only about 1800 classification nodes of Anseriformes, reducing memory usage to 1.2MB and database loading time to less than 0.1 seconds; it retains only Anseriformes species for specific host monitoring, enabling targeted group screening.

[0023] 3. Preparation of input files for nucleotide homology alignment The NCBI BLAST+ tool was used to align the assembled contiguous sequences with the NCBI nucleotide reference database for nucleotide homology. The output format was set to the Class 6 standard output format. The final input file field definitions are as follows: qseqid qlen sseqid stitle pident length evalue superkingdom staxidsLineage; Input sample fragment: TRINITY_xxx; 1650; srRNA;96.424;1650;0;Eukaryota;8849;Aves(class);Neognathae(infraclass);Galloanserae(superorder);Anseriformes(order);Anatidae(family);Anserinae(subfamily);Anser(genus);Anser caerulescens(species).

[0024] 4. Dual-layer pre-filter The tool reads the BLAST result file line by line, first performing field integrity verification and removing invalid records with empty species classification numbers or missing key fields; then it performs keyword pre-screening of mitochondrial marker genes, matching keywords in the stitle field, retaining any matching mitochondrial marker keywords, otherwise automatically filtering out the matching sequences. A total of 73 matching mitochondrial marker gene sequences were selected, including 17 sequences from the swan goose (Anser cygnoides), 12 sequences from the black swan (Cygnusatratus), and 5 sequences from the bar-headed goose (Anser indicus), etc. The remaining 31,018 non-mitochondrial gene sequences were automatically filtered out.

[0025] 5. Three-level threshold filtering This embodiment uses a custom threshold, with the similarity threshold lowered to 90%, allowing more closely related species to pass: Alignment length threshold: ≥500 base pairs; Sequence similarity threshold: ≥90%; Comparison coverage threshold: ≥60%; The 73 pre-screened records were filtered step by step: finally, 6 valid candidate comparison records were retained. The list of results after filtering is as follows: 2 records of snow goose (Anser caerulescens), 2 records of swan goose (Anser cygnoides), 1 record of black swan (Cygnusatratus), and 1 record of brown hard-tailed duck (Oxyura jamaicensis).

[0026] 6. Classification and lineage tracing For all 6 valid candidate sequences, the classification lineage was traced back one by one. Starting from the species classification number corresponding to each sequence, the parent node was searched upwards level by level. At each level of tracing, the name of the rank was added to the lineage string. The tracing stopped when the root node of the class Aves (species classification number = 8782) was reached, and the complete classification lineage information of each sequence was obtained.

[0027] 7. Host species determination and result output Valid sequences are grouped by species name. The number of matched marker gene types for each species is counted. When the number of matched marker gene types is different, the higher the number of types, the higher the ranking. Species with ≥2 matched marker gene types are directly identified as high-confidence hosts and ranked higher than species with only one matched marker gene. For species with the same number of matched marker gene types, a comprehensive score is calculated using the formula (Comprehensive Score = Number of Matching Sequences × Average Similarity). The higher the score, the higher the ranking. The top 5 candidate hosts are output by default for easy manual verification. The comprehensive score model is then used to complete the competitive ranking of multiple candidate species. Anser cygnoides: There are 2 valid sequences, which match COX1 and Cytb tags respectively, with an average similarity of 99.280%. It meets the condition of hitting two or more tags and is directly determined as a high-confidence dominant host in this library. Xueyan (Anser caerulescens): Matched COX1 marker, 2 valid sequences, average similarity 96.244%, overall score 192.488; Black Swan (Cygnus atratus): Matches Cytb marker, 1 valid sequence, similarity 97.537%, overall score 97.537; Brown hard-tailed duck (Oxyura jamaicensis): Matched D-loop marker, 1 valid sequence, similarity 93.443%, overall score 93.443; The remaining single-marker species, ranked from highest to lowest based on their overall scores, are: Snow Goose, Black Swan, and Brown Hardtailed Duck.

[0028] The species with the highest overall score was the Swan Goose (Anser cygnoides), which was determined to be the dominant host species for this mixed fecal sequencing library.

[0029] The tool generates standardized results, including 10 fields: query sequence ID, target sequence accession number, marker gene type, similarity, alignment length, coverage, E-value, species classification number, species name, and complete taxonomic lineage. It also outputs a host determination statistical summary, including minimum alignment length, minimum similarity, minimum coverage, root node taxid, number of valid candidate species, most likely host, number of marker genes, number of hit sequences, and average similarity, thus completing the entire analysis process.

[0030] Standardized evaluation form templates are shown in Table 1 below: Table 1: Standardized Evaluation Form Row Examples

[0031] This invention's entire analysis process is fully automated, requiring no manual screening, statistics, or sequence verification. While traditional manual processing of mixed samples of the same size takes 1-2 hours, this tool completes the entire analysis in just 0.9 seconds. This embodiment uses a lightweight bird classification library, reducing memory usage by 93.5% compared to the complete NCBI Taxonomy database. Utilizing the COX1+Cytb dual-label cross-validation mechanism, it effectively distinguishes between the Swan Goose and its close relative, the Snow Goose, avoiding the false positive problem common with single-label tools. Furthermore, the tool supports customizable filtering thresholds and customizable classification tracing levels, making it suitable for both general wildlife surveillance and specialized sample identification for single orders or classes, adapting to various pathogen monitoring scenarios.

[0032] Figure 1 This is a flowchart of a high-throughput sequencing library host rapid identification method provided in an embodiment of the present invention, as shown below. Figure 1 As shown, the identification method provided in this embodiment of the invention includes a pre-construction database automatic construction step and an automated analysis step consisting of S1 to S6. The overall process is divided into three stages: database preloading, automated sequence analysis, and host determination output. The pre-construction database automatic construction step must be completed before executing step S2. The specific execution flow is as follows: The pre-construction database automatically builds the following steps: The program selects animal classification nodes with a classification domain identifier of 5, builds a classification database with four types of memory mapping structures that reside in memory, and simultaneously builds a tag keyword library that includes all annotated variants of COX1, Cytb, and D-loop. After the database and keyword library are built, they are stored in memory for subsequent analysis steps.

[0033] The automated analysis steps include: S1, Obtain the alignment file with the field structure of the Class 6 standard output format, which is the output of the contiguous group sequence after nucleotide homology alignment; Specifically, after obtaining the contiguous group sequence from the de novo assembly of high-throughput sequencing data, the contiguous group sequence is compared with the NCBI nucleotide reference database for nucleotide homology. During the comparison, the standard output format file of category 6 is generated by using the command parameter -outfmt “6 qseqidsseqid pident length mismatch gapopen qstart qend sstart send evalue bitscoreqlen staxid stitle”. This file always contains the query sequence ID, target accession number, sequence similarity, alignment length, alignment expectation value, query sequence length, species classification number, and sequence title key fields, and serves as the unified input file for this method.

[0034] This program has a built-in error-tolerant filtering rule for the sixth type of standard output format files: The first layer verifies the integrity of the file fields. If the species classification number field in a row is empty or the comparison length / similarity value is missing, the comparison record is discarded directly. The second layer pre-screens mitochondrial markers. If the stitle field does not match any mitochondrial marker keywords, the row is filtered directly and does not enter the subsequent three-layer threshold filtering process.

[0035] The sixth type of standard output format field structure refers to the set of fields and the order of field arrangement defined by the BLAST tool - outfmt 6. It includes both the original alignment file directly exported by the BLAST software and the alignment file with the same field structure output by other alignment tools.

[0036] S2, read the comparison file, match the keywords corresponding to the mitochondrial DNA barcode markers of three types of animals: COX1, Cytb, and D-loop, and screen to obtain target gene matching sequences containing only the target marker; COX1, Cytb, and D-loop are universal mitochondrial DNA barcodes for animal identification, and their advantages are complementary: COX1 has broad annotation coverage and is applicable to the vast majority of species; Cytb is more effective in distinguishing closely related vertebrate species; and D-loop has high variability, enabling fine-grained species differentiation. Existing tools mostly recognize only a single marker, which is prone to missing sequences due to incomplete annotation or sample degradation, and single-gene identification is also prone to false positives in closely related species. This invention integrates all annotated variants of the three markers and significantly improves the accuracy of host identification through multi-marker cross-validation.

[0037] Specifically, the system automatically selects sequences that match three types of mitochondrial marker genes from a massive BLAST alignment result, filtering out other irrelevant sequences; the marker keyword library includes various annotation variants such as the full name and abbreviation of each gene to avoid missing sequences.

[0038] S3, perform a three-layer linkage filtering on the target gene matching sequence in sequence, which includes alignment length, sequence similarity and alignment coverage, to obtain a valid candidate sequence that simultaneously meets the three-layer threshold constraints; Specifically, the default parameters for the three-layer filtering are alignment length ≥ 500 base pairs, sequence similarity ≥ 97%, and alignment coverage ≥ 60%. The three threshold ranges can be adjusted by the user via command line within the ranges of 300~1000 base pairs, 95%~99%, and 50%~90%.

[0039] The alignment coverage described in this invention specifically refers to the coverage of contiguous group sequence queries. The denominator is the complete sequence length of the contiguous group sequence to be analyzed. Coverage is not calculated based on the sequence length of the NCBI reference mitochondrial gene (subject), which can avoid false positives caused by short conserved fragments.

[0040] S4. Extract the species classification number of the valid candidate sequence, verify the classification number to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back upwards step by step. Remove non-animal sequences and sequences with incomplete classification systems to obtain a complete classification sequence. Specifically, relying on the NCBI classification database, which stores only animal classification nodes, is pre-loaded and resident in memory, to complete the tracing of the complete classification lineage. The root species classification number can be configured as the tracing termination node to filter sequences without valid animal classification information or with broken lineages.

[0041] S5, group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity in each group; species with ≥2 marker gene types are identified as high-confidence hosts; when ranking candidate species, they are first arranged from highest to lowest according to the number of marker gene types hit, with those having more types ranked higher; when different species have the same number of marker gene types hit, these species are ranked from highest to lowest according to their comprehensive score; when all candidate species hit only one marker gene, all candidate species are directly ranked from highest to lowest according to their comprehensive score; comprehensive score = number of matching sequences × average similarity; Specifically, during grouped statistical analysis, all sequences within the same species are uniformly merged and statistically analyzed: Marker gene types: The marker gene types in all sequences of the species are counted after deduplication. If two or more of the following types—COX1, Cytb, and D-loop—are present simultaneously, it is considered a multi-marker match; Number of matching sequences: The total number of all valid sequences belonging to any marker gene within the species is summed; Average similarity: The sum of the similarity values ​​of all valid sequences in the species is divided by the number of matching sequences to obtain the average similarity. Multi-marker gene cross-validation reduces false positives in species identification. If no species simultaneously matches two or more markers, the product of the number of matching sequences and the average similarity is used as the comprehensive score for ranking.

[0042] Standardized calculation formula for composite score: Score = N×AvgPident; Where N represents the number of valid matching sequences for the species, and AvgPident represents the arithmetic mean of the similarity of all sequences for the species. Among candidate species with the same number of marker genes, the higher the comprehensive score, the higher the species’ host credibility ranking. If all candidate species only match a single marker gene, then the comprehensive score is used to rank all species.

[0043] S6 outputs multiple sets of candidate host information and generates a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

[0044] Specifically, a standardized host identification result table is output uniformly, and all comparison and classification information, such as query sequence ID, accession number, similarity, comparison expectation value, species classification number, and species name, is fully retained.

[0045] Figure 2 This is a flowchart of the classification database and tag keyword preloading process according to an embodiment of the present invention. This pre-processing is completed before step S2. The specific process is as follows: The first step is to download the taxdmp classification data package from the NCBI Taxonomy database and compress it into two core text files: nodes.dmp and names.dmp. Each record in nodes.dmp stores the species classification number, the parent node species classification number, the taxonomic rank, and the taxonomic domain identifier, which is used to build the classification hierarchy. Names.dmp records the species name corresponding to each species classification number, realizing a two-way mapping between species classification numbers and species names.

[0046] The taxonomic domain identifier is the boundary identifier generated by preprocessing the taxonomic nodes of the NCBI Taxonomy database in this invention. The numerical correspondence rules are: 1=bacteria, 2=archaea, 3=fungi, 4=plants, 5=animals. This method only parses animal taxonomic nodes with taxonomic domain identifier = 5 (i.e., the Metazoa taxonomic node of the animal kingdom, root classification number = 33208), removes other boundary data to reduce storage size, and simultaneously constructs four types of resident memory mapping structures to collaboratively complete verification and backtracking. taxid_ptaxid dictionary: key = species classification number, value = parent node species classification number, tracing the classification hierarchy level by level; taxid_lineage dictionary: Key = species classification number, Value = kingdom / phylum / class / order / family / genus level; The nodes_set collection stores all animal species classification numbers, allowing for quick filtering of non-animal sequences using a query method with algorithmic time complexity. taxid_name dictionary: key = species classification number, value = species name, used to concatenate pedigree text; The four types of structures reside entirely in memory, eliminating the need for repeated readings of local files, thus reducing memory usage by over 90% and increasing classification and tracing speed by tens of times. Taking the NCBI Taxonomy 2024 database as an example, the full set of classification nodes contains approximately 2.3 million species classification numbers. Loading all nodes to build the memory dictionary consumes approximately 280MB of memory, with the initial database load taking approximately 4.2 seconds. Using the on-demand loading mechanism of this invention, when loading only the animal kingdom (species classification number = 33208) classification nodes, approximately 1.2 million nodes are loaded, reducing memory usage to approximately 145MB, a reduction of approximately 48% compared to full loading; when loading only the bird class (species classification number = 8782) classification nodes, approximately 42,000 nodes are loaded, reducing memory usage to approximately 18MB, a reduction of approximately 93.5% compared to full loading; and for targeted research on Anseriformes birds (species classification number = 8838), only approximately 1,800 classification nodes need to be loaded, reducing memory usage to approximately 1.2MB, a reduction of approximately 99.6% compared to full loading. In terms of classification and tracing speed, the average time for classification and tracing a single sequence in the full database is approximately 1.2 milliseconds, while the average time is reduced to approximately 0.3 milliseconds when only Aves nodes are loaded, representing a speed improvement of approximately 4 times. Considering the savings in database preloading time, the overall analysis time for a single sample is reduced from approximately 5.1 seconds in full mode to approximately 0.9 seconds in on-demand mode, resulting in an overall analysis efficiency improvement of more than 5 times. For projects studying only specific order-level taxa, the preloading time is reduced from 4.2 seconds to less than 0.1 seconds, representing a loading speed improvement of more than 40 times.

[0047] The second step is to parse the names.dmp file, extract the Latin and Chinese names of all animal species corresponding to their classification numbers, and construct a dictionary mapping species classification numbers to species names.

[0048] The third step is to construct a mitochondrial marker keyword library, which includes all keyword variants corresponding to COX1, Cytb, and D-loop, adapts to different database annotation formats, and avoids missing matches of target sequences.

[0049] Once the classification database and the tagging keyword library are built, they reside in memory and do not need to be loaded repeatedly when executing the S1-S6 analysis process.

[0050] Figure 3 This is a flowchart of mitochondrial marker gene screening and three-level threshold filtering provided in an embodiment of the present invention, as follows: Figure 3 As shown, it includes: Read the nucleotide homology alignment results file in the 6th standard output format line by line, skipping the comment lines and lines with incomplete format, and extract the key fields of each alignment hit: query sequence number, target sequence accession number, similarity percentage, alignment length, alignment expectation, query sequence length, species classification number, and target sequence title.

[0051] The first step is to screen mitochondrial marker genes: Case-insensitive keyword matching is performed on the title field of each matched target sequence. For the COX1 gene, the keywords include cox1, coi, cytochromecoxidasesubuniti, cytochromeoxidasesubunit1, and co1; for the Cytb gene, the keywords include cytb, cytochromeb, and cob; and for the D-loop region, the keywords include d-loop, controlregion, displacementloop, and dloop. Only alignments containing any one of these keywords in the title are retained; numerous irrelevant alignments not related to the target gene are directly filtered out.

[0052] The second step is to perform alignment length threshold filtering: the alignment length can be adjusted from 300 to 1000 base pairs, with a default threshold of 500 base pairs. Short sequences below the threshold are directly rejected. 500 base pairs conforms to the standard fragment length of DNA barcodes, which can avoid misidentification of closely related species caused by short conserved regions. For severely degraded samples (such as environmental samples and old specimens), the threshold can be lowered through the --min-length parameter.

[0053] The third step is to perform sequence similarity threshold filtering: extract the similarity percentage field `pident`. If the similarity is less than the user-defined similarity threshold, the match is filtered out. The adjustable similarity range is 95%~99%, with a default threshold of 97%. This threshold is set based on the DNA barcoding gap principle: in animal species, the mitochondrial sequence difference between individuals of the same species is usually less than 3%, while the sequence difference between different species is usually greater than 3%. Therefore, 97% is an empirical threshold for animal species identification. For identification of closely related species, users can increase the threshold to 98% or 99% using the `--min-ident` parameter, and decrease it to 95% for preliminary screening.

[0054] The fourth step involves comparison coverage threshold filtering: The comparison coverage is calculated based on the comparison length and the query sequence length. If the coverage is less than the user-defined coverage threshold, the match is filtered out. The default coverage threshold is 60% to avoid false positives caused by matching short, conservative regions. Only sequences that pass all three levels of threshold filtering are retained as valid candidate sequences and proceed to the subsequent classification and tracing process.

[0055] Figure 4 This is a flowchart of the automatic tracing process for classification lineages provided in an embodiment of the present invention, such as... Figure 4 As shown, it includes: For each valid candidate sequence that passes the three-level threshold filtering, its corresponding species classification number is first extracted. It is then determined whether the species classification number exists in the preloaded animal kingdom node set nodes_set. If it does not exist, it is determined to be a non-animal sequence (such as bacterial, fungal, or plant sequences) and is removed.

[0056] For species classification numbers existing in the Animal Kingdom set, initiate a taxonomic tracing loop: initialize the current species classification number as the sequence's staxid, and initialize the tracing string to empty; starting from the sequence species classification number, query the parent node level by level through the taxid_ptaxid dictionary, and at each level upward, add the species name and taxonomic rank of the current level to the beginning of the tracing string in the fixed format "species name (taxonomic rank);"; the user can configure the root classification number as the tracing termination node, the default root species classification number = 33208 (Animal Kingdom). 33208 is the built-in root standard species classification number for the animal kingdom in the NCBI Taxonomy database. Users can customize and input species classification numbers at the class, order, and family levels as the tracing termination node according to their research needs. The program stops tracing when it reaches the custom node, retaining only species within the target higher-order unit and automatically filtering out other taxa sequences. After reaching the specified node, tracing stops, and a complete phylogenetic string is generated. After each tracing, it checks whether the current species classification number is equal to the root node species classification number specified by the user (the default is the animal kingdom species classification number = 33208). If they are equal, tracing stops, root node information is added to the phylogenetic header, and a complete taxonomic phylogenetic string is generated. If a species classification number is missing or the root node cannot be reached during the tracing process, the phylogenetic is considered incomplete, and the sequence is removed.

[0057] This classification and tracing algorithm is implemented entirely in memory through dictionary lookups, without requiring network access to call the NCBI API. The tracing time for a single sequence is less than 1 millisecond, which is hundreds of times faster than online classification queries, and can quickly process tens of thousands of valid sequences.

[0058] Figure 5 This is a flowchart of the host species determination and statistical ranking output provided in the embodiments of the present invention, such as... Figure 5 As shown, it includes: After all valid sequences have been classified and traced, the valid candidate sequences are grouped by species name, and three core indicators are calculated for each candidate species: the number of marker gene species that hit, the number of matching sequences, and the average similarity.

[0059] When ranking candidate species: first, sort them from highest to lowest by the number of types of marker genes hit, with those having more types ranked higher; species with ≥2 types of marker genes hit are considered high-confidence hosts; when different species have the same number of types of marker genes hit, they are ranked from highest to lowest by their overall score; when all candidate species hit only one type of marker gene, all candidate species are directly ranked from highest to lowest by their overall score; overall score = number of matching sequences × average similarity; the species with the highest score is selected as the most likely host.

[0060] The species with the highest overall ranking is identified as the optimal host for the sequencing library. The top 5 candidate hosts and all statistical indicators are also output for researchers to manually review. Finally, a standardized identification table is generated, which includes the query sequence ID, target sequence accession number, type of hit marker gene, similarity percentage, alignment length, alignment coverage, alignment expectation, species classification number, species name, and complete taxonomic lineage. The output data can be directly used for virus tracing and host interaction analysis.

[0061] Figure 6 This is a block diagram of the high-throughput sequencing library host rapid identification system provided in the embodiments of the present invention, such as... Figure 6 As shown, the identification system used to implement the aforementioned identification method consists of seven interconnected functional modules: database construction module, data input module, gene screening module, threshold filtering module, classification and tracing module, host determination module, and result output and statistics module. The division of labor and execution flow of each module are as follows: The database construction module is used to perform the pre-database automatic construction step. This pre-database automatic construction process needs to be completed before the gene screening module performs the sequence matching operation. The program generates a classification database containing four types of memory mappings and resident in memory by screening animal classification nodes with classification domain identifier = 5, and simultaneously constructs a tag keyword library that includes all annotation variants of COX1, Cytb, and D-loop. The data input module is used to execute the acquisition of alignment files with the nucleotide homology alignment output of contig sequences and the field structure of the sixth type of standard output format. The gene screening module is used to read the alignment file, match keywords corresponding to the mitochondrial DNA barcode markers of three types of animals (COX1, Cytb, and D-loop), and screen to obtain target gene matching sequences containing only the target marker. The threshold filtering module is used to perform a three-layer linkage filtering on the target gene matching sequence in sequence, which considers alignment length, sequence similarity, and alignment coverage, to obtain an effective candidate sequence that simultaneously meets the three-layer threshold constraints. The classification and tracing module is used to extract the species classification number of the valid candidate sequence, verify the species classification number to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back upwards step by step, eliminating non-animal sequences and sequences with incomplete system, and obtaining a complete system analysis sequence. The host determination module is used to group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity in each group; species with ≥2 marker gene types are identified as high-confidence hosts; when ranking candidate species, they are first arranged from highest to lowest by the number of marker gene types, with those having more types ranked higher; when different species have the same number of marker gene types, these species are ranked from highest to lowest by their comprehensive score; when all candidate species have only one marker gene, all candidate species are directly ranked from highest to lowest by their comprehensive score; comprehensive score = number of matching sequences × average similarity; The result output and statistics module is used to output multiple sets of candidate host information and generate a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

[0062] Each module operates automatically and collaboratively, requiring no manual intervention throughout the entire process, fully implementing all the steps of the identification method of this invention, and corresponding one-to-one with the program execution logic recorded in the identification system.

[0063] This invention also provides an electronic device for rapid host identification of high-throughput sequencing libraries. This device has a built-in memory and processor. The memory stores a computer program that can be called and executed by the processor. When the processor executes the computer program, it can implement the aforementioned rapid host identification method, sequentially completing the entire process of database preloading, mitochondrial marker sequence screening, three-level threshold filtering, classification lineage tracing, host statistical determination, and standardized table output. Figure 6 The seven functional modules work together to complete the automated host identification process.

[0064] The present invention further provides a computer-readable storage medium, which includes any medium with data storage capability such as a USB flash drive, hard disk, solid-state drive, ROM, RAM, optical disk, etc., and stores a computer program on the medium; when the computer program is read and executed by the processor of an electronic device, the aforementioned method for rapid identification of high-throughput sequencing library hosts is realized, which can complete the identification and statistical output of sequencing library hosts in batches without manual operation.

Claims

1. A method for rapid host identification of high-throughput sequencing libraries, characterized in that, It includes an automatic pre-database construction step and an automated analysis step consisting of S1-S6; the automatic pre-database construction step is completed before S2 is executed. The program generates a classification database containing four types of memory-mapped structures and resident in memory by filtering animal classification nodes with classification domain identifier = 5, and simultaneously constructs a tag keyword library that includes all annotation variants of COX1, Cytb, and D-loop. The automated analysis steps include: S1, Obtain the alignment file with the field structure of the Class 6 standard output format, which is the output of the contiguous group sequence after nucleotide homology alignment; S2, read the comparison file, match the keywords corresponding to the mitochondrial DNA barcode markers of three types of animals: COX1, Cytb, and D-loop, and screen to obtain target gene matching sequences containing only the target marker; S3, perform a three-layer linkage filtering on the target gene matching sequence in sequence, which includes alignment length, sequence similarity and alignment coverage, to obtain a valid candidate sequence that simultaneously meets the three-layer threshold constraints; S4. Extract the species classification number of the valid candidate sequence, verify the classification number to which the species classification number belongs to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back upwards step by step. Remove non-animal sequences and sequences with incomplete classification systems to obtain a complete classification sequence. S5, group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity in each group; species with ≥2 marker gene types are identified as high-confidence hosts; when ranking candidate species, they are first arranged from highest to lowest according to the number of marker gene types hit, with those having more types ranked higher; when different species have the same number of marker gene types hit, these species are ranked from highest to lowest according to their comprehensive score; when all candidate species hit only one marker gene, all candidate species are directly ranked from highest to lowest according to their comprehensive score; comprehensive score = number of matching sequences × average similarity; S6 outputs multiple sets of candidate host information and generates a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

2. The method for rapid host identification of a high-throughput sequencing library according to claim 1, characterized in that, Before executing step S2, the pre-construction database automatic construction step is completed, constructing a classification database and a tag keyword library; the classification database only stores animal classification nodes with classification domain identifier = 5, including four types of memory mapping structures: species classification number and parent node, species classification number and classification rank, animal species classification number, and species classification number and species name; the classification database and tag keyword library are resident in memory after construction; when tracing the classification lineage upwards step by step in step S4, the classification database is retrieved to complete the species classification number attribution verification, and the root species classification number is configured as the tracing termination node, with the root species classification number having a default value of 33208; the tag keyword library includes various annotation variants such as COX1, Cytb, and D-loop.

3. The method for rapid host identification of a high-throughput sequencing library according to claim 2, characterized in that, In the keyword library, the keywords corresponding to COX1 include cox1, coi, co1, cytochromecoxidasesubuniti, and cytochromeoxidasesubunit1; the keywords corresponding to Cytb include cytb, cytochromeb, and cob; and the keywords corresponding to D-loop include d-loop, dloop, controlregion, and displacementloop.

4. The method for rapid host identification of a high-throughput sequencing library according to claim 1, characterized in that, In step S3, the alignment length is adjustable in the range of 300 to 1000 base pairs, the sequence similarity is adjustable in the range of 95% to 99%, and the alignment coverage is adjustable in the range of 50% to 90%. The alignment coverage is the percentage obtained by dividing the alignment length by the length of the complete sequence of the contiguous group.

5. The method for rapid host identification of a high-throughput sequencing library according to claim 4, characterized in that, The default alignment length is 500 base pairs, the default sequence similarity is 97%, and the default alignment coverage is 60%.

6. The method for rapid host identification of a high-throughput sequencing library according to claim 1, characterized in that, The standardized identification form described in step S6 includes the query sequence ID, target sequence accession number, type of hit marker gene, sequence similarity, alignment length, alignment coverage, alignment expectation, species classification number, species name, and complete taxonomic lineage.

7. A rapid host identification system for high-throughput sequencing libraries, characterized in that, It includes a database construction module, a data input module, a gene screening module, a threshold filtering module, a classification and tracing module, a host determination module, and a result output and statistics module; The database construction module is used to perform the automatic construction of the front-end database. By filtering animal classification nodes with classification domain identifier = 5, it generates a classification database containing four types of memory mappings and resident in memory, and simultaneously constructs a tag keyword library that includes all annotation variants of COX1, Cytb, and D-loop. The data input module is used to acquire the alignment file with the field structure of the sixth type of standard output format, which is the output of the contiguous group sequence after nucleotide homology alignment. The gene screening module is used to read the alignment file, match keywords corresponding to the mitochondrial DNA barcode markers of three types of animals (COX1, Cytb, and D-loop), and screen to obtain target gene matching sequences containing only the target marker. The threshold filtering module is used to perform three-layer linkage filtering on the target gene matching sequence in sequence, namely, alignment length, sequence similarity and alignment coverage, to obtain effective candidate sequences that simultaneously meet the three-layer threshold constraints. The classification and tracing module is used to extract the species classification number of the valid candidate sequence, verify the classification number to which the species classification number belongs to the classification database, and generate a complete hierarchical classification system of kingdom, phylum, class, order, family, genus and species by tracing back level by level, eliminating non-animal sequences and sequences with incomplete system, and obtaining a complete system analysis sequence. The host determination module is used to group the complete phylogenetic sequence by species, and count the number of marker gene types, the number of matching sequences, and the average similarity of each group; species with ≥2 marker gene types are determined as high-confidence hosts; when ranking candidate species, they are first arranged from high to low according to the number of marker gene types, with those having more types ranked higher; when different species have the same number of marker gene types, these species are arranged from high to low according to their comprehensive score; when all candidate species hit only one marker gene, all candidate species are directly arranged from high to low according to their comprehensive score. Overall score = number of matching sequences × average similarity; The result output and statistics module is used to output multiple sets of candidate host information and generate a standardized identification table containing sequence, alignment parameters, species, and complete taxonomic lineage.

8. An electronic device for rapid host identification of high-throughput sequencing libraries, characterized in that, The device stores a computer program, which executes a method for rapid host identification of a high-throughput sequencing library according to any one of claims 1 to 6 when the program is run.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the method for rapid host identification of a high-throughput sequencing library as described in any one of claims 1 to 6.

Images