Virus identification method in NGS sequence based on multi-dimensional feature fusion

Abstract

The present application relates to the technical field of virus identification, in particular to a virus identification method in NGS sequence based on multi-dimensional feature fusion, comprising the following steps: S1, non-parametric assembly analysis and non-parametric feature extraction; S2, sample-specific reference sequence library construction; S3, parametric coverage analysis and feature extraction; S4, machine learning model design and pathogen accurate determination; S5, virus identification in NGS sequence. The present application performs non-parametric assembly analysis on NGS sequence, constructs a reference sequence for parametric comparison according to the results obtained by non-parametric assembly, carries out three kinds of parametric coverage analysis based on the constructed reference sequence library, finally designs a machine learning algorithm, constructs a training label based on qPCR verification data, screens sequence features and trains a model, uses the trained model as the final determination basis, and realizes accurate virus identification.

Description

Technical Field

[0001] This invention relates to the field of virus identification technology, and in particular to a method for virus identification in NGS sequences based on multidimensional feature fusion. Background Technology

[0002] Waterfowl farming is a crucial pillar of my country's agricultural economy, and its healthy development is vital for food safety and industrial economic stability. In core farming areas nationwide, high-density, open-air polyculture is widely adopted. While this model suits large-scale production, it presents significant disease risks: waterfowl have wide activity ranges and share water bodies, making pathogen spread during drinking and feeding; furthermore, the farming environment often overlaps with migratory bird routes, and insufficient biosecurity measures further accelerate cross-species viral transmission, seriously threatening the sustainable development of the industry. Under high-density polyculture, waterfowl diseases exhibit diverse and overlapping characteristics. Epidemiologically, infectious diseases are predominant, encompassing multiple pathogens such as viruses, bacteria, and mycoplasma, with viruses being the primary causative agent. These viruses exhibit diversity, rapid mutation, and the coexistence of multiple serotypes. More significantly, atypical infections, mixed infections, and secondary infections are common in clinical practice. The synergistic effect of multiple pathogens leads to complex symptoms, difficult diagnosis, and increased mortality. Simultaneously, the prevalence of immunosuppressive viroids and enhanced cross-host transmission capabilities further amplify the difficulty of prevention and control, rendering traditional "single-pathogen targeted prevention and control" approaches inadequate. Therefore, timely and accurate identification of various pathogens in the aquaculture environment and clarification of infection combinations are crucial for disease prevention and control. Strengthening the detection process can not only enable early detection and early warning, buying time for prevention and control measures, but also provide comprehensive data support for research on epidemic patterns and variation trends, promoting a shift in prevention and control efforts from passive response to proactive prevention.

[0003] Currently, there are three main methods for pathogen detection. The first relies on traditional PCR / qPCR, which uses primers designed based on known conserved regions of pathogens. This method cannot detect variants, rare or unknown pathogens, and its single detection range is limited, making multi-stage experiments cumbersome, time-consuming, and labor-intensive. The second method involves pathogen isolation and culture, but its clinical application is limited due to its complexity, long cycle, and low detection efficiency for low-abundance pathogens. The third method uses high-throughput sequencing, parameter-free assembly, and parameterized coverage analysis for pathogen discovery and the mining of unknown pathogens. The main related technical solutions currently include: 1. Chinese Patent Application No. 201610408646.2, this method targets the mixed infections that often occur in waterfowl. By designing specific primers for the conserved regions of different viruses VP1 in the same reaction system, it solves the limitation that a single experiment can only detect a single target and simplifies the cumbersome multi-group experimental process. 2. Chinese Patent Application No. 202010914568.X, this patent optimizes primers and probes for clinical low-abundance samples, avoiding cross-reaction with other common waterfowl viruses (such as goose reovirus, avian Escherichia coli, etc.), and greatly improves the detection limit of known targets.

[0004] 3. Chinese Patent Application No. 202210732531.4 addresses the pain points of host background interference in metagenomic sequencing (mNGS) and the ease with which single parameterless or parameterized assembly can lead to sequence breaks or omissions of variation information. It proposes to integrate a dual-track strategy of "parameterless assembly + parameterized assembly" to successfully construct a complete genome sequence that is free from host contamination and contains both point mutation and structural variation information, thus achieving a technological upgrade from "massive fragmented short sequences" to "highly accurate identification and discrimination of pathogenic microorganisms".

[0005] 4. The paper “The evolutionary history of vertebrate RNA viruses” (Zhang Yongzhen et al., Nature, 2018) mainly utilizes metatranscriptome sequencing (mNGS) technology combined with parameterless assembly algorithm to achieve large-scale mining and identification of 214 novel and unknown RNA viruses in a large number of lower vertebrate samples.

[0006] However, based on the above technical solutions, existing technologies still face problems such as strong background noise, extremely low target abundance, and high false positive rates. While Chinese Patent Application No. 202210732531.4 also has pathogen detection capabilities, its parameterless viral assembly method is too simplistic, resulting in a high false positive rate. If only a parameterized coverage method is used, the parameterized coverage analysis relies on a reference sequence, and the quality of the selected reference sequence directly affects the virus identification results. Chinese Patent Application No. 202010914568.X targets low-abundance clinical samples, but in complex waterfowl samples, the background nucleic acid abundance of the host and symbiotic microorganisms is extremely high, while the viral content is very low, exhibiting a significant difference in "abundance scale," making the screening and discovery of low-abundance unknown pathogens extremely difficult.

[0007] Therefore, avoiding the limitations of existing methods, reducing false positives, and improving the reliability of virus detection are crucial for the field of virus identification in complex waterfowl samples. Summary of the Invention

[0008] This invention provides a method for virus identification in NGS sequences based on multidimensional feature fusion. Its purpose is to solve the core problems of existing technologies such as "strong background noise, extremely low target abundance, and high false positive rate", and to provide a new solution for highly reliable identification of pathogens in complex samples.

[0009] To achieve the above objectives, this invention provides a method for virus identification in NGS sequences based on multidimensional feature fusion, comprising the following steps: S1. Parameterless assembly analysis and parameterless feature extraction: Through parameterless assembly and homology annotation, the original sequencing data is deeply cleaned, redundancy is removed and spliced, and the counting features and quality features of the sequences are extracted as parameterless features. S2. Construction of Sample-Specific Reference Sequence Library: Based on the parameterless assembly and homology annotation results of step S1, for suspected animal pathogenic viruses, the corresponding contigs are extracted and compared with the viral genome sequences in the nucleic acid sequence database. Based on the homology and query coverage indicators, the full-length viral genome or core gene sequence that is closest in evolutionary relationship and has the highest matching degree with the parameterless assembly sequence is selected from the database to construct the sample-specific reference sequence library. S3. Parametric Coverage Analysis and Feature Extraction: Based on the sample-specific reference sequence library constructed in S2, multi-dimensional feature extraction is performed to obtain parametric multi-dimensional features; the multi-dimensional features include one or more of the following: Reads dimension feature extraction, Contigs dimension feature extraction, or fusion dimension feature extraction. S4. Machine learning model design and accurate pathogen identification: Integrate non-parametric features and parametric multi-dimensional features, and calculate derived features based on quality features. Integrate non-parametric features, parametric multi-dimensional features and quality features to construct a multi-dimensional quantization matrix, and build a machine learning model through supervised learning. S5. Based on the machine learning model constructed in step S4, identify viruses in NGS sequences.

[0010] Preferably, the deep cleaning, redundancy removal, and splicing of the raw sequencing data in step S1 specifically includes: Data quality control and cleaning: To ensure downstream assembly quality, data quality control software such as Trimmomatic was used to perform quality control on the raw data. A sliding window filtering strategy was adopted to remove low-quality bases and linker sequences. The sliding window size was set to a range of 3 to 5 bases, and the average quality cutoff value ranged to 15 to 25. Short fragments with a length less than a preset minimum length threshold, ranging from 20 bp to 50 bp, were also removed to obtain clean reads. Host sequence stripping: Clean reads were aligned to the host reference genome using sequence alignment software such as Bowtie2, and unmapped reads were extracted using end-to-end alignment mode to eliminate host background interference; Low-abundance sequence capture and redundancy removal: Unmapped reads were assembled from scratch using assembly software such as MEGAHIT, employing a multiple k-mer iteration strategy with k-mer values ​​ranging from 21 to 119. Contigs with a length not less than a preset length threshold ranging from 50 bp to 200 bp were retained. Subsequently, redundancy removal was performed by clustering using tools such as CD-HIT-EST with a set sequence consistency threshold ranging from 0.90 to 0.99, constructing a globally non-redundant candidate set. The counting and quality features of the synchronously extracted sequences specifically include: Species annotation and feature quantification: BLASTx alignment was performed using homology sequence retrieval tools such as DIAMOND. The upper limit of the expected value (E-value) was set to range from 1e-5 to 1e-15, the lower limit of identity was set to range from 30% to 70%, and the lower limit of alignment length was set to range from 30 bp to 100 bp. Species phylogenetics were reconstructed using taxonomic tools such as TaxonKit, removing bacteriophages and plant and fungal viruses while retaining candidate sequences of animal pathogenic viruses. During this process, the system groups samples by name and virus name, and automatically extracts the following parameter-free features: Counting features include one or more of the following: Total Contigs, Correct Contigs, or Noise Contigs. Quality features: For a subset of correctly annotated data, aggregate and calculate one or more of the following: maximum / mean homology (identity_max / identity_mean), maximum / mean / sum score (BitScore_max / BitScore_mean / BitScore_sum), maximum / sum alignment length (AlignLength_max / AlignLength_sum), average query sequence length (qlen_mean), minimum expected value (Evalue_min), or average mismatches (Mismatches_mean).

[0011] Preferably, the sequence alignment in step S2 is performed using NCBI's online BLASTn tool.

[0012] Preferably, the Reads dimension feature extraction in step S3 specifically includes: aligning Clean reads to the reference library, parsing the output results and extracting the initial aggregation parameters, including one or more of the following: Reads_Count_reads, Coverage_bp_reads, Coverage_Percent_reads, Avg_Coverage_Depth_x_reads, Mean_BaseQ_reads, Mean_MapQ_reads, or Histo_Max_Bin_Percent_reads, which represents the maximum frequency interval of the coverage depth histogram. Contig dimensional feature extraction specifically includes: aligning non-redundant contigs to the same reference library, and similarly extracting the corresponding 7 parameters: total number of contig alignment reads (Reads_Count_contigs), total length of bases covered by the contig (Coverage_bp_contigs), percentage of genome coverage by the contig (Coverage_Percent_contigs), average sequencing depth of the contig (Avg_Coverage_Depth_x_contigs), average alignment quality value of the contig (Mean_MapQ_contigs), and the percentage of the largest interval in the contig coverage histogram (Histo_Max_Bin_Percent_contigs), to assess the macroscopic continuity of the assembled sequence; The fusion dimension feature extraction specifically includes: merging the BAM files of Reads and Contigs, calculating the comprehensive coverage under the superimposed state, and extracting the corresponding 7 feature parameters of the fusion sequence: including the total number of aligned reads (Reads_Count_fused), the total length of bases covered by the fusion sequence (Coverage_bp_fused), the coverage percentage of the fusion sequence (Coverage_Percent_fused), the average coverage depth of the fusion sequence (Avg_Coverage_Depth_x_fused), the average base quality value of the fusion sequence (Mean_BaseQ_fused), the average alignment quality value of the fusion sequence (Mean_MapQ_fused), and the maximum interval percentage of the coverage histogram of the fusion sequence (Histo_Max_Bin_Percent_fused), reflecting the comprehensive mapping strength of the global dimension.

[0013] Preferably, the calculation formula for derived features based on quality features in step S4 is as follows:

[0014] AlignLength_max=max{AlignLength_i|Right_Annotated_Contigs_i =1} qlen_mean=(1 / n)×Σ(qlen_i) for i where Right_Annotated_Contigs_i=1 in, Align_Qlen_Ratio_Max: Optimal coverage ratio; AlignLength: The maximum alignment length among all correctly annotated contigs within the group; i iterates through all correct comments (Contig) within the group; Biological meaning: Alignment length refers to the number of base pairs that successfully match between the query sequence (contig) and the reference sequence (viral reference genome). It reflects: The length of the homologous region between this contig and the reference virus; A longer AlignLength usually means that the contig covers a larger segment of the reference viral genome; In virus testing, the larger the AlignLength, the higher the reliability of the test results.

[0015] qlen: The arithmetic mean of the query sequence lengths of all correctly annotated contigs within the group; n is the total number of correctly annotated contigs in this group.

[0016] Biological meaning: The query sequence length refers to the total length (number of bases) of the contig itself. qlen_mean reflects: The average size of all correctly annotated contigs in this sample-virus pair. A larger qlen_mean usually indicates better contig assembly quality and more complete fragments. A smaller qlen_mean may indicate the presence of short fragments, fragmented assembly, or insufficient sequencing depth. (Note: A small value is introduced in the denominator to avoid division by zero error).

[0017] At this point, a multidimensional feature matrix containing 35 feature variables has been fully constructed.

[0018] Preferably, step S4, which involves calculating and integrating derived features from non-parametric features, parametric multi-dimensional features, and quality features to construct a multi-dimensional quantization matrix, and then constructing a machine learning model through supervised learning, specifically includes: Data preprocessing and labeling: qPCR validation results were introduced as the gold standard, and binary labels were assigned: positive = 1 and negative = 0. The feature matrix was subjected to outlier removal, median imputation and standard deviation normalization in sequence, and then the training set and test set were divided in a 7:3 ratio. Model training and evaluation: Random forest (RF), support vector machine (SVM), k-nearest neighbors (KNN), logistic regression (LR), decision tree (DT), and gradient boosting tree (GBDT) were used to construct classification models in parallel. Grid search combined with 5-fold cross-validation was used to optimize hyperparameters.

[0019] Preferably, in step S5, the machine learning model constructed in step S4 is used to identify viruses in NGS sequences. The input is the virus sample to be detected, and the output is the target pathogen of the sample to be detected.

[0020] Preferably, the virus identification method is used for identifying low-abundance viruses in complex waterfowl samples.

[0021] The above-described solution of the present invention has the following beneficial effects: (1) This invention provides a machine learning algorithm for virus identification that integrates non-parametric assembly features and parametric coverage features, for accurate detection of viruses from NGS high-throughput sequencing samples. Non-parametric assembly analysis is performed using NGS sequences. Reference sequences with reference alignment are constructed based on the results of non-parametric assembly. Three types of parametric coverage analysis are carried out based on the constructed reference sequence library, including reads coverage analysis, contigs coverage analysis, and multi-dimensional feature fusion coverage analysis. Finally, a machine learning algorithm is designed. This algorithm is based on qPCR validation data to construct training labels, screen sequence features and train the model. The trained model is used as the final judgment criterion to achieve accurate virus identification.

[0022] (2) The virus identification method in NGS sequences based on multidimensional feature fusion provided by this invention is constructed based on multidimensional features. Even when faced with novel clinical samples, it still maintains extremely high pathogen discrimination stability and robustness. This method completely eliminates the dependence of traditional single NGS analysis processes on a large number of downstream qPCR validations, and realizes fully automated intelligent identification from sequencing data to accurate pathogen identification, which has extremely high clinical application and industrial transformation value.

[0023] (3) The method of this invention achieves accurate detection of viruses from high-throughput sequencing samples by integrating parameterless assembly features and parameterized coverage features for virus identification. This avoids false negatives caused by parameterless assembly alone, improves the accuracy of virus identification, and is a highly reliable identification scheme suitable for low-abundance viruses in complex samples of high-abundance waterfowl. The experimental results above show that the method of this invention achieves a recall rate of 97.27% and an AUC value of 0.9877 for virus identification. Attached Figure Description

[0024] Figure 1 This is a flowchart of the virus identification method in NGS sequences based on multidimensional feature fusion according to the present invention; Figure 2 This is a comparison chart of the receiver operating characteristic (ROC) curves of the multiple machine learning classification model in Embodiment 1 of the present invention; Figure 3 This is a comparison chart of the F1 scores of the internal test set of the multiple machine learning classification model in Embodiment 1 of the present invention; Figure 4 This is a comparison chart of the average scores of the 5-fold cross-validation (CV) of the multiple machine learning classification model in Embodiment 1 of the present invention; Figure 5 This is the confusion matrix of the internal test set prediction of the optimal classification model (random forest) in Embodiment 1 of the present invention; Figure 6 This is the confusion matrix of the external independent test set prediction for the optimal classification model (random forest) in Embodiment 1 of the present invention. Figure 7 This is the receiver operating characteristic (ROC) curve of the external independent test set for the optimal classification model (random forest) in Embodiment 1 of the present invention. Detailed Implementation

[0025] To make the technical problems, solutions, and advantages of this invention clearer, a detailed description will be provided below with reference to the accompanying drawings and specific embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

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

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

[0028] Furthermore, the technical features involved in the different embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0029] Example 1: A method for virus identification in NGS sequences based on multidimensional feature fusion includes the following steps: S1. Parameterless assembly analysis and parameterless feature extraction: Through parameterless assembly and homology annotation, the original sequencing data is deeply cleaned, redundancy is removed and spliced, and the counting features and quality features of the sequences are extracted as parameterless features. Data quality control and cleaning: To ensure the quality of downstream assembly, Trimmomatic (v0.39) was used to perform quality control on the raw data. A sliding window filter (SLIDINGWINDOW:4:20) was used to remove low-quality bases and adapter sequences (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10), and short fragments less than 30 bp in length (MINLEN:30) were removed to obtain clean reads.

[0030] Host sequence stripping: Using Bowtie2 (v2.4.5), clean reads were aligned to the host reference genome (such as Anser cygnoides) of waterfowl, and unmapped reads were extracted using end-to-end alignment mode to eliminate host background interference.

[0031] Low-abundance sequence capture and redundancy removal: Unmapped reads were assembled de novo using MEGAHIT (v1.2.9) with a multiple k-mer iteration strategy (k=21, 33, 55, 77) to retain contigs ≥100 bp in length. Subsequently, redundancy removal was performed using CD-HIT-EST (v4.8.1) with a consistency threshold of 0.95, constructing a globally non-redundant candidate set.

[0032] Species annotation and feature quantification: BLASTx alignment was performed using DIAMOND (v2.0.15) (E-value ≤1e-10, Identity ≥ 50%, Alignment Length ≥ 50 bp). Species lineage was reconstructed using TaxonKit, removing bacteriophages and plant and fungal viruses while retaining candidate sequences of animal pathogenic viruses.

[0033] During this process, the system groups samples by name and virus name, and automatically extracts the following parameter-free features: Counting features include the total number of contigs (Count_Total_Contigs), the number of correctly annotated contigs (Count_Correct_Contigs), and the number of contigs representing noise (Count_Noise_Contigs).

[0034] Quality characteristics: For a subset of correctly annotated data, aggregate and calculate the maximum / mean of homology (identity_max / identity_mean), the maximum / mean / sum of scores (BitScore_max / BitScore_mean / BitScore_sum), the maximum / sum of alignment lengths (AlignLength_max / AlignLength_sum), the average length of the query sequence (qlen_mean), the minimum expected value (Evalue_min), and the average number of mismatches (Mismatches_mean).

[0035] S2. Construction of a Sample-Specific Reference Sequence Library: Based on the parameter-free assembly and homology annotation results of step S1, for suspected animal pathogenic viruses, corresponding contigs are extracted, and sequence alignment is performed using NCBI's online BLASTn tool. Based on homology (Identity) and query coverage metrics, the full-length viral genome or core gene sequences that are most closely related to and have the highest matching degree to the parameter-free assembled sequence are precisely selected from the database. The FASTA files of the selected sequences are downloaded, and a "sample-specific reference sequence library" adapted to the current clinical samples is constructed as the absolute benchmark for subsequent multidimensional coverage quantification.

[0036] S3. Parametric Coverage Analysis and Feature Extraction: Based on the sample-specific reference sequence library constructed in S2, multi-dimensional feature extraction is performed to obtain parametric multi-dimensional features; the multi-dimensional features include one or more of the following: Reads dimension feature extraction, Contigs dimension feature extraction, or fusion dimension feature extraction. To address false positives caused by complex background interference and missed detections due to viral mutations, this step builds a "Reads-Contigs dual-dimensional quantization" and its fusion dimension coverage detection system based on the sample-specific reference sequence library constructed by S2, combined with Minimap2 and SAMtools.

[0037] Reads Dimension Feature Extraction: Align Clean reads to the reference library, parse the output results and extract the initial aggregation parameters: number of aligned reads (Reads_Count_reads), number of effective covered bases (Coverage_bp_reads), percentage of reference sequence coverage (Coverage_Percent_reads), average sequencing coverage depth (Avg_Coverage_Depth_x_reads), average base quality (Mean_BaseQ_reads), average sequence alignment quality (Mean_MapQ_reads), and percentage of the maximum frequency interval of the coverage depth histogram (Histo_Max_Bin_Percent_reads).

[0038] Contig dimensional feature extraction: Non-redundant contigs are aligned to the same reference library, and the corresponding 7 parameters are extracted: total number of contig alignment reads (Reads_Count_contigs), total length of bases covered by the contig (Coverage_bp_contigs), percentage of genome coverage by the contig (Coverage_Percent_contigs), average sequencing depth of the contig (Avg_Coverage_Depth_x_contigs), average alignment quality value of the contig (Mean_MapQ_contigs), and the percentage of the largest interval in the contig coverage histogram (Histo_Max_Bin_Percent_contigs) to assess the macroscopic continuity of the assembled sequence.

[0039] Fused (fusion) dimension feature extraction: Merge the BAM files of Reads and Contigs, calculate the overall coverage under the superimposed state, and extract the corresponding 7 feature parameters: total number of aligned reads in the fused sequence (Reads_Count_fused), total length of bases covered by the fused sequence (Coverage_bp_fused), percentage of coverage by the fused sequence (Coverage_Percent_fused), average coverage depth of the fused sequence (Avg_Coverage_Depth_x_fused), average base quality value of the fused sequence (Mean_BaseQ_fused), average alignment quality value of the fused sequence (Mean_MapQ_fused), and the maximum proportion of the coverage histogram of the fused sequence (Histo_Max_Bin_Percent_fused), which reflect the overall mapping strength of the global dimension.

[0040] S4. Machine learning model design and accurate pathogen identification: Integrate non-parametric features and parametric multi-dimensional features, and calculate derived features based on quality features. Integrate non-parametric features, parametric multi-dimensional features and quality features to construct a multi-dimensional quantization matrix, and build a machine learning model through supervised learning. Derived features are calculated based on quality characteristics, using the following formula:

[0041] AlignLength_max=max{AlignLength_i|Right_Annotated_Contigs_i =1} qlen_mean=(1 / n)×Σ(qlen_i) for i where Right_Annotated_Contigs_i=1 in, Align_Qlen_Ratio_Max: Optimal coverage ratio; AlignLength: The maximum alignment length among all correctly annotated contigs within the group; i iterates through all correct comments (Contig) within the group; Biological meaning: Alignment length refers to the number of base pairs that successfully match between the query sequence (contig) and the reference sequence (viral reference genome). It reflects: The length of the homologous region between this contig and the reference virus; A longer AlignLength usually means that the contig covers a larger segment of the reference viral genome; In virus testing, the larger the AlignLength, the higher the reliability of the test results.

[0042] qlen: The arithmetic mean of the query sequence lengths of all correctly annotated contigs within the group; n is the total number of correctly annotated contigs in this group.

[0043] Biological meaning: The query sequence length refers to the total length (number of bases) of the contig itself.

[0044] qlen_mean reflects: The average size of all correctly annotated contigs in this sample-virus pair; A larger qlen_mean usually means that the contig assembly quality is better and the fragments are more complete; A smaller qlen_mean may indicate the presence of short fragments, fragmented assembly, or insufficient sequencing depth; (Note: A small value is introduced in the denominator to avoid division by zero error).

[0045] At this point, a multidimensional feature matrix containing 35 feature variables has been fully constructed.

[0046] Data preprocessing and labeling: qPCR validation results were used as the gold standard and assigned binary labels (positive = 1, negative = 0). The feature matrix was then subjected to outlier removal, median imputation, and standard deviation normalization. The training and test sets were then divided in a 7:3 ratio.

[0047] Model training and evaluation: Random Forest (RF), Support Vector Machine (SVM), K-Nearest Neighbors (KNN), Logistic Regression (LR), Decision Tree (DT), and Gradient Boosting Tree (GBDT) algorithms were used to construct classification models in parallel. Grid Search combined with 5-fold Cross-Validation (CV) was employed to optimize hyperparameters.

[0048] This embodiment constructs a multi-class machine learning model based on the extracted 35-dimensional pathogen assessment feature matrix (including 14-dimensional non-parametric assembly features and 21-dimensional parametric mapping features), and conducts full-process performance verification and optimal model selection. The specific implementation steps are as follows: 1. Dataset Construction and Hierarchical Partitioning To ensure the rigor of the model generalization ability verification, this embodiment divides the 102 clinical samples (covering 510 sick geese) into the original modeling set and the external independent test set, with no sample overlap between the two.

[0049] The original modeling set contained 29 independent disease samples, from which 120 valid combination records were extracted (71 positive RT-qPCR records and 49 negative records, with a positive rate of approximately 59.2%). To avoid data bias, stratified random sampling was used to divide the set into an internal training set (84 records) and an internal validation set (36 records) in a 7:3 ratio, ensuring that the ratio of positive to negative samples in the two subsets remained at approximately 6:4.

[0050] External independent test set: contains 73 novel clinical samples that were not used in training, with a total of 179 valid records (110 positive and 69 negative), used only for final validation of generalization ability.

[0051] 2. Multi-algorithm modeling and hyperparameter optimization Using the standardized 35-dimensional feature matrix as input, six mainstream supervised learning algorithms were selected for simultaneous modeling: Random Forest (RF), Support Vector Machine (SVM), K Nearest Neighbors (KNN), Logistic Regression (LR), Decision Tree (DT), and Gradient Boosting Tree (GBDT).

[0052] To address the imbalance between positive and negative samples, class weight balancing (class_weight='balanced') is globally enabled in the RF, SVM, LR, and DT models. Based on the internal training set, a grid search algorithm combined with 5-fold cross-validation (CV) is used to iterate parameters with the F1 score as the core optimization objective. Parameter iteration optimization is completed within the internal training set, and the optimization range of the core hyperparameters for each model is shown in Table 1 below. Table 1 Machine Learning Models and Core Hyperparameters

[0053] 3. Model performance grading and evaluation and optimal model establishment Based on the above optimal hyperparameter configuration, the performance of the six models was evaluated on the internal validation set. Quantitative evaluation showed that different underlying algorithmic structures had significant differences in their ability to fit the high-dimensional feature matrix of the macrotranscriptome.

[0054] Optimal Model Establishment: Considering all indicators, the Random Forest algorithm, by constructing multiple decision trees and introducing a random feature selection mechanism, can handle the nonlinear mapping relationship between heterogeneous variables very well. With AUC = 0.984 and F1 = 0.957 on the test set, it is established as the optimal classification model in this embodiment.

[0055] The comparison results of multiple machine learning classification models are attached. Figure 2-5 As shown; Judgment and Analysis: Using precision, recall, F1 score, and AUC as core evaluation metrics, the optimal judgment model—Random Forest—is established. The optimal model after training will serve as the basis for judging unknown samples. Simultaneously, by outputting feature importance scores, the contribution of each dimension of features to pathogen identification decisions is analyzed, enhancing the model's interpretability.

[0056] To avoid the risk of model overfitting and to verify its universality in real clinical scenarios, this embodiment uses the aforementioned optimal classification model (Random Forest) and conducts actual data input and classification prediction on an external independent test set that was not involved in model training. The complete application process and verification results are as follows: 1. Data Input and Test Set Sources Test set source: This validation uses an external independent test set, which comes from 73 brand new independent clinical disease samples. After a full process of parameterless assembly (SO1) and parameterized coverage analysis (SO3), a total of 179 valid "sample-target virus" records were extracted.

[0057] Data structure: Each record is quantized into a 1×35 high-dimensional feature vector. RT-qPCR gold standard validation results are used as true alignment labels, including 110 true positive records and 69 true negative records.

[0058] 2. Example of single-sample data input and output (taking a single call as an example) In practical clinical pathogen detection applications, data from a single unknown sample is input into this random forest model for automated analysis: Input: Taking a blind test sample suspected of being infected with avian mumps virus type 1 in the test set as an example, the system automatically extracts its 35-dimensional feature vector matrix. Its core feature measured values ​​include: non-parametric quality feature identity_max = 99.3000, non-parametric counting feature Count_Correct_Contigs = 1.0000, average sequencing coverage depth Avg_Coverage_Depth_x_reads = 34.8000, and derived feature Align_Qlen_Ratio_Max = 0.3103, etc.

[0059] Model computation: The above 1×35 dimensional feature vector is input into the optimal random forest model with fixed hyperparameters. Multiple decision trees (n_estimators=50, max_depth=5) within the model are computed in parallel and voted on.

[0060] Output: The model outputs a probability value of 1.0000 for this sample to be the target pathogen (greater than the default threshold of 0.5). Final classification result: Label 1 (positive). This output result is completely consistent with the subsequent independent qPCR validation results, successfully achieving accurate identification of the pathogen.

[0061] 3. Full validation results of external independent test sets All 179 records from the external independent test set were batch-input into the optimal random forest model. The output full-validation confusion matrix and core evaluation metrics are shown in the appendix. Figure 6 , 7 As shown; Application validation shows that the random forest model built based on 35-dimensional features maintains extremely high pathogen discrimination stability and robustness when faced with novel clinical samples. This method completely eliminates the dependence of traditional single NGS analysis workflows on a large number of downstream qPCR validations, realizing fully automated intelligent identification from sequencing data to accurate pathogen identification, and possesses extremely high clinical application and industrial transformation value.

Claims

1. A method for virus identification in NGS sequences based on multidimensional feature fusion, characterized in that, Includes the following steps: S1. Parameterless assembly analysis and parameterless feature extraction: Through parameterless assembly and homology annotation, the original sequencing data is deeply cleaned, redundancy is removed and spliced, and the counting features and quality features of the sequences are extracted as parameterless features. S2. Construction of Sample-Specific Reference Sequence Library: Based on the parameterless assembly and homology annotation results of step S1, for suspected animal pathogenic viruses, the corresponding contigs are extracted and compared with the viral genome sequences in the nucleic acid sequence database. Based on the homology and query coverage indicators, the full-length viral genome or core gene sequence that is closest in evolutionary relationship and has the highest matching degree with the parameterless assembly sequence is selected from the database to construct the sample-specific reference sequence library. S3. Parametric Coverage Analysis and Feature Extraction: Based on the sample-specific reference sequence library constructed in S2, multi-dimensional feature extraction is performed to obtain parametric multi-dimensional features; the multi-dimensional features include one or more of the following: Reads dimension feature extraction, Contigs dimension feature extraction, or fusion dimension feature extraction. S4. Machine learning model design and accurate pathogen identification: Integrate non-parametric features and parametric multi-dimensional features, and calculate derived features based on quality features. Integrate non-parametric features, parametric multi-dimensional features and quality features to construct a multi-dimensional quantization matrix, and build a machine learning model through supervised learning. S5. Based on the machine learning model constructed in step S4, identify viruses in NGS sequences.

2. The virus identification method as described in claim 1, characterized in that, Step S1, which involves deep cleaning, redundancy removal, and splicing of the raw sequencing data, specifically includes: Data quality control and cleaning: To ensure the quality of downstream assembly, Trimmomatic data quality control software was used to perform quality control on the raw off-machine data; a sliding window filtering strategy was adopted to remove low-quality bases and linker sequences, with the sliding window size set to a range of 3 to 5 bases and the average quality cutoff value ranging from 15 to 25; short segments with a length less than a preset minimum length threshold, ranging from 20 bp to 50 bp, were also removed to obtain Cleanreads; Host sequence stripping: Clean reads were aligned to the host reference genome using Bowtie2 sequence alignment software, and unmapped reads were extracted using end-to-end alignment mode to eliminate host background interference; Low-abundance sequence capture and redundancy removal: Unmapped reads were assembled from scratch using MEGAHIT assembly software, employing a multiple k-mer iteration strategy with k-mer values ​​ranging from 21 to 119. Contigs with a length not less than a preset length threshold ranging from 50 bp to 200 bp were retained. Subsequently, the CD-HIT-EST tool was used to perform clustering and redundancy removal with a set sequence consistency threshold ranging from 0.90 to 0.99, constructing a globally non-redundant candidate set. The counting and quality features of the synchronously extracted sequences specifically include: Species annotation and feature quantification: BLASTx alignment was performed using the DIAMOND homology sequence retrieval tool. The upper limit of the expected value (E-value) was set to range from 1e-5 to 1e-15, the lower limit of identity was set to range from 30% to 70%, and the lower limit of alignment length was set to range from 30 bp to 100 bp. The species phylogeny was reconstructed based on the TaxonKit taxonomy tool, removing bacteriophages and plant and fungal viruses, and retaining candidate sequences of animal pathogenic viruses. During this process, the system groups samples by name and virus name, and automatically extracts the following parameter-free features: Counting features include one or more of the following: Total Contigs, Correct Contigs, or Noise Contigs. Quality features: For a subset of correctly annotated data, aggregate and calculate one or more of the following: maximum / mean homology (identity_max / identity_mean), maximum / mean / sum score (BitScore_max / BitScore_mean / BitScore_sum), maximum / sum alignment length (AlignLength_max / AlignLength_sum), average query sequence length (qlen_mean), minimum expected value (Evalue_min), or average mismatches (Mismatches_mean).

3. The virus identification method as described in claim 1, characterized in that, The sequence alignment described in step S2 is performed using NCBI's online BLASTn tool.

4. The virus identification method as described in claim 1, characterized in that, The Reads dimension feature extraction in step S3 specifically includes: aligning Clean reads to the reference library, parsing the output results, and extracting the initial aggregation parameters, including one or more of the following: Reads_Count_reads (alignment read length), Coverage_bp_reads (effective coverage bases), Coverage_Percent_reads (reference sequence coverage percentage), Avg_Coverage_Depth_x_reads (average sequencing coverage depth), Mean_BaseQ_reads (average base quality), Mean_MapQ_reads (average sequence alignment quality), or Histo_Max_Bin_Percent_reads (percentage of the maximum frequency interval in the coverage depth histogram). Contig dimensional feature extraction specifically includes: aligning non-redundant contigs to the same reference library, and similarly extracting the corresponding 7 parameters: total number of contig alignment reads (Reads_Count_contigs), total length of bases covered by the contig (Coverage_bp_contigs), percentage of genome coverage by the contig (Coverage_Percent_contigs), average sequencing depth of the contig (Avg_Coverage_Depth_x_contigs), average alignment quality value of the contig (Mean_MapQ_contigs), and the percentage of the largest interval in the contig coverage histogram (Histo_Max_Bin_Percent_contigs), to assess the macroscopic continuity of the assembled sequence; The fusion dimension feature extraction specifically includes: merging the BAM files of Reads and Contigs, calculating the comprehensive coverage under the superimposed state, and extracting the corresponding 7 feature parameters of the fusion sequence: including the total number of aligned reads (Reads_Count_fused), the total length of bases covered by the fusion sequence (Coverage_bp_fused), the coverage percentage of the fusion sequence (Coverage_Percent_fused), the average coverage depth of the fusion sequence (Avg_Coverage_Depth_x_fused), the average base quality value of the fusion sequence (Mean_BaseQ_fused), the average alignment quality value of the fusion sequence (Mean_MapQ_fused), and the maximum interval percentage of the coverage histogram of the fusion sequence (Histo_Max_Bin_Percent_fused), reflecting the comprehensive mapping strength of the global dimension.

5. The virus identification method as described in claim 1, characterized in that, The calculation formula for derived features based on quality features in step S4 is as follows: ； AlignLength_max=max{AlignLength_i|Right_Annotated_Contigs_i =1} qlen_mean=(1 / n)×Σ(qlen_i) for i where Right_Annotated_Contigs_i=1 in, Align_Qlen_Ratio_Max: Optimal coverage ratio; AlignLength: The maximum alignment length among all correctly annotated contigs within the group; i iterates through all correct comments (Contig) within the group; qlen: The arithmetic mean of the query sequence lengths of all correctly annotated contigs within the group; n is the total number of correctly annotated contigs in this group.

6. The virus identification method as described in claim 1, characterized in that, Step S4, which involves calculating and integrating derived features from non-parametric features, parametric multi-dimensional features, and quality features to construct a multi-dimensional quantization matrix, and then building a machine learning model through supervised learning, specifically includes: Data preprocessing and labeling: qPCR validation results were introduced as the gold standard, and binary labels were assigned: positive = 1, negative = 0. The feature matrix was then subjected to outlier removal, median imputation, and standard deviation normalization, and then the training and test sets were divided. Model training and evaluation: Random forest (RF), support vector machine (SVM), k-nearest neighbors (KNN), logistic regression (LR), decision tree (DT), and gradient boosting tree (GBDT) were used to construct classification models in parallel. Grid search combined with 5-fold cross-validation was used to optimize hyperparameters.

7. The virus identification method as described in claim 1, characterized in that, In step S5, the machine learning model constructed in step S4 is used to identify viruses in NGS sequences. The input is the virus sample to be detected, and the output is the target pathogen of the sample to be detected.

8. The virus identification method as described in claim 1, characterized in that, The virus identification method is used for identifying low-abundance viruses in complex waterfowl samples.

Images