A pathogen gene sequence detection and analysis method based on nanopore sequencing

Abstract

The application discloses a pathogen gene sequence detection and analysis method based on nanopore sequencing, and relates to the field of customs quarantine.The method comprises the following steps: nanopore sequencing is performed on collected samples, original sequencing data is acquired, quality screening is performed, and high-quality sequences are obtained; species classification and mutation recognition are completed; pathogen types and mutation characteristics are extracted; a transmission map is constructed based on classification information, mutation information, sampling time and location, and a dynamic transmission structure diagram is formed; a graph diffusion modeling calculation is performed to calculate the transmission path score between samples; a Bayesian causal model is constructed by combining the mutation characteristics and the transmission score, and the posterior probability of a pathogen causing a high-risk event is inferred; a risk level label is generated according to the posterior probability, a pathogen analysis report and early warning information are output, and integrated analysis of rapid detection and risk assessment is realized.

Description

TECHNICAL FIELD

[0001] The present specification relates to the field of customs quarantine, and more particularly, a pathogen gene sequence detection and analysis method based on nanopore sequencing. BACKGROUND

[0002] Currently, in the scenarios of customs quarantine, port prevention and control, and entry pathogen monitoring, traditional pathogen detection and analysis methods mainly rely on fluorescence PCR, specific primer amplification, enzyme-linked immunosorbent assay, or high-throughput sequencing technology. These methods have the following limitations:

[0003] 1. Long detection period, difficult to deploy on site: High-throughput sequencing relies on laboratory environment and centralized sequencing platform, which cannot meet the requirements of real-time detection and rapid response in mobile scenarios such as ports and high-speed channels.

[0004] 2. Only limited target recognition, unable to monitor the full spectrum: Traditional PCR and antigen detection require pre-set detection targets, which cannot achieve the discovery of unknown pathogens and parallel monitoring of multiple types of mutations, and have limited recognition ability for variant strains.

[0005] 3. Lack of transmission map modeling mechanism: Most existing solutions focus on pathogen detection or genome annotation, lack of structured modeling of potential transmission relationships between samples, and are difficult to restore their transmission path or relevance.

[0006] 4. Risk identification method is relatively simple and cannot evaluate the transmission potential: Most current evaluation methods are based on pathogen species or specific mutation markers to evaluate risks, and fail to systematically model the transmission mechanism, transmission location, and mutation synergistic effects.

[0007] Therefore, there is an urgent need for a pathogen detection and analysis solution that supports rapid detection, structural mapping, intelligent modeling, and risk identification. SUMMARY

[0008] A series of simplified concepts are introduced in the summary section, which will be further described in detail in the specific embodiments section. The summary section of the present application does not mean to attempt to limit the key features and necessary technical features of the claimed technical solution, nor to attempt to determine the protection scope of the claimed technical solution.

[0009] In a first aspect, the present application proposes a pathogen gene sequence detection and analysis method based on nanopore sequencing, comprising:

[0010] Performing nanopore sequencing on the collected sample to obtain raw sequencing data;

[0011] Performing quality screening on the above raw sequencing data to obtain screened sequencing data;

[0012] performing a species classification operation and a mutation point monitoring operation on the sequencing data after the screening to obtain pathogen classification information and mutation information;

[0013] constructing a transmission graph based on the pathogen classification information and the mutation information to obtain a dynamic transmission structure graph of the pathogen;

[0014] performing a graph diffusion modeling operation on the dynamic transmission structure graph to obtain potential transmission path score information of the pathogen;

[0015] constructing a Bayesian causal model based on the mutation information and the potential transmission path score information to obtain a posterior probability of the pathogen causing a high-risk event;

[0016] generating a risk level label according to the posterior probability, and outputting a pathogen analysis report and early warning information.

[0017] In a feasible implementation, the constructing a transmission graph based on the pathogen classification information and the mutation information to obtain a dynamic transmission structure graph of the pathogen includes:

[0018] extracting metadata information of each sample, wherein the metadata information includes the pathogen classification information, the mutation information, sample collection time, sample collection location, and a sequencing timestamp;

[0019] grouping samples with the same or similar species into the same transmission graph construction unit according to the pathogen classification information;

[0020] generating a mutation feature vector for each sample based on the mutation information, wherein the mutation feature vector is encoded by a whole genome mutation site;

[0021] representing each sample as a graph node, wherein the graph node includes the pathogen classification information, the mutation information, the sample collection time, the sample collection location, and the sequencing timestamp;

[0022] calculating a genomic similarity score between sample nodes based on distance measurement between the mutation feature vectors;

[0023] establishing a transmission edge connection between the graph nodes according to the similarity score and sampling time and geographical location information, under the condition that at least one preset condition is met, wherein the preset condition includes that the similarity score is higher than a set threshold, the sample collection time interval is less than a preset window, and the sample collection location is in a potential transmission path;

[0024] setting an edge weight value in the transmission edge connection, wherein the edge weight value is calculated from the similarity score, the time difference, and the geographical distance;

[0025] The dynamic propagation structure graph is established according to the graph node, the propagation edge connection and the edge weight through a time sliding window mechanism.

[0026] In an available implementation, the graph diffusion modeling operation is performed on the dynamic propagation structure graph to obtain the potential propagation path score information of the pathogen, including:

[0027] The graph node in the dynamic propagation structure graph is encoded to construct a graph node feature matrix.

[0028] A time series modeling mechanism is introduced into the dynamic propagation structure graph, and the graph structure is time-sequenced based on sample collection time to form a time-aware dynamic propagation graph.

[0029] A graph diffusion model is used to perform graph diffusion modeling according to the time-aware dynamic propagation graph and the graph node feature matrix to obtain a graph node propagation embedding vector.

[0030] The propagation potential is extracted according to the graph node propagation embedding vector and the time-aware dynamic propagation graph to obtain a node potential propagation vector.

[0031] The similarity information between each two node potential propagation vectors of the propagation path is calculated to obtain the potential propagation path score information.

[0032] In an available implementation, a Bayesian causal model is constructed based on the mutation information and the potential propagation path score information to obtain the posterior probability of the pathogen causing a high-risk event, including:

[0033] A high-risk mutation feature identification operation is performed on the mutation feature vector to construct a first causal variable of the Bayesian causal model.

[0034] A propagation potential identification operation is performed on the potential propagation path score information to construct a second causal variable of the Bayesian causal model.

[0035] The Bayesian causal graph model is constructed based on the first causal variable and the second causal variable, wherein the Bayesian causal model includes a posterior variable node representing a high-risk event and a directed edge from the causal variable node to the posterior variable node.

[0036] According to the prior conditional probability and the joint distribution inference rule, the conditional probability of the posterior variable node in the Bayesian causal graph model is calculated to obtain the posterior probability of the pathogen causing a high-risk event.

[0037] In an implementable embodiment, the operation of identifying the high-risk mutation feature of the mutation feature vector to construct the first causal variable of the Bayesian causal model comprises:

[0038] obtaining a high-risk mutation site reference library;

[0039] comparing the mutation feature vector and the high-risk mutation site reference library to identify a high-risk mutation related site to obtain a comparison result;

[0040] obtaining statistical information of the number of mutation feature matches according to the comparison result and a preset mutation risk determination rule;

[0041] determining the first causal variable of the Bayesian causal model according to the statistical information.

[0042] In an implementable embodiment, the determination of the preset mutation risk determination rule comprises:

[0043] constructing a sliding time window sample set based on a mutation risk score, wherein the sliding time window sample set is all pathogen samples received within a preset time window, and the mutation risk score is calculated from the mutation feature vector;

[0044] performing distribution statistics on the mutation risk scores in the sliding time window sample set, and calculating a high quantile value of the distribution as a dynamic threshold;

[0045] comparing the mutation risk score of the current sample with the dynamic threshold to determine whether to update the preset mutation risk determination rule.

[0046] In an implementable embodiment, the operation of identifying the propagation potential of the potential transmission path score information to construct the second causal variable of the Bayesian causal model comprises:

[0047] determining propagation potential determination threshold information;

[0048] obtaining transmission path number information, maximum transmission probability information, and structural centrality indicators according to the propagation potential determination threshold information and the potential transmission path score information;

[0049] determining the second causal variable of the Bayesian causal model according to the transmission path number information, the maximum transmission probability information, the structural centrality indicators, and a preset propagation potential determination rule.

[0050] In an implementable embodiment, the determination of the second causal variable of the Bayesian causal model according to the transmission path number information, the maximum transmission probability information, the structural centrality indicators, and the preset propagation potential determination rule comprises:

[0051] determine the propagation probability driver, the propagation depth adjustment term, the target node importance weighting term, and the spatial reachability adjustment term according to the above-mentioned propagation path quantity information, the above-mentioned maximum propagation probability information, and the above-mentioned structural centrality index;

[0052] determine the propagation potential score function according to the above-mentioned propagation probability driver, the above-mentioned propagation depth adjustment term, the above-mentioned target node importance weighting term, and the above-mentioned spatial reachability adjustment term;

[0053] determine the second causal variable of the Bayesian causal model according to the above-mentioned propagation potential score function and the above-mentioned preset propagation potential determination rule.

[0054] In a feasible implementation, the generating of the risk level label according to the posterior probability, and the output of the pathogen analysis report and the early warning information, comprises:

[0055] obtaining dynamic level threshold information;

[0056] determining the risk level label according to the posterior probability and the dynamic level threshold information;

[0057] outputting the pathogen analysis report and the early warning information according to the risk level label and the sample analysis result.

[0058] In a feasible implementation, the specific determination of the dynamic level threshold information comprises:

[0059] constructing a historical pathogen sample set containing the posterior probability value and the corresponding true risk label;

[0060] constructing a receiver operating characteristic curve of the historical sample set with the posterior probability as the prediction variable and the risk label as the true classification standard;

[0061] calculating the sensitivity and specificity indexes corresponding to different posterior probability thresholds, and forming a Youden index based on the sum of the sensitivity and the specificity minus one;

[0062] taking the posterior probability value with the maximum Youden index as the dynamic judgment threshold of the first risk level division;

[0063] performing multi-threshold optimization on the posterior probability value in combination with the macro-average AUC index and the classification accuracy evaluation index, to form the probability interval corresponding to the high, medium, medium-high, and low risk levels.

[0064] In summary, the method proposed in the embodiments of the present application adopts nanopore sequencing technology, has the characteristics of plug-and-play, read length support, and on-the-fly analysis, supports on-site deployment at border ports, and can complete pathogen whole genome sequencing and analysis within a few hours. Through species classification and mutation monitoring algorithms based on reference library or depth model, multiple species, multiple subtypes and their mutation combinations can be identified at the same time, realizing rapid tracing of imported pathogens. The sequencing sample is modeled as a graph node, and a transmission graph is constructed based on mutation similarity, sampling space-time, and transmission reachability, which can reconstruct the potential transmission relationship of pathogens between different individuals / regions, support transmission path level modeling and visualization. Through graph diffusion algorithm, the transmission graph is modeled to calculate the transmission path score matrix between sample nodes, realize the joint evaluation of transmission potential from structure and behavior, and solve the problem of weak structure expression ability of traditional scoring model. The pathogen mutation characteristics and transmission potential are introduced as causal variables, the posterior probability of causing high-risk events is calculated through Bayesian network, an interpretable risk reasoning path is constructed, and quantitative modeling and reasoning identification of potential risk samples are supported. Combined with posterior probability and dynamic decision threshold, risk level labels are generated, and structured reports are output for access by border inspection, flow analysis and other systems, forming an early warning integrated mechanism from detection to disposal.

[0065] Other advantages, objects, and features of the present application will be apparent to those skilled in the art from the following detailed description of the preferred embodiments. BRIEF DESCRIPTION OF DRAWINGS

[0066] Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The accompanying drawings are intended to illustrate preferred embodiments of the present application, and should not be considered limiting of the present application. Indeed, the drawings can illustrate only one preferred embodiment of the present application and therefore should not be considered to narrow the scope of the present application. Like reference numerals are used to refer to like parts throughout the several views of the drawings.

[0067] Figure 1 A flowchart of a pathogen gene sequence detection and analysis method based on nanopore sequencing provided by the embodiments of the present application. DETAILED DESCRIPTION

[0068] The terms "first", "second", "third", "fourth" etc. (if any) in the description and claims of the present application and the above drawings are used to distinguish similar objects, and do not have to be necessarily used to describe a specific order or sequence. It should be understood that the data thus used can be interchanged under appropriate circumstances, so that the embodiments described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device that includes a series of steps or units does not have to be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to these processes, methods, products or devices. The technical solutions in the embodiments of the present application will be described in detail below in conjunction with the drawings of the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all.

[0069] Please refer to Figure 1 A flowchart of a pathogen gene sequence detection and analysis method based on nanopore sequencing provided by the embodiments of the present application is shown, which can specifically include:

[0070] S110, nanopore sequencing is performed on the collected sample to obtain raw sequencing data;

[0071] Exemplarily, suspicious samples are collected from inbound personnel, the surface of goods, transportation carriers, logistics chains or cold chain transportation links, and nanopore sequencing technology is used for real-time sequencing to obtain raw signal data containing complete pathogen genome sequences. Nanopore sequencing has the advantages of long read length, online analysis, and suitability for rapid detection in the field, and is particularly suitable for deployment in scenarios such as ports, expressways and airports.

[0072] S120, quality screening is performed on the above raw sequencing data to obtain screened sequencing data;

[0073] Exemplarily, to improve the accuracy of downstream analysis, the raw sequencing data needs to be processed, and the quality screening can include but is not limited to removing low-quality reads, filtering out short sequences, correcting sequencing errors (such as base calling errors), and selecting high-quality fragments (usually evaluated by Q value or Phred value). The sequencing data obtained after screening should have high confidence and coverage, and be suitable for subsequent alignment, mutation identification and other operations.

[0074] S130, species classification operation and mutation point monitoring operation are performed on the above screened sequencing data to obtain pathogen classification information and mutation information;

[0075] Exemplarily, a species identification is performed on the sequencing fragments by using a comparison-based classification algorithm (such as Minimap2+Kraken2, etc.) or a deep learning model, and corresponding pathogen classification information (genus, species, subtype, etc.) is extracted. Meanwhile, base differences in the sequence from a reference genome are detected, and mutation points, including SNPs, InDels, and variant segments, are recorded to construct a mutation feature vector, thereby providing a basis for subsequent graph modeling and risk analysis.

[0076] S140, a transmission map is constructed based on the pathogen classification information and the mutation information, so as to obtain a dynamic transmission structure diagram of the pathogen.

[0077] Exemplarily, metadata information (species category, mutation information, sampling time, sampling location, etc.) of each sample is first extracted, each sample is represented as a graph node, gene similarity between samples is calculated, and time / space reachability conditions are introduced to establish a transmission edge connection. A sliding time window mechanism is used to update the graph structure at each time period, and a dynamic transmission map G t is constructed, which depicts a pathogen transmission network.

[0078] S150, graph diffusion modeling is performed on the dynamic transmission structure diagram, so as to obtain potential transmission path score information of the pathogen.

[0079] Exemplarily, the transmission graph G t and the node feature matrix X t are input into a graph diffusion modeling algorithm, and a propagation potential learning is performed through a graph neural network structure to generate a propagation embedding vector of each sample node. Further, a potential transmission path score between any two nodes is calculated through a concatenation score function or a similarity function to form a transmission path score matrix, which is used for subsequent causal analysis.

[0080] S160, a Bayesian causal model is constructed based on the mutation information and the potential transmission path score information, so as to obtain a posterior probability of a high-risk event triggered by the pathogen.

[0081] Exemplarily, the mutation risk feature and the propagation potential score are jointly used as input variables of a causal graph, and a condition probability inference is performed through a Bayesian network with clear structure. The model output is a posterior probability of a sample triggering a high-risk event (such as large-scale transmission and high-pathogenicity outbreak), which reflects a potential threat of a pathogen represented by the sample to public safety.

[0082] S170, a risk level label is generated according to the posterior probability, and a pathogen analysis report and early warning information are output.

[0083] Exemplarily, in combination with the dynamic level determination threshold, the posterior probability is mapped to a risk level label (such as high, medium, and low). Subsequently, the sample identification, mutation information, transmission score, and risk inference result are integrated to output a structured pathogen analysis report for automatic collection and early warning use by the border inspection, flow survey, and epidemic monitoring departments.

[0084] To sum up, the method proposed in the embodiments of the present application adopts nanopore sequencing technology, has the characteristics of plug-and-play, read length support, and on-the-fly analysis, supports the deployment at the border port or on the spot, and can complete the whole genome sequencing and analysis of the pathogen within a few hours. Through the species classification and mutation monitoring algorithm based on the reference library or the depth model, multiple species, multiple subtypes, and their mutation combinations can be identified at the same time, and the rapid tracing of the imported pathogen is realized. The sequencing sample is modeled as a graph node, and a transmission graph is constructed based on the mutation similarity, sampling space-time, and transmission reachability, which can reconstruct the potential transmission relationship of the pathogen between different individuals / regions, support the modeling and visualization of the transmission path level. Through the graph diffusion algorithm, the transmission graph is modeled, the transmission path score matrix between sample nodes is calculated, the transmission potential is evaluated from the structure and behavior, and the problem of weak structure expression capability of the traditional scoring model is solved. The pathogen mutation characteristics and transmission potential are introduced as causal variables, the posterior probability of causing a high-risk event is calculated through the Bayesian network, an interpretable risk inference path is constructed, and the quantitative modeling and reasoning identification of the potential risk sample are supported. In combination with the posterior probability and the dynamic determination threshold, the risk level label is generated, and a structured report is output for the access of the border inspection, flow survey, and other systems, forming an early warning integrated mechanism from detection to disposal.

[0085] In a feasible implementation, the transmission graph is constructed based on the pathogen classification information and the mutation information to obtain a dynamic transmission structure graph of the pathogen, which includes:

[0086] Metadata information of each sample is extracted, wherein the metadata information includes the pathogen classification information, the mutation information, sample collection time, sample collection location, and sequencing timestamp;

[0087] According to the pathogen classification information, samples with the same or similar species are classified into the same transmission graph construction unit;

[0088] Based on the mutation information, a mutation feature vector is generated for each sample, wherein the mutation feature vector is encoded by the whole genome mutation site;

[0089] Each sample is represented as a graph node, wherein the graph node includes the pathogen classification information, the mutation information, the sample collection time, the sample collection location, and the sequencing timestamp;

[0090] a genomic similarity score between the sample nodes is calculated based on distance measurement between the mutation feature vectors;

[0091] According to the similarity score and the sampling time and geographical location information, a propagation edge connection is established between the graph nodes under the condition that at least one preset condition is met, wherein the preset condition includes that the similarity score is higher than a set threshold, the sample collection time interval is less than a preset window, and the sample collection location is in a potential propagation path;

[0092] An edge weight is set in the propagation edge connection, wherein the edge weight is calculated from the similarity score, the time difference, and the geographical distance;

[0093] A dynamic propagation structure graph is established by a time sliding window mechanism according to the graph nodes, the propagation edge connection, and the edge weight.

[0094] For example, the embodiment of the application constructs a dynamic propagation graph for describing the propagation relationship of pathogen samples by combining the classification information and mutation information of the pathogen samples. The propagation graph reflects the potential propagation relationship of different pathogen samples in the time, space, and genetic level, and provides a structural basis for subsequent propagation path modeling and risk reasoning.

[0095] First, the metadata information of each sample is extracted from the data of completed nanopore sequencing. The metadata information includes the species classification result of the pathogen, the whole genome mutation information, the collection time and collection location of the sample, and the timestamp of the completed sequencing, and other basic attributes.

[0096] Based on the classification information of the sample, the samples belonging to the same species or variant subtype are classified into the same propagation graph construction unit to ensure the rationality and comparability of the propagation relationship. Then, according to the mutation information, a mutation feature vector is constructed for each sample. The feature vector can be generated by encoding the mutation sites in the whole genome, and is usually expressed in the form of Boolean encoding, site index, or sparse vector, which is convenient for subsequent similarity analysis.

[0097] Subsequently, each sample is abstracted as a node in the graph. Each node contains its classification label, mutation feature vector, collection time and location, and other attributes. Then, based on the distance measurement of the mutation feature vectors between the nodes, the genomic similarity score between the samples is calculated. The similarity measurement can be in the form of Jaccard distance, Mash distance, or SNP site overlap rate.

[0098] On this basis, whether to establish a propagation edge between two sample nodes is judged according to the similarity score, the time difference of collection and the propagation accessibility of geographical location. Specifically, if the similarity score between the samples is higher than a set threshold, or the time interval of collection is within a certain preset time window, or the sample collection sites are in the transmissible path (such as located in the same port or adjacent areas), a propagation edge is established between the two nodes.

[0099] Further, a weight is set for the propagation edge to quantify the credibility of the propagation path. The edge weight can comprehensively consider factors such as sample similarity, time interval and geographical distance, and is generated by using a weighted calculation method. Finally, a sliding time window mechanism is used to update and reconstruct the propagation graph structure in different time periods to form a time-aware dynamic propagation graph sequence, which is used to reflect the propagation evolution trend of the pathogen in the real scene.

[0100] Through the above steps, the completed dynamic propagation graph not only retains the multi-dimensional information of the pathogen in genetic variation, time distribution and geographical distribution, but also provides structured data support for subsequent graph diffusion modeling and causal reasoning.

[0101] In a feasible implementation, the above-mentioned graph diffusion modeling operation is performed on the above-mentioned dynamic propagation structure graph to obtain potential propagation path score information of the pathogen, including:

[0102] The above-mentioned graph nodes in the above-mentioned dynamic propagation structure graph are encoded to construct a graph node feature matrix;

[0103] A time series modeling mechanism is introduced into the above-mentioned dynamic propagation structure graph, and the graph structure is time-ordered based on sample collection time to form a time-aware dynamic propagation graph;

[0104] A graph diffusion model is used to perform graph diffusion modeling according to the above-mentioned time-aware dynamic propagation graph and the above-mentioned graph node feature matrix to obtain graph node propagation embedding vectors;

[0105] According to the above-mentioned graph node propagation embedding vectors and the above-mentioned time-aware dynamic propagation graph, the propagation potential is extracted to obtain node potential propagation vectors;

[0106] The similarity information between each two of the above-mentioned propagation path node potential propagation vectors is calculated to obtain the above-mentioned potential propagation path score information.

[0107] Illustratively, based on the constructed dynamic propagation structure graph, a graph diffusion modeling operation is performed to obtain potential propagation path score information of the pathogen.

[0108] For each node in the propagation graph , three types of information are extracted to form an initial feature vector :

[0109] ;

[0110] where, is the mutation feature vector of the pathogen; is the sampling timestamp; is the sampling location code (e.g. region number, latitude-longitude vector).

[0111] Combine the feature vectors of all nodes into a graph node feature matrix:

[0112] ;

[0113] Based on the sampling time sequence, the propagation graph Introduce a time sliding window , construct a time segmented graph sequence:

[0114] ;

[0115] where, is the node subset within the time window ; is the edge set that satisfies the propagation condition within the window.

[0116] Use a graph diffusion model (such as GAT, TGN, DCRNN) to input the time-aware graph structure and the node feature matrix into the model to obtain the propagation representation vector of each node:

[0117] ;

[0118] is the graph diffusion modeling function, which represents the embedding features of each sample node in the propagation structure.

[0119] For any two nodes and in the graph, based on their propagation embedding vectors, calculate the propagation path score. The scoring method uses a concatenation neural scoring model:

[0120] ;

[0121] where: is the L2 norm; is the vector concatenation; is the sigmoid activation function; is the model trainable parameter, is the graph diffusion embedding vector of node , is the graph diffusion embedding vector of node .

[0122] combining the score values between all pairs of nodes into a propagation score matrix

[0123]

[0124] The matrix represents the potential propagation possibility between any two samples in the graph, which is used for subsequent risk inference and traceability judgment.

[0125] The node feature encoding considers gene mutations, sampling time and geographical location, making the propagation modeling more comprehensive and enhancing the model's ability to depict complex propagation paths.

[0126] The embodiment introduces a time sliding window mechanism to construct a dynamic graph, and uses a time-aware graph diffusion model to effectively reflect the time evolution characteristics in the pathogen propagation process. By representing the propagation potential through graph embedding vectors and calculating the propagation score matrix between node pairs, the propagation possibility between any samples can be quantified, improving the accuracy of propagation chain mining and traceability analysis.

[0127] In one possible implementation, based on the mutation information and the potential propagation path score information, a Bayesian causal model is constructed to obtain the posterior probability of the pathogen causing a high-risk event, including:

[0128] Performing a high-risk mutation feature recognition operation on the mutation feature vector to construct a first causal variable of the Bayesian causal model;

[0129] Performing a propagation potential recognition operation on the potential propagation path score information to construct a second causal variable of the Bayesian causal model;

[0130] Based on the first causal variable and the second causal variable, a Bayesian causal graph model is constructed, wherein the Bayesian causal model includes a posterior variable node representing a high-risk event and a directed edge from the causal variable node to the posterior variable node;

[0131] According to the prior conditional probability and the joint distribution inference rule, the conditional probability of the posterior variable node in the Bayesian causal graph model is calculated to obtain the posterior probability of the pathogen causing a high-risk event.

[0132] For example, first, for each sample's mutation feature vector, a high-risk mutation recognition operation is performed. This operation includes the following sub-steps: obtaining a high-risk mutation reference library containing mutation sites known to have enhanced transmission or immune escape ability. Compare the sample's mutation feature vector with the reference library to identify whether it contains a high-risk mutation; according to the comparison result, judge whether the sample hits a key risk mutation, output the first causal variable​​ If the number of hits of high-risk mutation sites exceeds a set threshold, or a specific combination of mutations is hit, it is considered to have a high-risk mutation feature, set ; otherwise set .

[0133] Subsequently, based on the potential propagation ability of the sample in the propagation score graph, a second causal variable is constructed. The specific operation is as follows: the following indicators of the sample are extracted from the propagation path score matrix, for example: the number of outward propagation paths, the maximum propagation score value, and the structural centrality index (such as out-degree, betweenness centrality) in the propagation graph. According to the above indicators and the preset rules, the propagation potential is judged, if the propagation ability of the sample exceeds the set threshold, set , otherwise set .

[0134] According to the above two causal variables, a simplified Bayesian causal graph structure is constructed:

[0135] ;

[0136] Among them: is a high-risk mutation feature (Boolean variable); is a propagation potential feature (Boolean variable); is whether to trigger a high-risk event (Boolean posterior variable). The graph represents that the occurrence of the high-risk event is jointly affected by the high-risk mutation and the propagation potential.

[0137] According to the Bayesian inference rule, the conditional probability of the posterior variable is calculated, that is:

[0138] ;

[0139] This can be achieved by the following steps: based on training samples or historical data, a joint probability table or CPT (conditional probability table) is constructed. The maximum likelihood estimation or Bayesian estimation is used to obtain the prior probability. For the current sample, the value of is substituted into the corresponding ;

[0140] The posterior probability of the sample triggering a high-risk event is output as the basis for risk assessment.

[0141] Through this embodiment, a double-causal explanation of the pathogen sample from the mutation level and the propagation level can be realized, and the Bayesian graph modeling method is used to complete the quantitative inference of the high-risk event, providing a basis for subsequent early warning classification and key monitoring.

[0142] In an implementable embodiment, the high-risk mutation feature identification operation on the mutation feature vector is performed to construct the first causal variable of the Bayesian causal model, including:

[0143] A high-risk mutation site reference library is obtained.

[0144] The mutation feature vector and the high-risk mutation site reference library are compared to identify high-risk mutation related sites to obtain a comparison result.

[0145] According to the comparison result and a preset mutation risk determination rule, statistical information of a mutation feature matching number is obtained.

[0146] According to the statistical information, the first causal variable of the Bayesian causal model is determined.

[0147] For example, a high-risk mutation site reference library is constructed or loaded, which can be derived from high pathogenicity and high transmissibility mutation sites reported in domestic and foreign public databases, and can also be derived from epidemiological statistical analysis results, and can also be derived from laboratory function verification data, etc.

[0148] The reference library may be in the form of:

[0149] ;

[0150] wherein each represents a key mutation site, including site number, gene name, risk level, and other meta information.

[0151] The mutation feature vector of the target sample is compared with the reference library bit by bit. If the sample has a mutation at the site, it is considered to hit a high-risk mutation. The matching situation between the sample and the reference library is output, for example, the number of hits , and the hit combination can also be recorded.

[0152] Based on the comparison result, the following information is counted: the total number of hit mutations , the risk level weighted sum of hit mutations , whether a specific mutation combination (such as “E484K+N501Y” joint mutation) is hit, and whether a known domain sensitive site (such as the RBD region) is hit.

[0153] The risk score can be represented as:

[0154] ;

[0155] wherein, is the risk weight of the mutation site in the reference library.

[0156] According to the preset mutation risk determination rule, the sample is determined as "whether it has high-risk mutation characteristics" as the first causal variable .

[0157] For example, the determination rule can include:

[0158] 1. If , set .

[0159] 2. If , set , otherwise, set . Wherein, is the minimum hit number threshold (such as 3), is the risk score threshold (such as the 90th percentile or the empirical boundary value).

[0160] In a feasible implementation, the determination step of the above-mentioned preset mutation risk determination rule includes:

[0161] Based on the mutation risk score, a sliding time window sample set is constructed, wherein the sliding time window sample set is all pathogen samples received within a preset time window, and the mutation risk score is calculated from the mutation feature vector;

[0162] The distribution of the mutation risk score in the sliding time window sample set is statistically analyzed, and the high quantile value of the distribution is calculated as a dynamic threshold;

[0163] The mutation risk score of the current sample is compared with the dynamic threshold to determine whether to update the preset mutation risk determination rule.

[0164] For example, in a feasible implementation, in order to improve the flexibility and timeliness of high-risk mutation identification, a sliding time window mechanism can be introduced to realize dynamic updating and adaptive optimization of the mutation risk determination rule, which includes the following steps:

[0165] A sliding time window (such as 7 days, 14 days, 1 month) is set, denoted as , which contains all pathogen sample data received by the system in this time period.

[0166] For each sample in the window, its corresponding mutation risk score is calculated through its mutation feature vector , for example:

[0167] ;

[0168] Wherein: is the sample the number of hit high-risk mutation points; the risk weight of each point in the reference library; is the total mutation risk score of the current sample.

[0169] constitute a window sample score set :

[0170] ;

[0171] perform distribution statistics on all risk scores in the sliding time window sample set , and calculate the high quantile (such as 90 quantile, 95 quantile, etc.) as a dynamic judgment threshold:

[0172] ;

[0173] wherein: is the quantile point setting parameter; is the dynamic risk score threshold of the current time period.

[0174] The dynamic threshold reflects the high abnormality of the current mutation risk score in the whole sample distribution, which helps to adaptively cope with the rising trend of new variant strains.

[0175] For the current new sample , calculate its mutation risk score , and compare it with the current dynamic threshold . If , it is considered that the sample has high-risk mutation characteristics, and the first causal variable is set. If it is lower than the threshold, set . At the same time, the preset judgment rule for the next round of identification can be updated accordingly.

[0176] The mechanism proposed in this embodiment can continuously adjust the definition boundary of “high-risk mutation” of the model, adapt to the mutation spectrum changes in different stages and different regions. The model no longer relies on fixed empirical thresholds, and can automatically adjust the judgment rules according to the virus evolution trend. The sliding time window design makes the model respond to sample risk changes in real time, establishes a judgment standard through statistical distribution, reduces the influence of individual extreme values, does not require frequent manual intervention, and is conducive to long-term operation of the system.

[0177] In a feasible implementation manner, the operation of identifying the transmission potential of the potential transmission path score information to construct the second causal variable of the Bayesian causal model comprises:

[0178] determining transmission potential judgment threshold information;

[0179] According to the above propagation potential determination threshold information and the above potential propagation path score information, propagation path quantity information, maximum propagation probability information and structural centrality indicators are obtained;

[0180] According to the above propagation path quantity information, the above maximum propagation probability information, the above structural centrality indicators and a preset propagation potential determination rule, a second causal variable of the Bayesian causal model is determined.

[0181] For example, a set of indicator thresholds for determining "high propagation potential" is set or dynamically generated. The propagation path quantity threshold is the minimum requirement for the effective path quantity of the sample outward diffusion. The maximum propagation probability threshold is the lower limit of the maximum single path score corresponding to the sample. The structural centrality threshold is the lower limit of the judgment of the graph structure centrality (such as out-degree, betweenness centrality, PageRank). The above threshold can be a static empirical value, or can be generated in combination with a dynamic sliding window mechanism (such as adaptive adjustment using 90th percentile).

[0182] For the target sample , the following indicators are extracted from the propagation path score matrix :

[0183] Propagation path quantity :

[0184] ;

[0185] Among them, is the path score determination threshold, is the indicator function.

[0186] Maximum propagation probability :

[0187] ;

[0188] Structural centrality : optional indicators include out-degree centrality, PageRank value, K-shell value, etc., reflecting the "strategic position" of the node in the propagation graph.

[0189] According to the extracted indicators and the preset or dynamically generated threshold set , the following strategy is used to determine:

[0190] If:

[0191] ;

[0192] , it is considered that the sample has high propagation potential, and the second causal variable is set.

[0193] Otherwise set The above causal variables are used to model the risk spreading ability of pathogenic samples at the level of the transmission path.

[0194] Through this embodiment, the diffusion intensity and centrality position of samples in the transmission graph can be systematically analyzed. The transmission score results are quantified as logical variables , which facilitates the inclusion of Bayesian causal graph modeling. This improves the identification efficiency of "potential super-spreading nodes" or "network key transmitters", and provides a model basis for accurate tracing and key early warning.

[0195] In a feasible implementation, the second causal variable of the Bayesian causal model is determined according to the transmission path quantity information, the maximum transmission probability information, the structural centrality index, and a preset transmission potential determination rule, comprising:

[0196] According to the transmission path quantity information, the maximum transmission probability information, and the structural centrality index, a transmission probability driving term, a transmission depth adjustment term, a target node importance weighting term, and a spatial accessibility adjustment term are determined;

[0197] According to the transmission probability driving term, the transmission depth adjustment term, the target node, the transmission potential score function, and the preset transmission potential determination rule, the second causal variable of the Bayesian causal model is determined.

[0198] For example, in a feasible implementation, to realize the quantitative evaluation of the potential diffusion ability of pathogenic samples in the transmission network, the present application constructs a transmission potential score function based on transmission path quantity information, maximum transmission probability information, and structural centrality index. And combined with a preset transmission potential determination rule, the second causal variable in the Bayesian causal model is determined. Specifically, the following steps are included:

[0199] (1) Constructing a transmission capacity sub-item index: for each pathogenic sample node , the following three types of core indicators are extracted from the transmission graph and score matrix:

[0200] 1. Transmission path quantity : represents the total number of transmission paths (i.e. the number of edges that exceed the set threshold) sent by the node ;

[0201] 2. Maximum transmission probability : represents the edge with the highest score in the path sent by the node ;

[0202] 3. Structural centrality index : node The centrality measure in the propagation graph, such as PageRank value or betweenness centrality, etc.

[0203] In addition, the spatial correlation information can also be introduced as a "spatial accessibility adjustment term", such as a score weighted based on the number of areas covered by the propagation path or the geographical distance.

[0204] (2) Build a propagation potential score function: based on the extracted multiple indicators, build a fusion type propagation potential score function to comprehensively evaluate the propagation ability of the node. The score function is expressed as follows:

[0205] ;

[0206] Wherein: : represents the Sigmoid activation function, which ensures that the final score falls in the interval ; : is the weight parameter of the five features; is the exponential coefficient for adjusting the non-linear enhancement effect;

[0207] In the above score function is the propagation depth adjustment term, which is used to measure the number of paths that the node can reach in the propagation graph, reflecting the breadth or coverage ability of the propagation. The more the number of paths, the wider the range of the propagation, and the higher the score.

[0208] is the propagation probability driven term, which represents the probability value of the highest scored path in the outward propagation path of the node, reflecting the ability of the potential "strongest propagation" path. The higher the score, the stronger the single-point high-efficiency propagation potential.

[0209] is the target node importance weighting term, wherein represents the structural centrality of the node, such as PageRank value or betweenness centrality. This term is used to describe the importance of the node in the graph structure. The higher the centrality, the more strategic the node is in a key position, and the more strategic the node is in a key position.

[0210] is the breadth and strength cross term, which means that the node has both sufficient path quantity (breadth) and strong high-score path (strength), and the product of the two enhances its ability to spread widely but accurately in the transmission network.

[0211] is the structure and strength cross term, which means that the node not only has strong path, but also is in an important position in the propagation graph, which is a further strengthening of the "super propagation point" feature. Nodes with high strength and centrality will get higher scores.

[0212] In an implementable embodiment, the risk level label is generated according to the posterior probability, and the pathogen analysis report and warning information are output, including:

[0213] The dynamic level threshold information is obtained.

[0214] The risk level label is determined according to the posterior probability and the dynamic level threshold information.

[0215] The pathogen analysis report and warning information are output according to the risk level label and the sample analysis result.

[0216] For example, the posterior probability values of multiple pathogen samples output based on the Bayesian causal model are combined with known historical label information (such as aggregation and transmission events, case pathogenicity level, etc.) to construct the judgment standard of risk level division.

[0217] Specifically, the posterior probability value corresponding to the maximum Youden index can be selected as the first-level risk level judgment threshold by constructing a receiver operating characteristic curve (ROC); in the scenario where multiple risk levels need to be divided, the posterior probability intervals corresponding to multiple risk levels are dynamically determined based on the macro-average AUC index or distribution statistics (such as 90th percentile, 75th percentile, 50th percentile).

[0218] The level interval can be set as, for example, high risk for posterior probability ≥ 0.85; medium-high risk for 0.70 ≤ posterior probability < 0.85; medium risk for 0.50 ≤ posterior probability < 0.70; and low risk for posterior probability < 0.50. The threshold value can be updated in real time based on the sample data distribution within a sliding time window.

[0219] For the current sample to be evaluated, the posterior probability value output by the Bayesian causal model is matched with the currently determined dynamic level threshold to assign a corresponding risk level label. For example, if the posterior probability is 0.91, the label “high risk” is assigned; if the posterior probability is 0.68, the label “medium risk” is assigned. The above label information will serve as the grading basis for downstream pathogen management strategies.

[0220] Combined with the metadata information obtained in the previous processing process of the sample (such as classification results, mutation characteristics, transmission score vector, transmission potential sub-items, causal variable values, etc.), a structured output template is constructed to form a pathogen analysis report and warning summary information containing the following contents: sample number, collection time and sampling location; classification species and its confidence; hit mutation site information and risk mutation identification; obtained transmission potential score and sub-item score; Bayesian posterior probability and key causal variable values; automatically generated risk level label; warning suggestions generated in combination with rules (such as review, isolation, reporting, warning prompts, etc.).

[0221] The pathogen analysis report can be output as a graphic PDF by a visualization component, or output as JSON, XML, table file, etc. for access to the border inspection monitoring system, flow survey platform, emergency response system, etc.

[0222] In a feasible implementation, the specific determination of the dynamic level threshold information includes:

[0223] A historical pathogen sample set containing the above posterior probability values and corresponding true risk labels is constructed.

[0224] A receiver operating characteristic curve is constructed for the historical sample set with the above posterior probability as the prediction variable and the above risk label as the true classification standard.

[0225] The sensitivity and specificity indicators corresponding to different posterior probability thresholds are calculated, and the Youden index formed based on the sum of sensitivity and specificity minus one is calculated.

[0226] The posterior probability value with the maximum Youden index is taken as the dynamic judgment threshold for the first risk level division.

[0227] The above posterior probability value is optimized with multiple thresholds in combination with the macro-average AUC index and the classification accuracy evaluation index to form the probability interval corresponding to the high, medium, medium-high and low risk levels.

[0228] For example, first, based on the completed pathogen sample data, a historical sample set is constructed. Each sample in the sample set includes a posterior probability value calculated by a Bayesian causal model and a true risk label corresponding to the sample (such as whether it has triggered a transmission event, whether it is highly pathogenic, etc.). The sample set is used to train and optimize the correspondence between the posterior probability and the risk level.

[0229] The posterior probability is taken as the prediction variable, and the historical true risk label is taken as the classification standard. The model performance indicators including sensitivity (Sensitivity) and specificity (Specificity) are calculated under different probability thresholds. Based on the above results, a receiver operating characteristic curve (ROC) is drawn to comprehensively evaluate the classification ability of the model under each probability threshold.

[0230] For each candidate posterior probability interval, the corresponding Youden index is calculated, defined as:

[0231] ;

[0232] wherein Sensitivity is the sensitivity indicating the ability of the model to identify positive samples (e.g. high-risk pathogens), and Specificity is the specificity indicating the ability of the model to identify negative samples (e.g. low-risk pathogens).

[0233] The posterior probability interval that maximizes the Youden index is selected as the dynamic judgment threshold for the classification of the first-level risk grade. This value has the optimal balance and resolution capability under the current state of sample distribution.

[0234] When finer classification of high, medium, high, low risk grades is performed, further multi-classification model evaluation indicators such as Macro-AUC, weighted accuracy, F1-score, etc. are used. Through the analysis of the distribution density of posterior probability and the classification performance in different intervals, the probability boundary points corresponding to multiple risk grades are determined. The interval division can be regenerated according to the latest samples in the sliding window in each round of system update, thereby forming an adaptive grade judgment mechanism.

[0235] The dynamic classification mechanism proposed in the embodiment makes the risk grade classification more timely and scenario-adaptive. The introduction of ROC and Youden index judgment methods ensures the optimality of the classification threshold in the statistical sense. Through multi-level classification optimization, the risk management strategy can be refined to provide more accurate decision support for the subsequent early warning response system.

[0236] The above embodiments are only used to illustrate the technical solutions of the present application, but not limit them; although the present application has been described in detail with reference to the foregoing embodiments, it should be understood by those skilled in the art that the technical solutions recorded in the foregoing embodiments can be modified, or some technical features can be replaced by equivalents; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present application.

Claims

1. A pathogen gene sequence detection analysis method based on nanopore sequencing, characterized by, The method comprises the following steps: nanopore sequencing of collected samples to obtain raw sequencing data; quality screening of the raw sequencing data to obtain screened sequencing data; species classification and mutation point monitoring operations on the screened sequencing data to obtain pathogen classification information and mutation information; constructing a transmission map based on the pathogen classification information and mutation information to obtain a dynamic transmission structure diagram of the pathogen; performing graph diffusion modeling on the dynamic transmission structure diagram to obtain potential transmission path score information of the pathogen; constructing a Bayesian causal model based on the mutation information and the potential transmission path score information to obtain the posterior probability of the pathogen triggering a high-risk event; generating a risk level label according to the posterior probability, and outputting a pathogen analysis report and early warning information; the method of constructing a transmission map based on the pathogen classification information and mutation information to obtain a dynamic transmission structure diagram of the pathogen comprises: extracting metadata information of each sample, wherein the metadata information includes the pathogen classification information, mutation information, sample collection time, sample collection location, and sequencing timestamp; grouping samples with the same or similar species into the same transmission map construction unit according to the pathogen classification information; generating a mutation feature vector for each sample based on the mutation information, wherein the mutation feature vector is encoded by the whole genome mutation site; representing each sample as a graph node, wherein the graph node includes the pathogen classification information, the mutation information, the sample collection time, the sample collection location, and the sequencing timestamp; calculating the genomic similarity score between sample nodes based on the distance measurement between the mutation feature vectors; establishing a transmission edge connection between the graph nodes under the condition that at least one preset condition is met according to the similarity score and sampling time and geographical location information, wherein the preset conditions include that the similarity score is higher than a set threshold, the sample collection time interval is less than a preset window, and the sample collection location is in the potential transmission path; setting an edge weight value in the transmission edge connection, wherein the edge weight value is calculated from the similarity score, the time difference, and the geographical distance; establishing the dynamic transmission structure diagram according to the graph nodes, the transmission edge connection, and the edge weight value through a time sliding window mechanism; the method of performing graph diffusion modeling on the dynamic transmission structure diagram to obtain potential transmission path score information of the pathogen comprises: encoding the graph nodes in the dynamic transmission structure diagram to construct a graph node feature matrix; introducing a time series modeling mechanism into the dynamic transmission structure diagram, time sorting the graph structure based on sample collection time to form a time-aware dynamic transmission graph; performing graph diffusion modeling using a graph diffusion model based on the time-aware dynamic transmission graph and the graph node feature matrix to obtain a graph node transmission embedding vector; extracting transmission potential based on the graph node transmission embedding vector and the time-aware dynamic transmission graph to obtain a node potential transmission vector; Calculate similarity information between every two potential propagation path nodes to obtain the potential propagation path score information.

2. The pathogen gene sequence detection analysis method based on nanopore sequencing according to claim 1, characterized in that, Based on the mutation information and the potential propagation path score information, a Bayesian causal model is constructed to obtain the posterior probability of a pathogen triggering a high-risk event, including: Performing a high-risk mutation feature recognition operation on the mutation feature vector to construct a first causal variable of the Bayesian causal model; Performing a propagation potential recognition operation on the potential propagation path score information to construct a second causal variable of the Bayesian causal model; Based on the first causal variable and the second causal variable, the Bayesian causal graph model is constructed, wherein the Bayesian causal model includes a posterior variable node representing a high-risk event and a directed edge from the causal variable node to the posterior variable node; According to the prior conditional probability and the joint distribution inference rule, the conditional probability of the posterior variable node in the Bayesian causal graph model is calculated to obtain the posterior probability of the pathogen triggering a high-risk event.

3. The nanopore sequencing-based pathogen gene sequence detection analysis method according to claim 2, characterized in that, The high-risk mutation feature recognition operation on the mutation feature vector to construct the first causal variable of the Bayesian causal model includes: Obtain a high-risk mutation point reference library; Perform an alignment operation on the mutation feature vector and the high-risk mutation point reference library to identify high-risk mutation related positions to obtain an alignment result; According to the alignment result and a preset mutation risk determination rule, statistical information of the number of mutation feature matches is obtained; According to the statistical information, the first causal variable of the Bayesian causal model is determined.

4. The pathogen gene sequence detection analysis method based on nanopore sequencing according to claim 3, characterized in that, The determination step of the preset mutation risk determination rule includes: Based on the mutation risk score, a sliding time window sample set is constructed, wherein the sliding time window sample set is all pathogen samples received within a preset time window, and the mutation risk score is calculated from the mutation feature vector; Distribute statistics of the mutation risk score in the sliding window sample set, and calculate the high quantile value of the distribution as a dynamic threshold; Compare the mutation risk score of the current sample with the dynamic threshold to determine whether to update the preset mutation risk determination rule.

5. The nanopore sequencing-based pathogen gene sequence detection analysis method according to claim 3, characterized in that, The propagation potential recognition operation on the potential propagation path score information to construct the second causal variable of the Bayesian causal model includes: Determine the propagation potential determination threshold information; According to the propagation potential determination threshold information and the potential propagation path score information, obtain propagation path quantity information, maximum propagation probability information and structural centrality index; According to the propagation path quantity information, the maximum propagation probability information, the structural centrality index and the preset propagation potential determination rule, the second causal variable of the Bayesian causal model is determined.

6. The nanopore sequencing-based pathogen gene sequence detection analysis method according to claim 5, characterized in that, The determination of the second causal variable of the Bayesian causal model according to the propagation path quantity information, the maximum propagation probability information, the structural centrality index and the preset propagation potential determination rule includes: According to the propagation path information, the maximum propagation probability information, and the structural centrality index, determine a propagation probability driving item, a propagation depth adjustment item, a target node importance weighting item, and a spatial reachability adjustment item; According to the propagation probability driving item, the propagation depth adjustment item, the target node importance weighting item, and the spatial reachability adjustment item, determine a propagation potential score function; According to the propagation potential score function and the preset propagation potential determination rule, determine a second causal variable of the Bayesian causal model.

7. The nanopore sequencing-based pathogen gene sequence detection analysis method according to claim 1, characterized in that, According to the posterior probability, generate a risk level label, and output a pathogen analysis report and early warning information, including: Obtaining dynamic level threshold information; According to the posterior probability and the dynamic level threshold information, determine the risk level label; According to the risk level label and the sample analysis result, output the pathogen analysis report and early warning information.

8. The nanopore sequencing-based pathogen gene sequence detection analysis method according to claim 7, characterized in that, The specific determination steps of the dynamic level threshold information include: Constructing a historical pathogen sample set containing the posterior probability value and the corresponding true risk label; For the historical sample set, construct a receiver operating characteristic curve with the posterior probability as the prediction variable and the risk label as the true classification standard; Calculate the sensitivity and specificity indexes corresponding to different posterior probability thresholds, and form the Youden index based on the sum of sensitivity and specificity minus one; The posterior probability value with the maximum Youden index is taken as the dynamic judgment threshold of the first risk level division; Combine the macro average AUC index and the classification accuracy evaluation index to optimize the posterior probability value with multiple thresholds to form the probability interval corresponding to the high, medium and low risk levels.

Images