Abnormity identification traceability method based on coronavirus high-throughput detection

Abstract

The invention relates to the technical field of virus traceability, and particularly discloses an anomaly identification traceability method based on coronavirus high-throughput detection, comprising the following steps: S1, acquiring high-throughput original sequencing data of a target sample; s2, inputting the original sequencing data into a dynamic learning type recognition model, and outputting each virus pedigree and the credibility of each pedigree; s3, performing noise perception variation detection on the original sequencing data; s4, aiming at each virus lineage, determining a variation point; s5, constructing a Bayesian network to determine a genetic relationship among the variation sites, and generating a virus haplotype sequence based on the genetic relationship; s6, calculating a genetic distance between the abnormal variation site and the reference database, and deducing a propagation path through a maximum likelihood method; and outputting a traceability report, wherein the traceability report comprises the geographic position and the propagation timeline of the abnormal variation source. A dynamic learning type recognition model is adopted, and continuous optimization of the recognition capability of new variants is ensured through incremental learning mechanism training.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of virus tracing, and particularly relates to an abnormal identification and tracing method based on coronavirus high-throughput detection. BACKGROUND

[0002] With the development of neural network technology, in the coronavirus high-throughput detection, a preset identification model (such as BERT) is gradually used to determine the virus lineage, but the following problems still exist at present.

[0003] 1. The training data may not cover emerging or rare virus variants, there are model bias or data imbalance problems, for example, if the training data is focused on known lineages, the identification accuracy of unknown variant viruses (such as newly emerging viruses) may be significantly reduced, in the outbreak scenario, the model may misjudge the lineage, leading to subsequent analysis errors.

[0004] 2. The variant site detection is based on the comparison of the preset virus reference genome, but the reference genome may be incomplete or outdated, in high-throughput sequencing, the coverage depth of certain genome regions is insufficient, which may lead to missed detection of variants. Moreover, the linkage relationship is calculated based on the co-occurrence frequency (such as the distribution of variant site sets in segmented data), but high-throughput data often contains sequencing noise, if the error correction algorithm is not integrated, only the simple frequency threshold is relied on, false positive linkage (such as non-genetically associated random co-occurrence) may be caused.

[0005] 3. When matching the variant site, it is generally assumed that the expected frequency of each virus lineage can be preset or predicted, but in practice, the frequency distribution is complex when mixed infection or multiple lineages coexist, if the measured frequency of the candidate variant site is different from the expected value by more than a threshold value, the site is discarded, which may cause important low-frequency variants to be ignored (such as drug resistance related mutations), especially in the dynamic evolution scenario of virus quasispecies.

[0006] 4. The linkage relationship is only inferred by the co-occurrence frequency (such as the co-occurrence of variant site sets in segmented data), but the recombination in virus genetics is often ignored, which may cause false association of non-linkage sites. The scheme relies on high-frequency variant site screening haplotype, but low-frequency linkage (such as <5% frequency) may reflect the key evolutionary path, but is excluded due to abundance filtering. SUMMARY

[0007] In order to solve the problems in the prior art, the present application provides an abnormal identification and tracing method based on coronavirus high-throughput detection.

[0008] The technical scheme adopted by the present application is:

[0009] An abnormal identification and tracing method based on coronavirus high-throughput detection, comprising the following steps:

[0010] S1, obtaining high-throughput raw sequencing data of a target sample, the raw sequencing data being composed of segment data, each segment data representing a short sequence fragment in a nucleic acid sequence of the target sample;

[0011] S2, inputting the raw sequencing data into a dynamic learning type recognition model, the model being trained through an incremental learning mechanism, and outputting each virus lineage of each virus contained in the target sample and a credibility of each virus lineage;

[0012] S3, performing noise-aware variant detection on the raw sequencing data to determine each variant site contained in the raw sequencing data and a variant frequency of each variant site, and to determine a linkage relationship between the variant sites, the linkage relationship being used to reflect a genetic mode of whether part of the variant sites are transmitted as a whole to offspring;

[0013] S4, for each virus lineage, determining a variant site matched with a preset specific variant site corresponding to the virus lineage from the variant sites as a variant site belonging to the virus lineage;

[0014] S5, constructing a Bayesian network to determine a genetic relationship between the variant sites according to the virus lineage to which each variant site belongs, the variant frequency of each variant site, and the linkage relationship between the variant sites, and generating a virus haplotype sequence based on the genetic relationship;

[0015] S6, marking a variant site satisfying a threshold of a deviation of a variant frequency from a historical distribution of a global virus database or a preset risk threshold of a confidence of a linkage relationship as an abnormal variant site; calculating a genetic distance between the abnormal variant site and a reference database, inferring a transmission path through a maximum likelihood method; and outputting a tracing report, the tracing report including a geographical location and a transmission timeline of an abnormal variant source.

[0016] Preferably, in step S1, the high-throughput raw sequencing data is from an Illumina Novaseq 6000 or MGI DNBSEQ-T7 sequencing platform, and the sequencing mode is PE150.

[0017] Preferably, in step S2, the incremental learning mechanism triggers fine-tuning based on a data distribution offset degree of the raw sequencing data, the offset degree being calculated by KL divergence according to a data distribution of the raw sequencing data and a data distribution of training sequencing data, and when the offset degree exceeds an offset degree threshold, global abnormal virus variant data is injected to update model parameters.

[0018] Preferably, in step S3, the noise-aware variant detection uses a convolutional autoencoder to denoise the segment data, a loss function is composed of a reconstruction error and a sparsity constraint, and a hidden Markov model interpolation is applied to a low-coverage region with a depth less than a threshold.

[0019] Preferably, in step S3, the loss function L is:

[0020] L = ‖X noisy -X clean ‖2+λ·KL(q(z|X noisy )‖p(z))

[0021] wherein X noisy is the original sequencing data aligned segment data sequence, X clean is the denoised segment data sequence, ‖X noisy -X clean ‖2is the Euclidean distance between X noisy and X clean , q(z|X noisy ) is the posterior distribution of the input segment data sequence in the latent variable space, p(z) is the preset latent variable prior distribution, and λ is the weight of the translation reconstruction error and distribution constraint.

[0022] Preferably, in step S3, the interpolation process of the hidden Markov model integrates the recombination rate correction, and the state transition probability formula is:

[0023] P adj (s t |s t-1 )=P(s t |s t-1 )×(1+α·r recomb )

[0024] wherein P adj (s t |s t-1 ) is the corrected transition probability of site s t-1 to s t , P(s t |s t-1 ) is the standard state transition probability of site s t-1 to s t , r recomb is the preset virus recombination rate, and α is the virus type coefficient.

[0025] Preferably, in step S5, in the construction of the Bayesian network, the edge weight W ij is calculated according to the following formula:

[0026]

[0027] wherein d ij is the genomic distance of sites i and j, β is 2.0, and is used to enhance the detection sensitivity of abnormal linkage.

[0028] Preferably, in step S6, the genetic distance calculation adopts the Tamura-Nei model, the transmission path inference combines the spatiotemporal metadata, and the abnormal threshold is set as genetic distance > 0.05 and confidence < 75%.

[0029] The beneficial effects of the present application are:

[0030] 1. The dynamic learning type recognition model is adopted, so that the system can quickly adapt to emerging viruses, and the incremental learning mechanism is trained, without the need to retrain the entire network, significantly shortening the response period and ensuring the continuous optimization of the recognition ability of new variants.

[0031] 2. The noise perception variation detection technology effectively solves the problem of sensitivity to sequencing noise, significantly reduces the interference of sequencing errors, and improves the accuracy of variation site detection;

[0032] 3. The dynamic update retains low-frequency but important variation sites, ensuring that critical signals are not filtered out;

[0033] 4. The recombination-aware Bayesian network is constructed, which integrates the virus recombination rate parameter and the genomic distance decay coefficient, accurately quantifies the genetic linkage strength between variation sites, effectively distinguishes real genetic association from random noise, and significantly reduces false positive results;

[0034] 5. By calculating the genetic distance of the virus sequence and the global database, combined with the spatiotemporal metadata of the sample (such as geographical location and sampling time), a transmission path report can be generated. BRIEF DESCRIPTION OF DRAWINGS

[0035] Figure 1 The method flowchart in the embodiments of the present application is shown. DETAILED DESCRIPTION

[0036] The present application will be described in detail below in conjunction with the drawings and embodiments.

[0037] EMBODIMENT

[0038] S1, obtaining high-throughput raw sequencing data of a target sample, the raw sequencing data being composed of segmented data, each segmented data representing a short sequence fragment in the nucleic acid sequence of the target sample.

[0039] The high-throughput raw sequencing data comes from the Illumina Novaseq 6000 or MGI DNBSEQ-T7 sequencing platform, and the sequencing mode is PE150.

[0040] S2, inputting the raw sequencing data into a dynamic learning type recognition model, the model being trained by an incremental learning mechanism, and outputting the virus spectrum of each virus belonging to the target sample and the credibility of each spectrum.

[0041] The incremental learning mechanism triggers fine-tuning based on a data distribution offset of the original sequencing data, the offset being calculated by KL divergence D between a data distribution of the original sequencing data and a data distribution of the training sequencing data KL (P new ||P old ), wherein P new is the data distribution of the original sequencing data, and P old is the data distribution of the training sequencing data, when the offset D KL > 0.3, inject global abnormal virus variant data (such as GISAID database) to update the model parameters.

[0042] S3, noise-aware variant detection is performed on the original sequencing data to determine each variant site contained in the original sequencing data and the variant frequency of each variant site, and to determine the linkage relationship between each variant site, the linkage relationship being used to reflect whether part of the variant sites in the variant sites are transmitted as a whole to the genetic mode of the offspring.

[0043] The noise-aware variant detection uses a convolutional autoencoder to denoise the segmented data, and a loss function is composed of a reconstruction error and a sparsity constraint, and a hidden Markov model interpolation is applied to a low coverage region with a depth less than a threshold.

[0044] The loss function L is:

[0045] L = ‖X noisy -X clean ‖2+λ·KL(q(z|X noisy )‖p(z))

[0046] wherein X noisy is the segmented data sequence of the original sequencing data after alignment, X clean is the segmented data sequence after denoising, ‖X noisy -X clean ‖2 is the Euclidean distance between X noisy and X clean , q(z|X noisy ) is the posterior distribution of the input segmented data sequence in the hidden variable space, p(z) is a preset prior distribution of the hidden variable, and λ is the weight of the translation reconstruction error and the distribution constraint.

[0047] The interpolation process of the hidden Markov model integrates the recombination rate correction, and the state transition probability formula is:

[0048] P adj (s t |s t-1 )=P(s t |s t-1 )×(1+α·r recomb )

[0049] where P(s adj (s t |s t-1 ) is the corrected transition probability from site s t-1 to s t , P(s t |s t-1 ) is the standard state transition probability from site s t-1 to s t , r recomb is the preset virus recombination rate, and a is a virus type coefficient (0.8 for RNA viruses).

[0050] S4. For each virus lineage, determine the variation sites that match the preset specific variation sites corresponding to the virus lineage from the variation sites, as the variation sites belonging to the virus lineage.

[0051] The specific variation sites are obtained from a dynamically updated global abnormal site library, and the site library is synchronized once every 48 hours to include emerging abnormal variation data.

[0052] S5. According to the virus lineage to which each variation site belongs, the variation frequency of each variation site, and the linkage relationship between the variation sites, a Bayesian network is constructed to determine the genetic relationship between the variation sites, and a virus haplotype sequence is generated based on the genetic relationship.

[0053] In the construction of the Bayesian network, the edge weight W ij is calculated by the following formula:

[0054]

[0055] where d ij is the genomic distance between sites i and j, and b is 2.0, which is used to enhance the detection sensitivity of abnormal linkage.

[0056] S6. Mark the variation sites that meet the threshold of variation frequency deviating from the historical distribution of the global virus database or the confidence of linkage relationship being lower than the preset risk threshold as abnormal variation sites; calculate the genetic distance of the abnormal variation sites from the reference database, and infer the transmission path by maximum likelihood method; output the tracing report, which includes the geographical location and transmission timeline of the abnormal variation source.

[0057] The genetic distance is calculated using the Tamura-Nei model, the transmission path is inferred combined with spatiotemporal metadata, and the abnormal threshold is set to genetic distance > 0.05 and confidence < 75%.

[0058] The above method has the following advantages for abnormal identification and tracing of coronavirus high-throughput detection:

[0059] 1. Dynamic improvement of the model

[0060] Introduce dynamic learning-based identification model, automatically trigger incremental learning mechanism by real-time monitoring of data distribution shift (such as KL divergence > 0.3). This mechanism only fine-tunes the parameters of the fully connected layer, while preserving historical knowledge and injecting global new variant data (such as GISAID database XBB.1.5 sequence). This design enables the system to adapt to viral evolution, such as when unknown mutations appear in sewage samples, the model can quickly update and improve the identification accuracy of emerging variants, avoiding the waste of computing resources caused by traditional retraining.

[0061] 2. Data reliability enhancement

[0062] Design a double-track error correction mechanism:

[0063] Convolutional autoencoder: learn noise distribution features through encoding-decoding structure, separate real signals and sequencing errors (such as Illumina platform base substitution errors) in hidden variable space.

[0064] Recombination correction HMM: For the recombination characteristics of viral genomes, integrate the recombination rate parameter (r recomb ) in the hidden Markov model to optimize the state transition probability of low coverage areas.

[0065] 3. Low-frequency signal integrity protection

[0066] Build a dynamic abnormal site library, synchronize global emerging variation data every 48 hours (such as Spike: K444T drug-resistant site), and set intelligent retention rules:

[0067] Variations in functionally conserved regions (such as polymerase coding regions) are marked as key sites even if the frequency is below the threshold;

[0068] Low-frequency variations with strong linkage to known drug resistance / virulence sites (confidence > 85%) are included in the analysis.

[0069] 4. Genetic analysis depth optimization

[0070] Integrate two types of key parameters in Bayesian networks:

[0071] Recombination rate r recomb : Quantify the correlation strength of sites based on historical recombination hotspot data (such as coronavirus S protein region);

[0072] Distance attenuation coefficient β: weaken the random co-occurrence interference of distant sites through genomic distance d ij .

[0073] The design reduces the false positive rate of genetic linkage and accurately locates the cross-lineage co-evolution module (such as Spike-RdRp linkage block) in the analysis of Delta-Omicron recombinant strains.

[0074] Finally, the genetic distance is calculated by the Tamura-Nei model, and the propagation network is constructed by combining the spatiotemporal metadata (sampling place / time) of the samples, so as to further generate a traceability report containing a propagation path diagram and a timeline.

[0075] The above embodiments only express specific implementations of the present application, which are described in detail and specifically, but cannot be understood as a limitation on the scope of the patent of the present application. It should be noted that for ordinary skilled persons in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which are within the protection scope of the present application.

Claims

1. An anomaly identification and tracing method based on high-throughput detection of coronaviruses, characterized in that, The method comprises the following steps: S1, obtaining high-throughput raw sequencing data of a target sample, the raw sequencing data being composed of segment data, each segment data representing a short sequence fragment in a nucleic acid sequence of the target sample; S2, inputting the raw sequencing data into a dynamic learning type identification model, the model being trained through an incremental learning mechanism, and outputting each virus spectrum of viruses contained in the target sample and a credibility of each spectrum; S3, performing noise-aware variant detection on the raw sequencing data to determine each variant site contained in the raw sequencing data and a variant frequency of each variant site, and determine a linkage relationship between the variant sites, the linkage relationship being used to reflect a genetic mode of whether part of the variant sites are transmitted as a whole to offspring; S4, for each virus spectrum, determining a variant site matched with a preset specific variant site corresponding to the virus spectrum from the variant sites as a variant site belonging to the virus spectrum; S5, constructing a Bayesian network according to the virus spectrum to which each variant site belongs, the variant frequency of each variant site and the linkage relationship between the variant sites to determine a genetic relationship between the variant sites, and generating a virus haplotype sequence based on the genetic relationship; S6, marking a variant site meeting a threshold of a deviation of a variant frequency from a historical distribution of a global virus database or a credibility of a linkage relationship being lower than a preset risk threshold as an abnormal variant site, calculating a genetic distance of the abnormal variant site from a reference database, inferring a transmission path through a maximum likelihood method, and outputting a tracing report, the tracing report including a geographical position of an abnormal variant source and a transmission timeline.

2. The method of claim 1, wherein the method further comprises: In step S1, the high-throughput raw sequencing data is from an Illumina Novaseq 6000 or MGI DNBSEQ-T7 sequencing platform, and the sequencing mode is PE150.

3. The method of claim 2, wherein the method further comprises: In step S2, the incremental learning mechanism triggers fine-tuning based on a data distribution offset degree of the raw sequencing data, the offset degree being calculated by KL divergence according to a data distribution of the raw sequencing data and a data distribution of training sequencing data, and when the offset degree exceeds an offset degree threshold, global abnormal virus variant data is injected to update model parameters.

4. The method of claim 3, wherein the method further comprises: In step S3, the noise-aware variant detection uses a convolutional autoencoder to denoise the segment data, and a loss function is composed of a reconstruction error and a sparsity constraint, and a hidden Markov model interpolation is applied to a low-coverage region with a depth less than a threshold.

5. The method of claim 4, wherein the method further comprises: In step S3, the loss function L is: L = ‖X noisy - X clean ‖2+λ·KL(q(z|X noisy )‖p(z)) where X noisy is the segmented data sequence after alignment of the original sequencing data, X clean is the segmented data sequence after denoising, ‖X noisy -X clean ‖2 is the Euclidean distance between X noisy and X clean , q(z|X noisy ) is the posterior distribution of the input segmented data sequence in the latent variable space, p(z) is the preset latent variable prior distribution, and λ is the weight of the translation reconstruction error and distribution constraint.

6. The method of claim 5, wherein the method further comprises: In step S3, the interpolation process of the hidden Markov model integrates a recombination rate correction, and a state transition probability formula is: P adj (s t |s t-1 )=P(s t |s t-1 )×(1+α·r recomb ) Among them, P adj (s t |s t-1 ) is the site s t-1 To s t The corrected transition probability, P(s) t |s t-1 ) is the site s t-1 To s t The standard state transition probability, r recomb α is the preset viral recombination rate, and α is the viral type coefficient.

7. The method of claim 6, wherein the method further comprises: In step S5, in the construction of the Bayesian network, the edge weight W ij The calculation formula is: where d ij is the genomic distance between sites i and j, and β is taken as 2.0 to enhance the sensitivity of detecting aberrant linkage.

8. The method of claim 7, wherein the method further comprises: In step S6, the genetic distance calculation adopts a Tamura-Nei model, the transmission path inference combines spatiotemporal metadata, and the abnormal threshold is set as genetic distance > 0.05 and credibility < 75%.

Images