Multi-dimensional data collaborative fusion multi-pathogenic bacterium invasion mechanism comparative analysis method

Abstract

The invention belongs to the technical field of multi-pathogenic-bacterium invasion mechanism comparative analysis methods, and particularly relates to a multi-dimensional data collaborative fusion multi-pathogenic-bacterium invasion mechanism comparative analysis method, which comprises the following steps of: eliminating noise and scale difference for genome, transcriptome and proteome multi-source heterogeneous data by adopting a combination standardization method; extracting corresponding basic features through a bioinformatics tool, and integrating and constructing a multi-modal original feature library; a double-branch attention mechanism is introduced, intra-modal features are enhanced respectively, inter-modal bidirectional association is established, and a cross-modal feature set is generated through weighted fusion; abstracting a pathogenic bacterium invasion process into a time sequence interaction network, embedding a gating circulation unit based on a time sequence diagram convolution network to establish a dynamic framework, and dynamically adjusting an edge weight to describe network evolution; extracting a pathogenic bacterium virulence gene expression time sequence curve and host immune response dynamic data, and supplementing to the time sequence interaction network.

Description

Technical Field

[0001] This invention belongs to the technical field of comparative analysis methods for invasion mechanisms of multiple pathogens, and particularly relates to a comparative analysis method for invasion mechanisms of multiple pathogens through multidimensional data synergistic fusion. Background Technology

[0002] In current comparative analyses of invasion mechanisms of multiple pathogens, the integration and effective utilization of multi-source heterogeneous data has become a core bottleneck. Genomic, transcriptomic, and proteomic data suffer from noise interference and significant scale differences. Traditional standardization methods struggle to completely eliminate data bias, and single-modal feature extraction ignores the intrinsic connections between molecular mechanisms at different levels, leading to incomplete feature expression. Furthermore, existing technologies lack efficient cross-modal fusion strategies; redundant features within a modality are not effectively suppressed, and bidirectional correlation mapping between modalities is insufficient, resulting in low biological relevance and discriminative power of fused features, failing to accurately support subsequent differential analysis.

[0003] Traditional analytical methods often employ static models to characterize the pathogen invasion process, failing to fully consider the temporal dynamics of invasion and making it difficult to capture the evolutionary patterns of molecular interaction networks at different stages. Furthermore, existing methods frequently focus only on single-dimensional differences, lacking a comprehensive analysis of the entire chain of characteristics, processes, and differences. This results in insufficient precision in quantifying differences, ambiguous interpretations of their biological significance, and an inability to systematically reveal the core differences in the invasion mechanisms of various pathogens. Consequently, this limits our in-depth understanding of pathogenic mechanisms and the progress of related prevention and control technologies and drug development. Summary of the Invention

[0004] The purpose of this invention is to address the aforementioned technical problems by providing a comparative analysis method for the invasion mechanisms of multiple pathogens through multi-dimensional data synergistic fusion.

[0005] In view of this, the present invention provides a comparative analysis method for invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation. The method uses a combined normalization method to eliminate noise and scale differences in multi-source heterogeneous data of genome, transcriptome, and proteome, extracts corresponding basic features through bioinformatics tools, and integrates and constructs a multimodal original feature library. A dual-branch attention mechanism is introduced to strengthen intra-modal features and establish bidirectional inter-modal correlations, and weighted fusion is used to generate a cross-modal feature set; The pathogen invasion process is abstracted into a temporal interaction network. A dynamic framework is built by embedding gated recurrent units into a temporal graph convolutional network, and the edge weights are dynamically adjusted to characterize the network evolution. Extract the expression time-series curves of pathogenic bacteria virulence genes and dynamic data of host immune response, supplement them into the temporal interaction network, highlight the contribution of time nodes through temporal attention, and realize the dynamic update of node status; Based on dynamic fusion characteristics and temporal network evolution results, a full-chain correlation model is constructed. It analyzes from three levels: characteristics, process, and differences, quantifies the differences in pathogen invasion mechanisms, and outputs an explanatory report.

[0006] The beneficial effects of this invention are as follows: This invention accurately eliminates noise and scale differences in multi-source heterogeneous data by combining standardization methods, and efficiently extracts multimodal basic features by combining bioinformatics tools. The constructed original feature library has both completeness and effectiveness. The introduction of the dual-branch attention mechanism not only strengthens core functional features and suppresses redundant information through intramodal self-attention, but also establishes full-dimensional bidirectional associations by using cross-modal cross attention. Combined with the adaptive weight adjustment mechanism, it achieves optimal feature fusion. The generated cross-modal feature set has stronger biological relevance and recognizability, and completely solves the pain points of incomplete feature expression and poor modality fusion effect of traditional methods.

[0007] Based on the dynamic framework of temporal graph convolutional networks and gated recurrent units, this invention transforms the pathogen invasion process into a quantifiable temporal interaction network evolution process. Through adaptive edge weight updates and dynamic node state updates, it accurately captures the molecular interaction patterns at each stage of invasion. The whole-chain association model systematically analyzes the differences in invasion mechanisms from three levels: features, processes, and differences. It achieves quantitative characterization and biological interpretation of the differences in invasion mechanisms. The output interpretation report is both accurate and readable, which not only greatly improves the scientificity and efficiency of comparative analysis of multiple pathogen invasion mechanisms, but also provides strong technical support for pathogenic mechanism research, prevention and control technology development, and targeted drug design. Attached Figure Description

[0008] Figure 1 This is a flowchart illustrating the comparative analysis method for the invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation, as presented in this invention. Detailed Implementation

[0009] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0010] Figure 1 This is a flowchart illustrating a method for comparative analysis of multiple pathogenic bacteria invasion mechanisms based on illumination-adaptive plane representation, provided in one embodiment of the present invention. The following is a detailed description of this method.

[0011] Step S110: Use a combined normalization method to eliminate noise and scale differences for multi-source heterogeneous data of genome, transcriptome, and proteome, extract corresponding basic features through bioinformatics tools, and integrate them to construct a multimodal original feature library; A combined strategy of modality-specific standardization and cross-modality unified scale was employed to process multi-source heterogeneous data: genomic data were first standardized using Z-scores to eliminate sequencing depth differences, and then median centering was used to correct sample bias; transcriptome data underwent TMM standardization to balance the influence of gene length, combined with quantile standardization to unify sample distribution; proteome data were first log2 transformed to compress the quantification range, and then quantile standardization was used to eliminate systematic errors in the mass spectrometry instrument. Bioinformatics tools used included FastQC (quality control preprocessing), Salmon (transcriptome quantification), MaxQuant (proteome quantification), and BWA (genomic alignment) to extract fundamental characteristics such as gene length, GC content, and sequencing coverage from the genome; FPKM value, transcript integrity, and exon utilization from the transcriptome; and peptide matching number, protein abundance, and modification site information from the proteome. The Python Pandas tool was used to standardize the feature data format (numerical features were retained to 4 decimal places, and categorical features were encoded as integers). Cross-modal feature associations were established through unique sample identifiers (such as strain ID-sample number). A data table containing 120 samples and 3800 basic features was constructed. After noise filtering, the feature CV values ​​(coefficient of variation) were all controlled within 20%. Finally, an original feature library covering sample source information, multimodal feature matrices, and feature function descriptions was formed.

[0012] Step S111: For the raw genome sequencing data, quality control methods such as adapter removal and low-quality base splicing are used. For the transcriptome RNA-seq data, adapter sequence and ploy-N tail removal and short sequence filtering are performed by detecting base quality distribution, GC content and sequence repetition rate. For the proteome mass spectrometry data, peptide identification and quantification are performed, and low-reliability protein entries are removed by database matching and false positive result filtering. For the raw genome sequencing data (Illumina HiSeq platform, 150bp paired ends), quality control was performed using the Trimmomatic tool: the adapter sequence was set to the Illumina TruSeq universal adapter (AGATCGGAAGAGC), and low-quality fragments were cut using parameters LEADING:3 (removing bases with a quality value <3 at the 5' end), TRAILING:3 (removing bases with a quality value <3 at the 3' end), SLIDINGWINDOW:4:15 (sliding window of 4bp, average quality value ≥15), and MINLEN:36 (retaining sequences of length ≥36bp). After quality control, the Q30 of the data was ≥90%. For transcriptome RNA-seq data, FastQC was used to detect base quality distribution (ensuring Phred values ​​≥20 at each position), and samples with abnormal GC content (deviation from the population mean ±10%) were directly removed. Picard MarkDuplicates was used to detect sequence repetition rate (retaining samples with repetition rate <30%), and the Cutadapt tool was used with parameters --adapter AGATCGGAAGAGC --trim-n to remove adapter sequences and poly-N tails, filtering short sequences <50bp in length. For proteome mass spectrometry data (LTQ-Orbitrap platform), MaxQuant software was used for peptide identification, matching to the Uniprot pathogen-specific database (containing 12,000 protein sequences), and false positive filtering criteria were set for peptide FDR ≤1% and protein FDR ≤1%. Low-reliability protein entries with peptide matching count <2 or quantitative value missing rate >30% were removed, ultimately obtaining a high-quality protein dataset.

[0013] Step S112: Sequence alignment correction is used for genomic data. The quality-controlled sequences are aligned to the corresponding pathogenic bacteria reference genome. Unaligned, multiple-aligned, and PCR repetitive sequences are removed, and base quality values ​​are corrected. For transcriptome data, a combination of expression level normalization and batch effect correction is used to achieve data scale uniformity. For proteome data, a combination of protein structure homology filtering and quantitative standardization is used to remove homologous redundant proteins, correct quantitative differences between samples, and fill in missing values. Genomic data correction employed the BWA-MEM tool to align quality-controlled sequences to a pathogenic bacterial reference genome (e.g., *E. coli* MG1655, GenBank accession number NC_000913.3), with parameters set to -k32-M to improve alignment accuracy. Unaligned sequences were filtered using Samtools view -F 4, and multiple alignment sequences (mapped quality value <20) were removed using -q20. PCR repetitive sequences were removed using Picard MarkDuplicates. The GATK BaseRecalibrator tool was used to correct base quality values ​​against a known SNP database (e.g., dbSNP), reducing the base error rate to below 0.01%. Transcriptome data underwent a combined TMM normalization and ComBat batch correction strategy: gene expression levels were first corrected using the TMM method in the edgeR package, then the ComBat function in the sva package was used to eliminate systematic errors from different sequencing batches, unifying the expression levels of all samples to FPKM values. After batch correction, the Pearson correlation coefficient between samples was ≥0.85. The proteomic data were clustered using the CD-HIT tool (sequence identity ≥90%) to remove homologous redundant proteins. The median centering method was used for quantitative standardization. For missing values, the k-nearest neighbor algorithm (k=5) was used to fill in the missing values ​​to ensure that the missing rate of protein quantitative values ​​for each sample was <15% and the standard deviation of protein abundance distribution among samples after correction was ≤0.3.

[0014] Step S113: Genomic feature extraction involves homology comparison between the corrected sequence and the virulence gene database to extract virulence gene-related sequences and functional domain features. Transcriptome feature extraction uses differential expression analysis combined with functional enrichment analysis to screen transcripts and extract their related features. Proteome feature extraction involves screening differentially expressed proteins, predicting their functional domains, and extracting domain-related features and quantitative and functional annotation information. The corrected sequences were aligned with VFDB (Virulent Factor Database, v2023) using the BLASTn tool, with screening criteria of E-value ≤ 1e-5 and sequence identity ≥ 80%. Features such as length, GC content, and homology score of virulence gene-related sequences were extracted. The SMART tool was used to predict the functional domains of virulence genes (such as serine protease domains and adhesion domains), recording features such as domain location and the number of conserved sites, resulting in 120 core virulence-related features. Transcriptome feature extraction: Differential expression analysis was performed using the DESeq2 package, screening for differentially expressed transcripts with |log2FC| ≥ 1 and FDR < 0.05. GO (gene ontology) and KEGG (metabolic pathway) functional enrichment analyses were performed using the clusterProfiler package, extracting features such as FPKM value, fold change, enriched pathway ID, and pathway enrichment score of differentially expressed transcripts. Structural features such as exon number and alternative splicing type were also retained, resulting in a total of 860 transcriptome features extracted. Proteome feature extraction: Differentially expressed proteins with p < 0.05 and |log2FC| ≥ 1.2 were screened using t-test. The functional domains of the proteins (such as ATP-binding domains and signal peptide domains) were predicted using the InterProScan tool. Features such as domain ID, peptide matching number, protein abundance value, and GO annotation information were extracted, resulting in a total of 680 core proteome features.

[0015] Step S114: Standardize the format of the basic features of each modality, associate and map the features of different modalities through sample identifiers, and construct a sample and feature data table; adopt a screening strategy of low-variation feature filtering and high-correlation feature removal to retain effective features; finally, construct an original feature library containing sample metadata, multimodal feature matrix and feature annotation information for subsequent cross-modal fusion and dynamic analysis.

[0016] Genomic features are named using the format `genome_feature type_gene ID` (e.g., `genome_domain length_VF0012`), transcriptome features using `transcriptome_feature type_transcript ID` (e.g., `transcriptome_FPKM_T0003`), and proteome features using `proteome_feature type_protein ID` (e.g., `proteome_abundance value_P01234`). All numerical features are uniformly retained to four decimal places, and categorical features are converted using one-hot encoding. Using strain ID-sample number as the unique association key, the genomic, transcriptomic, and proteomic feature data tables are left-joined using Python Pandas to construct an integrated data table containing sample identifiers and multimodal features. Feature selection employs a low-variance filtering and high-correlation removal strategy: the `varianceThreshold` tool in scikit-learn is used to remove low-variance features with a coefficient of variation (CV) < 0.1, and Pearson correlation analysis is used to remove redundant features with pairwise correlation coefficients (r > 0.85), ultimately retaining 2100 effective features. The constructed original feature library contains three core components: sample metadata (including 12 items such as strain type, culture conditions, and collection time), a multimodal feature matrix (200 samples × 2100 features), and feature annotation information (including feature source modality, functional description, extraction tools, and screening criteria). It is stored in HDF5 format (supporting efficient reading and writing and cross-platform access). This feature library can be directly used for subsequent training of cross-modal fusion models (such as attention mechanism fusion networks) and dynamic evolution analysis of pathogenic bacteria.

[0017] Step S120: Introduce a dual-branch attention mechanism to strengthen intra-modal features and establish bidirectional inter-modal correlations, and then weightedly fuse them to generate a cross-modal feature set; A two-branch attention fusion model is built based on the PyTorch framework. The intra-modal branch uses a self-attention mechanism to strengthen core features, while the inter-modal branch establishes bidirectional associations through cross-attention. Specifically, genomic, transcriptomic, and proteomic features are input into the intra-modal self-attention module, outputting enhanced features. Subsequently, genomic and transcriptomic, genomic and proteomic, and transcriptomic and proteomic features are input pairwise into the cross-modal cross-attention module to generate bidirectional associated features. A gated recurrent unit (GRU) is introduced as an adaptive weight adjustment mechanism, inputting the information entropy (quantifying contribution) of each modal feature and the inter-modal mutual information value (quantifying the degree of association), dynamically learning the fusion ratio (intra-modal weight range 0.4-0.6, inter-modal weight range 0.4-0.6). Element-wise weighted summation is used during weighted fusion, combined with L2 regularization (λ=0.001) to suppress overfitting, ultimately generating a cross-modal feature set with dimensions of 200 samples × 4500 features. The model was trained using the Adam optimizer (learning rate 0.0005, batch size=32), and converged after 50 iterations, with a core information retention rate of ≥92% for the feature set.

[0018] Step S121: Perform dimensionality unification and normalization on the genomic, transcriptomic, and proteomic features in the multimodal original feature library, respectively. Use a feature mapping algorithm to transform the heterogeneous features of each modality to the same high-dimensional feature space. Use a standardization method to eliminate the numerical scale differences of each modality feature to obtain the standardized feature matrix corresponding to each modality. The original feature dimensions for each modality are as follows: 2100 dimensions for the genome, 1800 dimensions for the transcriptome, and 1500 dimensions for the proteome. A multilayer perceptron (MLP) is used as the feature mapping algorithm, constructing a 3-layer fully connected network (input layer dimension corresponding to the number of features for each modality, hidden layer dimension 2048, output layer dimension 2048) to uniformly map heterogeneous features to a 2048-dimensional high-dimensional space. Standardization is performed using the Z-score method, implemented through scikit-learn's StandardScaler tool. The mean μ and standard deviation σ of each modality feature are calculated and transformed according to the formula x'=(x-μ) / σ, ensuring that the mean ≈ 0 and the variance ≈ 1 after transformation. Outliers (deviations from the mean by more than 3 standard deviations) are truncated during processing, retaining only valid data. Finally, standardized feature matrices for the genome, transcriptome, and proteome are obtained, each with a dimension of 200 samples × 2048 features, and the numerical scale is uniformly within the [-3,3] interval, meeting the numerical stability requirements for subsequent attention calculations.

[0019] Step S122: For the standardized feature matrix of each modality, construct the intramodal self-attention calculation module, generate query, key, and value matrices through linear transformation, obtain the intramodal attention weight matrix based on similarity calculation, perform weighted summation with the value matrix after normalization, strengthen the key core features of the function and suppress redundant features, and output the intramodal enhanced feature matrix of each modality. A self-attention module for modalities was built using PyTorch, taking a standardized feature matrix of 200×2048 as input for each modality. This was achieved through a 3-layer linear transformation layer with shared parameters. features =2048,out featuresGenerate query (Q), key (K), and value (V) matrices (all 200×2048) using a matrix of √2048 (bias=True). Initialize layer parameters using a Xavier normal distribution. Set a scaling factor of √2048≈45.25 to eliminate bias in high-dimensional similarity calculations. Calculate the matrix multiplication of Q and the transpose of K to obtain the original similarity matrix (200×200). After dividing by the scaling factor, generate an attention weight matrix using Softmax (dim=-1), with core feature weights ≥0.05 and redundant features ≤0.01. Multiply the weight matrix with the V matrix to obtain the preliminary enhanced feature matrix. Introduce a residual connection (add the original features to the preliminary enhanced features element-wise after linear transformation), and then normalize using LayerNorm (eps=1e-5). Finally, output the intra-modal enhanced feature matrix of 200×2048 for each modality. The core feature response intensity is improved by more than 35% compared to the original features.

[0020] Step S1221: For the single-modality normalized feature matrix, construct a three-layer linear transformation layer with shared parameters, which is used to generate query, key, and value matrices respectively. The parameters of the transformation layer are initialized using a parameter initialization method to ensure convergence. Set a scaling factor based on the feature dimension to eliminate the scale bias in similarity calculation under high-dimensional feature space. A three-layer linear transformation layer with shared parameters is constructed, implemented using PyTorch's nn.Linear. Each layer's parameters are set to in_features=2048, out_features=2048, and bias=True. The weights of all three layers are fully shared to reduce parameter redundancy. Xavier normal distribution is used to initialize parameters to ensure model training convergence. A scaling factor is set based on the feature dimension: since the standardized feature matrix has a dimension of 2048, the scaling factor is taken as the square root of the feature dimension, i.e., √2048≈45.25, to compensate for the numerical amplification problem caused by matrix multiplication in high-dimensional space. During training, the loss function (MSE loss) is monitored. If the loss decrease rate is <1% after 10 iterations, the parameters are reinitialized and the scaling factor is fine-tuned (±5%) to ensure the stability and accuracy of attention calculation.

[0021] Step S1222: Input the standardized feature matrix into the three linear transformation layers respectively, and generate a query matrix, key matrix and value matrix with uniform dimensions through matrix multiplication; perform dimension verification on the generated three types of matrices, remove feature rows with abnormal dimensions, and ensure that the matrix dimensions match and the data is complete; The standardized feature matrices of each mode (200×2048) are input into the three-layer shared linear transformation layer, and matrix multiplication is performed. Generate Q, K, and V matrices, each with an output dimension of 200×2048. Dimension validation employs NumPy's `shape` function and an outlier detection mechanism: First, verify that the row dimension of the Q, K, and V matrices is 200 and the column dimension is 2048. If the dimensions do not match, trace back to the feature mapping step for reprocessing. Second, perform dimension detection on each row of the matrix features. If the number of feature elements in a row is ≠ 2048, it is determined to be an outlier row. Set a threshold of ≤5% for outlier rows. If the threshold is not exceeded, the outlier row is directly removed; if the threshold is exceeded, standardization is re-executed. After validation, ensure that the dimensions of the Q, K, and V matrices are perfectly matched, the data integrity rate is ≥95%, and there are no missing or outlier dimensional features.

[0022] Step S1223: Using a scaled dot product attention mechanism, matrix multiplication is performed on the transpose of the query matrix and the key matrix to obtain the original feature similarity matrix; the original similarity matrix is ​​divided by the preset scaling factor to complete scale normalization, generating a scaled similarity matrix, and quantifying the association strength of the corresponding feature pairs. Feature similarity is calculated using a scaled dot product attention mechanism: First, the Q matrix (200×2048) is multiplied by the transpose of the K matrix (2048×200) to obtain a 200×200 original feature similarity matrix, with element values ​​ranging from -500 to 500. Each element in the original similarity matrix is ​​then divided by a preset scaling factor of 45.25 to normalize the scale, generating a scaled similarity matrix with element values ​​mapped to [-11.05, 11.05], effectively reducing numerical bias in high-dimensional space. The association strength of feature pairs is quantified: similarity values ​​≥ 2.0 are defined as strong association, 0.5-2.0 as moderate association, and < 0.5 as weak association. Subsequent weight allocation will prioritize strongly associated features to ensure the focus on core information within the modality.

[0023] Step S1224: Input the scaled similarity matrix into the Softmax activation function, normalize the elements of each row of the matrix, generate the intramodal attention weight matrix, and realize the differentiated weight allocation of core functional features and redundant features. The scaled 200×200 similarity matrix is ​​input into PyTorch's nn.Softmax(dim=-1) activation function. Each row of the matrix is ​​normalized so that the sum of all elements in each row is 1, generating an intra-modal attention weight matrix. Weight allocation follows the principle of high weight for core functional features and low weight for redundant features: differentiated allocation is achieved through preset thresholds. Weights ≥0.05 are assigned to core functional features related to virulence and metabolism, while weights ≤0.01 are assigned to redundant features derived from sequencing noise and repetitive sequences. For example, the weights corresponding to virulence gene domains in the genome are concentrated between 0.06 and 0.08, and the weights of low-expression redundant transcripts in the transcriptome are ≤0.008, ensuring that attention resources are tilted towards key features.

[0024] Step S1225: Perform matrix multiplication on the attention weight matrix and the value matrix to obtain the preliminary enhanced feature matrix; introduce the residual connection mechanism, perform element-wise addition on the original standardized feature matrix and the preliminary enhanced feature matrix after linear transformation, and then eliminate the numerical offset through layer normalization to output the intramodal enhanced feature matrix. A matrix multiplication operation is performed between the 200×200 attention weight matrix and the 200×2048 value matrix to obtain a preliminary 200×2048 enhanced feature matrix. A residual connection mechanism is introduced: a single linear transformation layer (in...) is constructed. features =2048,out features The original normalized feature matrix (normalized_shape=2048, eps=1e-5) is converted into dimension-matched residual features, which are then element-wise added to the initial enhanced feature matrix (torch.add(enhanced,residual)) to alleviate the vanishing gradient problem during deep network training. Subsequently, a LayerNorm layer (normalized_shape=2048, eps=1e-5) is used for normalization to eliminate numerical bias, ensuring that the mean of the enhanced features is approximately 0 and the variance is approximately 1. Finally, a 200×2048 intra-modal enhanced feature matrix is ​​output for each modality, improving the stability of the feature response by 40% compared to the original features.

[0025] Step S1226: Calculate the weight ratio and variance contribution rate of functional features in the enhanced feature matrix. If the weight increase of the core feature does not reach the preset threshold, fine-tune the linear transformation layer parameters or scaling factor, and repeat the above steps to screen high contribution dimensions in the enhanced feature matrix and remove low contribution redundant dimensions, and finally determine the intramodal enhanced feature matrix of each modality.

[0026] Calculate the functional feature weight proportions and variance contribution rates of the enhanced feature matrix: Statistical methods are used to calculate the sum of weights for the top 20% of core functional features, with a preset threshold of ≥60%; PCA analysis is used to calculate the feature variance contribution rate, with a preset cumulative contribution rate of ≥85% for the top 500 features. If the core feature weight increase does not reach the threshold (e.g., only 25% < 30%), the linear transformation layer parameters (learning rate adjusted to 0.0003) or scaling factor (adjusted to 42 or 48) are fine-tuned, and steps S1221-S1225 are re-executed. Simultaneously, low-contribution redundant dimensions are removed: features with a variance contribution rate <0.1% are defined as low-contribution dimensions, and after batch removal, the dimensions of each modality's enhanced feature matrix are optimized to 200×1500. Finally, the core feature weight proportions for the genomics, transcriptomics, and proteomics reach 68%, 72%, and 65%, respectively, and the variance contribution rates all meet the requirement of ≥85%.

[0027] Step S123: Construct a cross-modal cross-attention module. Take the enhancement feature matrices of any two modalities as input, generate mutual query matrix and mutual key matrix through linear transformation respectively, calculate the mutual information similarity between features of different modalities to generate a bidirectional cross-attention weight matrix, and perform weighted filtering on the enhancement features of another modality based on the weight matrix to obtain the inter-modal association feature matrix, realizing the bidirectional association mapping between genome and transcriptome, genome and proteome, and transcriptome and proteome. A cross-modal attention module is built using PyTorch, taking as input two 1500-dimensional enhancement feature matrices of any two modalities (e.g., genome G: 200×1500, transcriptome T: 200×1500). Two independent linear transformation layers are constructed (G, cross-interrogation / cross-key, T, cross-interrogation / cross-key) to generate Q. G K G Q T K T Matrix (all 200×1500), with a mutual information calculation window size set to 50. A similarity matrix (M) is generated through bidirectional mutual information calculation (I(G,T)=H(G)+H(T)-H(G,T)). GT : 200×200, M TG (200×200), and a bidirectional cross-attention weight matrix is ​​obtained through Softmax. The weight matrix is ​​then multiplied by the other modality enhancement feature matrix, and strongly correlated features (weight ≥ 0.03) are selected to obtain A. GT (200×1500), A TG A 200×1500 correlation feature matrix is ​​generated. Bidirectional correlation mappings are then performed sequentially for GT, GP (proteome), and TP to achieve accurate cross-modal feature correlation.

[0028] Step S1231: Construct two independent linear transformation layers corresponding to the bimodal, initialize the transformation layer parameters to ensure mapping stability, and set the mutual information calculation window to provide support for similarity measurement. For each set of bimodal inputs (e.g., GT), two independent linear transformation layers are constructed: the first set is used for mutual querying of genome enhancement feature generation (Q). G ) and interbonding (K) G The matrix, the second group is used for transcriptome enhancement feature generation mutual query (Q) matrix. T ) and interbonding (K) T ) matrix. The network parameters for each layer are set to in features =1500、out features =1500, bias=True, and Xavier uniform distribution initialization is used to ensure the stability of bimodal feature mapping. The mutual information calculation window is set to 50, meaning that the mutual information values ​​of 50 feature dimensions are calculated each time, accumulating to obtain the full-dimensional mutual information matrix. The window size is determined based on 1 / 30 of the feature dimension to balance computational efficiency and association accuracy. During training, the gradient norm of the weight matrix is ​​monitored to ensure that it is ≤1.0 to avoid gradient explosion.

[0029] Step S1232: Select the dual-modal enhanced feature matrix as the input pair, generate four sets of mutual query and mutual key matrices through linear transformation, and ensure that the query and the corresponding key matrix column dimensions match through dimensional consistency verification; Select a bimodal enhancement feature matrix as the input pair (e.g., G: 200×1500, T: 200×1500), and input them into the corresponding linear transformation layer respectively: G generates Q through the first set of transformation layers. G (200×1500) and K G (200×1500), T generates Q through the second set of transformation layers. T (200×1500) and K T (200×1500), a total of four sets of mutual query and mutual key matrices. Perform dimension consistency verification: via torch.equal(Q) G .shape,K T .shape) and torch.equal(Q T .shape,K G The `.shape` function verifies that all column dimensions are 1500. If a dimension mismatch exists (e.g., a modality has a dimension of 1480 due to feature selection), zero-padding is used to fill the lower-dimensional matrix to 1500 dimensions. After verification, it ensures that the dimensions of the four sets of matrices are completely matched, data integrity is ≥98%, and there are no dimension conflicts.

[0030] Step S1233: Based on the generated mutual query and mutual key matrix, construct a two-way mutual information calculation process, integrate the original association matrix and mutual information values ​​to generate a two-way mutual information similarity matrix, and quantify the association strength between the two modalities; Firstly, based on Q G With K T Calculate the original incidence matrix (Q) G @K T The mutual information value I(G,T) between G and T (range 0-1) is calculated using kernel density estimation. The elements of the original association matrix are multiplied by the corresponding mutual information values ​​to generate the bidirectional mutual information similarity matrix M between G and T. GT (200×200). Similarly, based on Q... T With K G The similarity matrix M between T and G is calculated. TG (200×200). Quantify the correlation strength between the two modes: Set M GT and M TG A correlation strength of ≥0.3 indicates a strong association, 0.1-0.3 indicates a moderate association, and <0.1 indicates a weak association. For example, the association strength between genomic virulence genes and differentially expressed transcripts in the transcriptome is ≥0.35, ensuring accurate matching of core features across modalities.

[0031] Step S1234: Input the bidirectional mutual information similarity matrix into the Softmax activation function for normalization to generate a bidirectional cross-attention weight matrix that reflects the contribution of the source mode to the target mode. M GT (200×200) and M TG The (200×200) bidirectional mutual information similarity matrix is ​​input into the nn.Softmax(dim=-1) activation function and row normalized to generate the bidirectional cross-attention weight matrix W. GT (G's correlation contribution to T) and W TG (The correlation contribution of T to G), the sum of the elements in each row of the matrix is ​​1. The weight matrix quantifies the correlation contribution of the source mode to the target mode: for example, W GT In the study, the weight of transcriptome functional features corresponding to genomic virulence gene characteristics was ≥0.04, while the weight of non-associated features was ≤0.005; W TG In this study, the weight of transcriptome metabolic pathway-related features corresponding to genomic gene sequence features is ≥0.035, ensuring that the contribution allocation is consistent with the biological function and avoiding interference from the weight of meaningless features.

[0032] Step S1235: Perform matrix multiplication operations between the bidirectional attention weight matrix and the other modality enhancement feature matrix respectively, filter strongly correlated features to obtain two sets of correlated feature matrices, which together constitute the intermodal correlated feature matrix of the bimodal pair; Feature weighting selection using a bidirectional attention weight matrix: W GT Performing matrix multiplication between the (200×200) transcriptome enhancement feature matrix T (200×1500) and the Transcriptome Enhancement Feature Matrix A (200×1500) yields the correlation feature matrix A of G→T. GT (200×1500); W TG Multiplying (200×200) by the genome enhancement feature matrix G (200×1500) yields the association feature matrix A of T→G. TG (200×1500). The two sets of matrices together constitute the intermodal correlation feature matrix (200×3000) of the GT bimodal pair. The selection criteria are to retain strong correlation features with a weight ≥0.03 and remove weak correlation features (weight <0.03). After processing, A GT and A TG The proportions of medium and strong correlation features are ≥60% and 58% respectively, effectively extracting key cross-modal correlation information.

[0033] Step S1236: Select the genome and transcriptome, genome and proteome, and transcriptome and proteome as bimodal input pairs in sequence, repeat the above steps to generate bidirectional association feature matrices for each group, and complete the bidirectional association mapping of the full cross-modal combination.

[0034] Following the order of genome, transcriptome → genome, proteome → transcriptome, proteome, three bimodal input pairs were selected sequentially, and steps S1231-S1235 were repeated. When processing the genome-proteome (GP), the inputs were G (200×1500) and P (200×1500), generating A. GP (G→P) and A PG (P→G) Association Feature Matrix; When processing transcriptome-proteome (TP), the inputs are T (200×1500) and P (200×1500), and the generated matrix is ​​A. TP (T→P) and A PT (P→T) Correlation Feature Matrix. All bimodal pairs maintain consistent processing parameters (linear transformation layer dimension, scaling factor, weight threshold, etc.) to ensure comparability of results. Ultimately, six sets of correlation feature matrices (200×1500 each) are obtained, completing the bidirectional correlation mapping across all modal combinations and covering the feature relationships between all modes.

[0035] Step S124: Introduce an adaptive weight adjustment mechanism based on gating units. Based on the information contribution of each modal feature and the degree of correlation between modalities, dynamically learn the fusion ratio of intramodal self-attention weights and intermodal cross-attention weights. Perform secondary normalization on the adjusted weight matrix to ensure the rationality of weight allocation and the balance of feature contributions. An adaptive weight adjustment mechanism is constructed by introducing an MLP-based gating unit. The gating unit consists of one hidden layer (128 dimensions, sigmoid activation function) and one output layer (1 dimension). Input features include the information contribution of each modality (calculated through feature entropy values; higher entropy values ​​indicate greater contribution) and the degree of correlation between modalities (mutual information ≥ 0.5 indicates close correlation). The output is a fusion ratio coefficient α (range 0-1), where α represents the proportion of self-attention weights within a modality, and 1-α represents the proportion of cross-attention weights between modalities. For example, when the genomic feature information contribution is 0.7 and the mutual information with the transcriptome is 0.65, α = 0.55, and the inter-modal weight is 0.45; when the proteome contribution is 0.6 and the mutual information with the genome is 0.52, α = 0.5, and the inter-modal weight is 0.5. The adjusted weight matrix is ​​L2 normalized (torch.norm(weight,p=2,dim=-1)) to ensure the rationality of weight allocation and the balance of feature contributions, and to avoid the weight of a single mode exceeding 70%.

[0036] Step S125: The intra-modal enhancement feature matrix and the corresponding inter-modal correlation feature matrix of each modality are summed element-wise according to the adaptively adjusted weights to obtain the preliminary fusion feature matrix of each modality; the preliminary fusion feature matrices of all modalities are concatenated to form a high-dimensional cross-modal candidate feature set; the risk of overfitting is reduced by regularization, and the cross-modal fusion feature set is finally output.

[0037] Preliminary fusion characteristics of the genome = α G ×Genomic enhancement features + (1-α) G )×(A GT +A GP ) / 2, transcriptome = α T × Transcriptome enhancement features + (1-α) T )×(A TG +A TP ) / 2, proteome = α P ×Proteome enhancement features + (1-α) P )×(A PG +A PT ) / 2, where α G α T α PThe fusion coefficients output by the gating unit are all in the range of 0.4-0.6. The three preliminary fusion feature matrices (all 200×1500) are concatenated column-wise to form a high-dimensional cross-modal candidate feature set of 200×4500. Dropout regularization (p=0.3) is introduced to reduce the risk of overfitting, while L2 regularization (λ=0.001) is used to constrain the parameter size. The final output cross-modal fusion feature set has been validated, showing feature redundancy ≤15% and core information coverage ≥93%, and can be directly used for subsequent tasks such as pathogen virulence prediction and drug resistance analysis.

[0038] Step S130: The pathogen invasion process is abstracted into a temporal interaction network. A dynamic framework is built by embedding gated recurrent units into a temporal graph convolutional network, and the edge weights are dynamically adjusted to characterize the network evolution. A dynamic evolutionary analysis framework based on PyTorch+PyTorch Geometric is built, abstracting the process of pathogenic bacteria (such as Staphylococcus aureus) invading the host as a temporal interaction network. Four time steps are divided according to the invasion biological stages (t0: adhesion phase 0h, t1: invasion phase 6h, t2: proliferation phase 12h, t3: diffusion phase 24h). A temporal graph convolutional network (T-GCN) is used to extract spatial features, and gated recurrent units (GRUs) are embedded to capture temporal dependencies. The core design consists of a graph convolutional layer with an output dimension of 1024, a GRU hidden layer with a dimension of 1024, and an edge weight update function that integrates node feature changes with biological interaction patterns. The framework is trained using the Adam optimizer (learning rate 0.0003, batch size). size =16), converged after 80 iterations, with a dynamic edge weight adjustment error ≤5%. This framework can characterize the interaction strength between pathogenic virulence genes and host immune genes over time, and accurately capture key invasion nodes (such as the interaction activation between t1 phase adhesion proteins and host receptors).

[0039] Step S131: Divide the time steps according to the pathogen invasion stage, define the types of network nodes and edges, collect the interaction data of each time step to construct a static network snapshot, arrange them in time sequence to form an initial temporal interaction network and label the relevant attributes; The time steps were divided according to the cellular experimental time points of pathogen invasion: t0 (0h, no invasion), t1 (6h, mucosal adhesion), t2 (12h, cell invasion), and t3 (24h, intracellular proliferation). Network node types were defined: pathogen nodes (80 core virulence genes, such as spa and agrA) and host nodes (120 immune-related genes, such as TLR2 and IL-6), totaling 200 nodes; edge types included protein-protein interaction edges, gene regulation edges, and metabolic pathway interaction edges. Data were collected at each time step: transcriptome RNA-seq expression levels, Co-IP-verified interaction relationships, and pathway activation status. Four static network snapshots were constructed (each snapshot containing a 200×200 adjacency matrix + a 200×4500 node feature matrix). The initial temporal interaction network is formed by arranging nodes in the order t0→t1→t2→t3, and the node attributes (gene ID, functional annotation, expression level) and edge attributes (interaction type, initial strength) are labeled to ensure that the biological scene and time step of each snapshot are accurately matched.

[0040] Step S132: Initialize nodes with cross-modal fusion of pathogenic and host features, assign initial weights to edges based on experimental data and literature information, and construct a time-step correlation matrix to quantify the correlation degree of node states; Nodes were initialized using the cross-modal fusion feature set output by S125: each of the 200 nodes was assigned a 4500-dimensional feature vector (including genomic virulence domains, transcriptomic differential expression, and proteomic abundance information). Initial edge weighting rules were applied: protein-protein interaction edges referenced STRING database interaction scores (normalized to 0-1); experimentally validated interactions (e.g., Co-IP positivity) were weighted at 0.8-1.0, and predictive interactions at 0.3-0.7. A time-step correlation matrix was constructed: the Pearson correlation coefficients of the same node at adjacent time steps were calculated, forming a 200×200 correlation matrix with elements ranging from -1 to 1, quantifying the temporal correlation of node states (e.g., the correlation between t0-t1 adhesion gene nodes ≥ 0.7). The correlation matrix, after Z-score normalization, served as auxiliary input for GRU temporal dependency calculation, ensuring the accuracy of node state correlation characterization.

[0041] Step S133: Construct a graph convolution module to extract network snapshot spatial features, embed gated recurrent units to integrate temporal dependency features, form a spatiotemporal fusion structure, and design a forward propagation path; A spatiotemporal fusion structure based on graph convolution and GRU is constructed: The graph convolution module uses the Kipf first-order approximation GCN operator, with the weight matrix dimension set to 4500×1024 and initialized using a Xavier normal distribution; the adjacency matrix processing function is A'=(A+I) / D (A is the original adjacency matrix, I is the identity matrix, and D is the degree matrix), eliminating node degree difference bias; the activation function is ReLU, and the dropout rate is set to 0.3. Each time step snapshot is decomposed into a 200×4500 node feature matrix and a 200×200 normalized adjacency matrix, input to the graph convolutional layer, and outputting a 1024-dimensional spatial feature sequence. The GRU hidden layer dimension is matched to 1024, and the reset and update gate weights are orthogonally initialized. The spatial features are converted to a 1×1024 temporal input format through an adaptation layer. A serial fusion structure is constructed, introducing residual connections (linear transformation 1024→1024) and LayerNorm, and adding a fully connected layer to map back to 4500 dimensions. The forward propagation path follows a sequence of snapshot input → graph convolution extraction → format conversion → GRU calculation → fusion output. After verification with 50 sample snapshots, the dimension matching rate is 100%, and the computational efficiency is improved by 30%.

[0042] Step S1331: Select the first-order approximate graph convolution as the operator. Define the graph convolution layer weight matrix based on the node feature dimension of the temporal interaction network snapshot and the preset spatial feature output dimension. Initialize the weight matrix using the parameter initialization method. Construct a self-loop correction adjacency matrix processing function to add an identity matrix to the original network adjacency matrix and perform normalization processing to eliminate feature extraction bias caused by node degree differences. Set the activation function and dropout rate of the graph convolution module to reduce the risk of overfitting. The Kipf first-order approximate graph convolution is selected as the core operator, and the weight matrix of the graph convolution layer is defined. (Input dimension is 4500 for node features, output dimension is 1024 for spatial features), initialized using a Xavier normal distribution (nn.init.xavier_normal). (W) ), bias term Initialize to 0.01. Construct the adjacency matrix processing function with self-loop correction: First, add the identity matrix I to the original adjacency matrix A (200×200) to obtain A. self =A+I, then calculate the degree matrix D (the diagonal elements are A). self Each row and), through Normalization is completed. The activation function is set to ReLU (inplace=True), and the dropout rate is 0.3 (nn.Dropout(0.3)). During training, some neurons are randomly deactivated to reduce the risk of overfitting. After parameter initialization, the gradient norm is checked (≤1.0) to ensure stable convergence of subsequent training.

[0043] Step S1332: Decompose the network snapshot of a single time step into a node feature matrix and a normalized adjacency matrix. Input the node feature matrix into the graph convolutional layer to perform matrix operations to complete spatial feature extraction and output the spatial feature matrix of that time step. Perform the above operation on the network snapshot of all time steps in sequence to obtain the spatial feature sequence, and perform dimension verification on the feature sequence to ensure that the output feature dimension of each time step is consistent. A single time-step network snapshot is decomposed into two parts: the node feature matrix. (200 nodes, each with 4500 cross-modal features), normalized adjacency matrix Input X into the graph convolutional layer, and apply the formula H=ReLU(A) norm Perform matrix operations (×X×W+b) and output the spatial feature matrix for that time step. The four snapshots are processed sequentially according to the time step order (t0→t1→t2→t3) to obtain a spatial feature sequence of length 4 [H0,H1,H2,H3]. Dimension verification uses the numpy.shape function to verify that the dimension of each Hi is (200,1024). If an abnormal dimension is found (e.g., H1 is (198,1024)), the missing node features are filled in by tracing back to the node filtering step to ensure that the output feature dimensions of all time steps are completely consistent, and the data integrity rate is ≥99%.

[0044] Step S1333: Set the GRU hidden layer dimension based on the spatial feature output dimension, initialize the GRU reset gate, update gate, weight matrix and bias term corresponding to the candidate hidden state; design the GRU input adaptation layer to convert the spatial feature matrix of each time step into a format that matches the GRU time sequence input requirements, forming the GRU input sequence. Based on the spatial feature dimension (1024) of the graph convolution output, the GRU hidden layer dimension is set to 1024 to ensure feature dimension matching. Initialize GRU core parameters: reset gate weights. Update gate weights Candidate hidden state weights All parameters are orthogonally initialized, with bias terms br, bu, and bc initialized to 0. The GRU input adaptation layer is designed as follows: a single linear transformation layer (in...) is constructed. features =1024,out features =1024), the spatial feature matrix of each time step Sum and compress into a vector by node dimension. The GRU input sequence is formed by concatenating the data in time step order. It meets the GRU's format requirements for time-series input.

[0045] Step S1334: Connect the output of the graph convolution module to the input of the GRU module to construct a serial fusion structure of "graph convolution + GRU"; introduce a residual connection mechanism to linearly transform the spatial features output by the graph convolution and add them element-wise to the hidden state of the current time step output by the GRU, and then perform layer normalization; add a fully connected layer after the fusion structure to map the dimension of the hidden state output by the GRU back to the dimension that matches the node features, forming a complete spatiotemporal fusion feature extraction module; The output of the graph convolution module is directly connected to the input of the GRU module to construct a serial fusion structure where graph convolution extracts spatial features → GRU integrates temporal features. A residual connection mechanism is introduced: a linear transformation layer (in...) is built. features =1024,out features =1024), transforming the spatial feature vector v output by graph convolution into a residual vector r, and combining it with the hidden state h at the current time step output by GRU. t Perform element-wise addition (h) t' =h t+r This alleviates the vanishing gradient problem. The summation result is passed through a LayerNorm (normalized) layer. shape =1024, eps=1e-5) normalization to eliminate numerical offset. A fully connected layer (in) is added after the fused structure. features =1024,out features =4500), mapping the 1024-dimensional hidden state output by GRU back to the 4500-dimensional features that match the original node features, forming a complete spatiotemporal fusion feature extraction module that combines spatial extraction, temporal integration, and dimensional restoration.

[0046] Step S1335: Define the complete data flow path of the forward propagation, which is as follows: the input layer receives network snapshots at each time step, the graph convolutional layer extracts spatial features, the adaptation layer converts the input format, the GRU layer calculates and outputs the hidden state in the order of time steps, and the fusion output layer outputs the spatiotemporal fusion features; define the computation priority and data dependency of each layer to ensure that the data flows unidirectionally in the order of time steps and has no reverse dependency; The complete forward propagation data flow path is defined as follows: 1) The input layer receives network snapshots at four time steps (A0, X0) → (A3, X3); 2) The graph convolutional layer processes the snapshots sequentially at each time step, outputting the spatial feature sequence [H0, H1, H2, H3]; 3) The adaptation layer compresses each Hi into a vector vi, forming the GRU input sequence V; 4) The GRU layer calculates in the order t0 → t1 → t2 → t3, receiving vi and updating the hidden state at each step, outputting the hidden state sequence [h0, h1, h2, h3]; 5) The fusion output layer outputs the spatiotemporal fusion feature sequence [F0, F1, F2, F3] (Fi ∈ ℝ) through residual connections and fully connected layers. (200×4500)Define the computation priority of each layer: input layer > graph convolutional layer > adaptation layer > GRU layer > fusion output layer, clarify the data dependencies (GRU computation depends on the hidden state of the previous time step and the current time step vi, with no reverse dependency), ensure that data flows unidirectionally according to time steps, and avoid computation conflicts.

[0047] Step S1336: Select a small number of temporal network snapshots and input them into the fusion module. Verify the feature dimension matching step by step. If dimension incompatibility issues occur, adjust the output dimension of the graph convolution or the dimension of the GRU hidden layer and reinitialize the weight parameters. Optimize the matrix operation method to improve computational efficiency and ensure that the module can stably output the spatiotemporal fusion feature sequence.

[0048] 160 time-series network snapshots (4 time steps per sample) from 20 samples were selected and input into the spatiotemporal fusion module for step-by-step verification: First, the output Hi dimension of the graph convolutional layer was verified (all must be 200×1024), then the output vi dimension of the adaptation layer was verified (must be 1×1024), and finally the output hi dimension of the GRU was verified (must be 1×1024). If dimension incompatibility occurred (e.g., GRU output was 1×512), the graph convolutional output dimension was adjusted to 512 and the weight parameters were reinitialized. Matrix operation method was optimized: PyTorch's batch matrix multiplication (torch.bmm) was used to replace loop operations, reducing the single-sample processing time from 0.8s to 0.2s. After three rounds of verification and optimization, the module can stably output spatiotemporal fusion feature sequences with 100% feature dimension consistency, and the computational efficiency meets the requirements for large-scale time-series data processing.

[0049] Step S134: Construct an edge weight update function based on the temporal characteristics of nodes, set a differentiated adjustment threshold in combination with biological interaction rules, and achieve adaptive update of edge weights as the invasion process is achieved by adjusting the coefficients. An edge weight update function is constructed based on the temporal feature changes of nodes: The spatiotemporal fusion feature vectors of the same node in adjacent time steps (e.g., vi for t0 and t1) are extracted, and the feature difference is calculated using cosine similarity (normalized to 0-1). A dual-input update function is constructed by combining GRU temporal association weights (outputting the weight coefficients of the hidden state). Edges are categorized into three types according to biological event type: pathogen-host interaction edges, pathogen intra-pathogen regulatory edges, and host intra-pathway edges. Differentiation thresholds are set using the KEGG pathway model (e.g., a pathogen-host interaction edge trigger threshold of 0.25, with weight adjustment ranges from 0.1 to 0.95). The network edges at the current time step are traversed, and the basic adjustment coefficient (β = 0.3 × difference + 0.7 × temporal weight) is calculated. After biological rule correction (e.g., a correction coefficient of 1.2 for virulence gene interaction edges) and truncation, the final adjustment coefficient is obtained. The initial edge weight is multiplied by the adjustment coefficient to obtain the dynamic edge weight, which is then updated synchronously to the corresponding snapshot after 3σ outlier removal and 3-step moving average filtering. After 10 rounds of iterative optimization, the weight update results are ≥88% consistent with the experimentally measured interaction strength.

[0050] Step S1341: Extract the spatiotemporal fusion feature vectors of the same node in adjacent time steps, and obtain standardized node temporal feature difference values ​​through feature difference degree calculation and normalization. Extracting adjacent time steps (e.g., t) i With t {i+1} Spatiotemporal fusion feature vectors of the same node and The feature difference is calculated using cosine similarity: The degree of difference d = 1 - cos sim (Value range 0-1). Perform Min-Max normalization on the dissimilarity of all nodes: d norm =(dd min ) / (d max -d min This ensures that the standardized node temporal feature differences all fall within the [0,1] interval. For example, the original feature difference value of the pathogenic bacteria adhesion protein node in the t1-t2 phase was 0.35, which, after normalization, became d... norm =0.62, which intuitively reflects the intensity of the functional state change of the node during the invasion and proliferation stage, and provides a quantitative basis for edge weight adjustment.

[0051] Step S1342: Construct an edge weight update function with node feature difference degree and GRU temporal association weight as dual inputs. The function is designed in two layers: basic adjustment coefficient calculation and biological rule correction, to ensure that the coefficients conform to the laws of biological interaction dynamics. Construct the edge weight update function as w new =w old ×(1+α×d norm +β×wgru ), where α is the feature difference weight coefficient (set to 0.4), β is the temporal correlation weight coefficient (set to 0.6), and w gru The time-series correlation weights output by the GRU (taken from the normalized coefficients of the hidden states). The function is designed in two layers: the first layer calculates the basic adjustment coefficient k. base =α×d norm +β×w gru The second layer involves biological rule correction: referencing the KEGG pathogen-host interaction pathway model, a correction coefficient of 1.2 is added to the virulence gene-host receptor interaction side, and a correction coefficient of 0.9 is added to the non-core regulatory side, ensuring that the adjustment coefficients conform to the laws of biological interaction dynamics (e.g., the interaction intensity increases faster during the critical period of invasion). By comparing with experimental interaction data, the prediction accuracy of the corrected adjustment coefficients is improved by 23% compared to the basic model.

[0052] Step S1343: Classify network edges according to biological event types, and combine experimental data and pathway models to set differentiated state change trigger thresholds and upper and lower limits for weight adjustment for different types of edges. Based on the type of biological event, network edges are divided into three categories: Type 1 (pathogen virulence protein-host cell receptor interaction edge), Type 2 (pathogen endovirulence gene regulation edge), and Type 3 (host immune pathway gene interaction edge). Combining in vitro infection experimental data and the Reactome pathway model, differentiated parameters were set for different edge types: Type 1 trigger threshold 0.25 (weight adjustment is initiated if the feature difference exceeds this value), with upper and lower limits of weight adjustment ranging from 0.1 to 0.95; Type 2 trigger threshold 0.18, with upper and lower limits of weight adjustment ranging from 0.2 to 0.9; Type 3 trigger threshold 0.22, with upper and lower limits of weight adjustment ranging from 0.15 to 0.92. The rationale for these settings is as follows: Type 1 is a key interaction for invasion and requires a more sensitive adjustment response; Type 2 has a more stable regulatory relationship and a higher threshold; Type 3 is affected by dynamic changes in the immune system, with thresholds falling between the two, ensuring that the parameter settings closely reflect biological realities.

[0053] Step S1344: Traverse the network edges at the current time step, calculate the basic adjustment coefficient and determine whether it exceeds the trigger threshold. Obtain the final adjustment coefficient through biological rule correction and truncation processing. Traverse all edges (1860 in total) in the network at the current time step (e.g., t2), and process each edge as follows: 1) Extract the d_norm and w_gru of the corresponding node pair, and calculate the basic adjustment coefficient k. base =0.4×d norm +0.6×w gru ;2) Determine k base1) Does it exceed the trigger threshold for this type of edge (e.g., for type 1 edges, check if it is ≥0.25); 2) If it exceeds the threshold, add a correction coefficient according to the biological rules (add 1.2 for type 1, 0.9 for type 2, and 1.05 for type 3) to obtain k. temp =k base × Correction factor; 4) For k temp Perform truncation to ensure it falls within the [0.05, 0.5] interval (to avoid excessive weight changes), and obtain the final adjustment coefficient k. final For example, k of a type 1 edge. base =0.3, k after correction temp =0.36, not exceeding the upper limit, final k final =0.36.

[0054] Step S1345: The dynamic edge weights are obtained by calculating the initial edge weights and the final adjustment coefficients. After outlier removal and smoothing filtering, they are synchronously updated to the network snapshot at the corresponding time step. According to formula w new =w old ×(1+k final ) Calculate the dynamic edge weights, where w old This represents the initial weight of the current edge (or the weight updated in the previous time step). The calculated w... new Outlier removal: Using the 3σ principle, weight values ​​exceeding the range [μ-3σ, μ+3σ] (where μ is the mean of all edge weights at the current time step, and σ is the standard deviation) are removed, with the outlier percentage controlled to ≤3%. A 3-step moving average filter is then applied to the weights after outlier removal: w smooth =(w new [i-1]+w new [i]+w new [i+1]) / 3 (bilateral filtering is used for boundary edges), and the smoothed weight fluctuation is ≤10%. The smoothed dynamic edge weights are synchronously updated to the adjacency matrix of the network snapshot at the current time step to ensure that the snapshot can reflect the temporal changes of edge weights in real time and provide accurate input for feature extraction in the next time step.

[0055] Step S1346: Use measured data to compare and evaluate the accuracy of weight updates, adjust relevant parameters in reverse, and form a stable adaptive edge weight update mechanism through multiple rounds of iterative optimization.

[0056] In vitro infection data from 10 pathogenic bacteria (Co-IP interaction strength detection results for each bacteria at 4 time steps) were collected and compared with the weight update results as the experimental data. The Pearson correlation coefficient between the updated edge weights and the measured interaction strength was calculated. If r < 0.8, the parameters were adjusted in reverse (e.g., increasing α to 0.5 and decreasing the type 1 edge trigger threshold to 0.22). Through multiple rounds of iterative optimization: in the first round r = 0.72, after adjusting α and the threshold, in the third round r = 0.85, and in the fifth round r = 0.89, a stable adaptive edge weight update mechanism was finally formed. The optimized mechanism still had a correlation coefficient ≥ 0.87 on the independent test set (5 pathogenic bacteria), indicating that it has good generalization ability and can accurately capture the dynamic changes in biological interactions.

[0057] Step S135: Input network snapshots into the fusion framework according to time steps, output dynamic edge weights and node states, record evolution trajectories to generate evolution graphs, and highlight key evolution nodes and stages through dimensionality reduction visualization; Four network snapshots were sequentially input into the spatiotemporal fusion framework in time step order (t0→t1→t2→t3). The framework output dynamic edge weight matrices (one 200×200 matrix per step) and node state features (one 200×4500 matrix per step) per time step. The core outputs for each time step were recorded: edge weight changes, node feature response intensity, and key interaction pairs (edges with the top 20 weights), generating pathogen invasion evolution trajectory data (CSV format, including time step, node ID, edge ID, weight / feature value). The t-SNE dimensionality reduction algorithm was used to reduce the 4500-dimensional node features to 2-dimensionality. A time-series evolution map was plotted using Matplotlib, with node size representing feature response intensity, edge thickness representing dynamic weight, and color indicating node type. The visualization clearly highlights key evolutionary nodes (such as adhesion protein genes with high response in t1) and core stages (t1→t2 being the critical transition period from invasion to proliferation), providing intuitive support for biological mechanism analysis.

[0058] Step S136: Collect experimental data and divide it into training and validation sets. Optimize the framework parameters using a specific loss function. Verify the effectiveness of the core modules through ablation experiments. Determine the optimal parameters to ensure the accuracy of evolutionary characterization and biological significance.

[0059] Temporal interaction data of 100 pathogenic bacteria were collected (4 time steps per strain, totaling 400 network snapshots), and divided into a training set (70 strains) and a validation set (30 strains) in a 7:3 ratio. The loss function selected was MSE loss and a biological rule regularization term: MSE loss quantifies the difference between predicted edge weights and measured values, while the regularization term (λ=0.001) constrains weight adjustment to conform to the pathway model. The Adam optimizer (learning rate 0.0003, weight decay 1e-5) was used to optimize the framework parameters. After 100 iterations, the training set loss converged to 0.021, and the validation set loss to 0.028. Ablation experiments were designed to verify the effectiveness of the core modules: 1) Removing the GRU module (retaining only graph convolution) increased the validation set loss to 0.045; 2) Removing the edge weight update module (fixed weights) increased the loss to 0.052, demonstrating the core role of GRU temporal integration and dynamic weight update. The optimal parameters were finally determined as follows: graph convolution output dimension 1024, GRU hidden layer 1024, α=0.4, β=0.6. Under this parameter combination, the evolutionary characterization accuracy of the framework is ≥91%, and the biological significance is ≥89%.

[0060] Step S140: Extract the expression time-series curve of pathogenic bacteria virulence genes and dynamic data of host immune response, supplement them into the temporal interaction network, highlight the contribution of time nodes through temporal attention, and realize the dynamic update of node status; A temporal attention fusion framework was built based on PyTorch, extracting qPCR expression time-series curves of 80 Staphylococcus aureus virulence genes (hla, psmα, agrA, etc.) and dynamic data of host (human alveolar epithelial cells) immune response (concentrations of cytokines such as TNF-α and IFN-γ detected by ELISA). These two types of data were supplemented into the S135 temporal interaction network to construct a temporal attention layer highlighting the contribution of key time nodes (e.g., the 6-hour invasion period). Specifically, virulence gene temporal data were mapped to pathogen nodes by node ID, and immune data were mapped to host nodes. The attention layer quantified the contribution of data at each time step, generating weight coefficients (0-1), which were then weighted and fused with the spatiotemporal fusion features of the nodes. The model was trained using the Adam optimizer (learning rate 0.0004), converging after 60 iterations. After node state updates, the agreement between the model and the protein expression trends validated by Western blot was ≥90%, successfully capturing the dynamic association between virulence gene activation and immune response.

[0061] Step S141: Collect time-series data of pathogenic bacteria virulence gene expression and dynamic data of host immune response; perform noise reduction and normalization on the raw data, align the two types of data according to the preset time step of pathogenic bacteria invasion, fill in the data of missing time steps, and generate a standardized time-series data sequence that corresponds one-to-one with the time step of the time-series interaction network. The expression levels of 80 virulence genes in 100 pathogenic bacteria (50 MRSA strains and 50 MSSA strains) at four time steps (0h, 6h, 12h, and 24h) were determined by qPCR, and the concentrations of 10 immune factors in the host cell supernatant at the corresponding time steps were detected by ELISA. Noise reduction was achieved using Savitzky-Golay filtering (window size 5, polynomial order 2) to eliminate detection noise. Normalization was performed using the Z-score method (x'=(x-μ) / σ) to ensure that the mean ≈ 0 and the variance ≈ 1. The two types of data were aligned according to the invasion time step, and missing time steps (e.g., missing 6h IFN-γ data for some samples) were filled using the k-nearest neighbor algorithm (k=3) to ensure data integrity at each time step. Finally, a standardized time-series data sequence (100 samples × 90 indicators × 4 time steps) was generated, corresponding one-to-one with the four time steps of the time-series interaction network, with data integrity ≥ 98%, meeting the requirements for dynamic node updates.

[0062] Step S142: Construct a node and dynamic data association dictionary, map the time-series data of virulence gene expression to the pathogenic virulence factor nodes in the time-series interaction network, and map the dynamic data of host immune response to the corresponding host-related nodes; bind a unique time-series dynamic data sequence to each node and generate a three-dimensional data table; verify the uniqueness and completeness of the mapping relationship, remove redundant dynamic data, and ensure that the core nodes are supported by corresponding time-series dynamic features. The key is the node ID (e.g., pathogen node VF). hla Host Node IM TNFα The value corresponds to the dynamic data identifier (e.g., qPCR). hla ELISA TNFα To ensure uniqueness of the mapping, virulence gene expression time-series data (80 genes × 4 time steps × 1 expression level value) were mapped to 80 pathogenic bacterial virulence factor nodes, and host immune response data (10 factors × 4 time steps × 1 concentration value) were mapped to corresponding functional nodes in 120 host nodes (e.g., TNF-α data mapped to the host immune regulation node). A three-dimensional data table (100 samples × 200 nodes × 4 time steps × 1-dimensional dynamic data) was generated, and the mapping relationship was verified using a Python dictionary: no duplicate mappings (100% uniqueness), 100% dynamic data coverage of core nodes (Top 50 high-contribution nodes), and removal of 3 redundant dynamic data (metabolite concentrations unrelated to node function), ensuring that each core node has dedicated time-series dynamic features.

[0063] Step S143: Design a temporal attention layer based on the contribution of time nodes, set the attention score calculation dimension and output dimension, initialize the parameters of the weight matrix and bias term of the attention layer; define the core logical framework of attention calculation, and clarify the input and output content; set the dropout rate of the attention module to reduce the risk of overfitting; Design a temporal attention layer based on the contribution of time nodes, implemented using PyTorch's nn.Module: The input dimension is dynamic data feature dimension (1) + node spatiotemporal fusion feature dimension (4500) = 4501, and the output dimension is 1 (attention weights). Attention layer weight matrix Bias terms It uses Xavier normal distribution for initialization; the core logical framework is input features → linear transformation → tanh activation → weight mapping → score output, i.e. , (w is the output weight, c is the bias). The dropout rate is set to 0.3 (nn.Dropout(0.3)), randomly deactivating some neurons during training to reduce the risk of overfitting. Gradient clipping (max... norm =1.0) ensures stable training, and the weight values ​​output by the attention layer are constrained to the [0,1] interval, which conforms to the contribution quantification logic.

[0064] Step S144: Extract relevant features of dynamic data at each time step node, quantify the contribution of a single time step node to the overall state of the node; calculate the attention score using an attention mechanism, and obtain the original attention score after processing the contribution features of time steps with the spatiotemporal fusion features of the node; normalize the original scores of all time steps to generate the temporal attention weights of each time step, highlighting the role of high-contribution time steps. Three types of features were extracted from virulence gene expression data: peak value, slope of change, and mean expression level. Two types of features were extracted from immune response data: peak concentration and rate of increase. A total of five types of time node contribution features were quantified. An additive attention mechanism was used to calculate the original attention score: the five types of contribution features were concatenated with the node spatiotemporal fusion features (4500 dimensions) to form a 4505-dimensional vector, which was then input into the attention layer to obtain the original score. Softmax normalization is performed on the raw scores of the four time steps to generate temporal attention weights α=[α0,α1,α2,α3], satisfying α0+α1+α2+α3=1. High-contribution time nodes are highlighted: for example, the weight α1 of the pathogen invasion period (6h, t1) is generally between 0.35 and 0.45, significantly higher than other time steps (α0≈0.15, α2≈0.25, α3≈0.15), accurately focusing on the contribution of key invasion time nodes.

[0065] Step S145: Perform a weighted operation on the temporal attention weights of each time step and the standardized dynamic data of the corresponding node to obtain the attention-weighted node dynamic feature sequence; concatenate the sequence with the original spatiotemporal fusion feature vector of the node and input it into the fusion layer to generate an intermediate feature vector; introduce a residual connection mechanism to fuse the intermediate feature vector with the original node feature vector to complete the dynamic update of the node state and output the updated node feature vector containing multiple features. The standardized dynamic data at each time step (e.g., hla gene expression level of 0.82 at time t1) is multiplied by the corresponding temporal attention weight (α1=0.4) to obtain the attention-weighted node dynamic feature sequence (e.g., [0.12, 0.33, 0.21, 0.10]). This 1×4 dimensional sequence is then concatenated column-wise with the original 4500-dimensional spatiotemporal fusion feature vector of the node to form a 4504-dimensional vector, which is then input into the fusion layer (one fully connected layer, in...). features =4504,out features An intermediate feature vector is generated (=4500). A residual connection mechanism is introduced: the intermediate feature vector is element-wise added to the original node feature vector (torch.add(middle,original)) to mitigate information loss after feature fusion and achieve dynamic updating of node state. The updated node feature vector with 4500 dimensions is output, and the response strength of core features (such as virulence domain + high-weight dynamic expression) is improved by 28% compared with that before the update.

[0066] Step S146: Compare the updated node state with known pathogen and host interaction experimental data to evaluate the biological rationality of the node state update; for nodes that do not meet the preset requirements, reversely adjust the relevant parameters and contribution quantification dimensions of the temporal attention module; repeat the attention weight calculation and node state update process until the updated state of the core nodes can accurately reflect the dynamic change pattern and form a stable node state dynamic update mechanism.

[0067] The comparison metric is the Pearson correlation coefficient between the node feature response intensity and the Hla protein expression level detected by Western blot, with a preset threshold of ≥0.85. If the requirement is not met (e.g., a correlation coefficient of 0.78 for a certain virulence gene node), the parameters of the temporal attention module are adjusted in reverse: the learning rate of the attention layer is reduced from 0.0004 to 0.0003, and a contribution quantification dimension is added (a new expression level variation coefficient feature is added). The attention weight calculation and node update process is repeated. After 3 iterations, the correlation coefficients of all core nodes (Top 50) are ≥0.88, and those of non-core nodes are ≥0.80. Finally, a stable dynamic update mechanism for node states is formed, and the update results can accurately reflect the dynamic changes in virulence gene expression and immune response.

[0068] Step S150: Construct a full-chain correlation model based on dynamic fusion features and temporal network evolution results, analyze from three levels: features, process, and differences, quantify the differences in pathogen invasion mechanisms, and output an explanation report.

[0069] A full-chain correlation model based on Scikit-learn, dynamically fusing features and temporal evolution results, is built, with a random forest regressor at its core (n_estimators=200, max_depth=10). The model is implemented at three levels: at the feature level, the top 50 core differential features (such as hla and TLR2) are extracted, and their contribution percentages are quantified; at the process level, four invasion stages are broken down, and evolutionary synchronicity is analyzed; at the difference level, a multi-dimensional evaluation system is constructed. The model inputs are an updated node feature matrix (100 samples × 200 nodes × 4500 features) and temporal evolution data (dynamic edge weights, node state trajectories), and the output is a quantitative result of the pathogen invasion mechanism differences. Training uses 5-fold cross-validation, achieving an average accuracy of ≥92%. The generated explanation report includes a difference heatmap, an evolutionary trajectory visualization, and annotations of core conclusions (such as the difference in virulence between MRSA and MSSA mainly stemming from the high expression activation of the hla gene in phase t1), consistent with biological explanation logic.

[0070] Step S151: Collect three types of core data and align them according to sample identifiers and time steps, unify the format, divide feature subsets and filter input variables, and initialize model parameters and three-level analysis index system; Updated node dynamic features (S145), time-series network evolution data (S135, dynamic edge weights, evolutionary trajectories) and standardized dynamic data (S141) are aligned and linked by sample ID and time step, and stored uniformly in HDF5 format. Feature subsets are divided: virulence feature subset (80 pathogenic bacterial node features), immune feature subset (120 host node features), and interaction feature subset (features derived from dynamic edge weights). Input variables are selected (features with a variance contribution rate ≥ 0.1%, retaining a total of 3800). Model parameters are initialized: learning rate of 0.001 for the random forest regressor, weight decay of 1e-5; a three-level analysis index system is constructed: feature level (contribution ratio, mutual information value), process level (synchronicity coefficient, stage duration), and difference level (feature difference degree, process similarity, functional difference score), ensuring that the indexes cover the core analysis dimensions.

[0071] Step S152: Mine the core differential features of pathogenic bacteria, quantify the contribution ratio of features, analyze the correlation and functional classification of features, and compare the differences in feature expression and activation mode of different pathogenic bacteria. SHAP (SHapley Additive exPlanations) analysis was used to identify core differential characteristics of pathogenic bacteria and quantify the contribution percentage of each characteristic: the cumulative contribution percentage of the top 20 core characteristics was ≥75%, among which the core differential characteristics between MRSA and MSSA were hla (12%), psmα (8%), and agrA (6%). Pearson correlation analysis was used to explore the associations between characteristics: the correlation coefficient between hla and psmα expression was r=0.78 (p<0.001), classifying them as "cell lysis function"; the correlation coefficient between TLR2 and TNF-α was r=0.82 (p<0.001), classifying them as "immune recognition". Comparison of the characteristic patterns of different pathogenic bacteria: the hla gene expression level in MRSA at t1 phase (log2FC=2.3) was significantly higher than that in MSSA (log2FC=0.8), and the activation time of the agrA gene in MRSA was 2 hours earlier. This difference was verified by a t-test (p<0.01) and had clear biological significance.

[0072] Step S153: Break down the evolution data by invasion stage, quantify the synchronicity and differences in the invasion process, construct an evolutionary correlation map, and mine key turning points and stage difference patterns in the process. Evolutionary data were broken down into invasion stages (adhesion phase t0, invasion phase t1, proliferation phase t2, and diffusion phase t3) to quantify the synchronicity and heterogeneity of the invasion process. Synchronicity was assessed using the temporal correlation coefficient of bimodal interaction strength. MRSA and MSSA showed the highest synchronicity at t0 (r=0.85) and the greatest heterogeneity at t1 (r=0.52). Heterogeneity was assessed using the change in stage weights. The change in virulence-immunity interaction weights for MRSA from t1 to t2 (0.35) was 1.9 times that of MSSA (0.18). An evolutionary correlation map was constructed using Cytoscape, with node size representing feature contribution and edge thickness representing interaction strength. Key turning points were identified: the t1 to t2 stage was the core turning point for differences in invasion mechanisms (weighted mutation rate ≥40%). The stage difference pattern was characterized by rapid invasion and strong proliferation for MRSA, and slow invasion and weak proliferation for MSSA, consistent with the colony counting results from in vitro infection experiments.

[0073] Step S154: Construct a multi-dimensional difference assessment system, integrate to obtain a comprehensive difference score, verify the significance of differences and eliminate false differences, and construct a difference matrix of pathogen invasion mechanisms; A multi-dimensional difference assessment system was constructed, including 3 primary indicators and 8 secondary indicators: characteristic differences (differences in characteristic expression and functional annotation), process differences (differences in stage synchronicity and inflection point timing), and functional differences (differences in virulence intensity and immunosuppressive capacity). The Analytic Hierarchy Process (AHP) was used to determine the indicator weights, and a weighted summation was performed to obtain a comprehensive difference score (values ​​from 0 to 100). The comprehensive difference score between MRSA and MSSA was 68.5. The significance of the differences was verified by a t-test (p < 0.05), and the Benjamini-Hochberg method was used for FDR correction to eliminate 3 spurious differences (contribution < 0.5%). A pathogen invasion mechanism difference matrix was constructed (2 pathogens × 8 secondary indicators), with matrix elements being standardized difference values ​​(0-1), visually presenting the strength distribution of various differences and providing quantitative support for mechanism analysis.

[0074] Step S155: Construct causal relationship paths for features, processes, and differences, quantify the strength of associations at each stage, adjust parameters through cross-validation, and optimize the model's fit to biological interaction patterns. A causal pathway based on feature-process-difference was constructed using Bayesian networks: The core pathway is HLA high expression (feature) → accelerated invasion in T1 phase (process) → amplified virulence differences (difference). The association strength (conditional probability) of each step was quantified: the probability that HLA high expression leads to accelerated invasion is 0.83, and the probability that accelerated invasion leads to virulence differences is 0.79. Five-fold cross-validation was used to adjust model parameters: the structural complexity parameter α of the Bayesian network was set to 0.01 to ensure the pathway closely matches the KEGG pathogen-host interaction pathway model. The optimization objective was "the consistency between the association strength and experimental validation results ≥ 85%". After three rounds of parameter adjustments, the consistency of the core pathway reached 89%, and the consistency of secondary pathways (such as PSMα expression → immunosuppression differences) reached ≥ 82%, significantly improving the model's fit to the laws of biological interactions.

[0075] Step S156: Identify the biological functions and impacts of the core differences, transform the quantitative results into biologically meaningful statements, write an explanatory report with visualizations according to a standardized structure, and indicate the reliability of the conclusions; The biological functions and impacts of the core differences were analyzed: Based on GO / KEGG annotations, high expression of hla in MRSA corresponds to perforin activity (GO:0019838), leading to a 40% increase in host cell lysis rate; premature activation of agrA corresponds to the quorum sensing pathway (KEGG:02024), promoting the synergistic expression of virulence factors. The quantitative results were translated into a biologically meaningful statement: During the invasion phase (t1), MRSA enhances cell lysis capacity through high hla gene expression and premature activation of the agrA pathway, while simultaneously inhibiting host TNF-α release, forming a 'strong invasion-weak immune response' invasion pattern. The report was written according to a standardized structure: abstract, data sources and methods, three-level analysis results (including differential heatmaps and evolutionary trajectory diagrams), core conclusions, limitations, and the reliability of the conclusions (e.g., conclusions based on validation of 100 samples, p<0.01 after FDR correction, indicating high reliability).

[0076] Step S157: Evaluate the model accuracy using an independent validation set, supplement the validation of disputed differences, adjust the model algorithm and parameters, optimize the feature selection and process decomposition logic, and improve the accuracy and rationality of the interpretation.

[0077] The model accuracy was evaluated using 30 independent pathogenic bacteria (15 MRSA and 15 MSSA) as a validation set: the overall difference prediction accuracy was ≥90%, and the core feature difference identification accuracy was ≥93%. For controversial differences (such as functional differences in the psmα gene), Western blot and Transwell experiments were added for verification: the expression level of Psmα protein in MRSA was 2.1 times that in MSSA, and MRSA showed significantly stronger cell migration ability (p<0.01), validating the rationality of the controversial conclusions. The model algorithm and parameters were adjusted: the feature selection threshold was reduced from variance contribution rate ≥0.1% to ≥0.08%, and the duration of stage differences was added as a process decomposition indicator. After optimization, the model's interpretability improved by 15%. The final model achieved a biological significance agreement of ≥91%, and the logic and reliability of the difference interpretation met the needs of clinical and research applications.

[0078] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of this application.

Claims

1. A comparative analysis method for invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation, characterized in that: Includes the following steps: We used a combined normalization method to eliminate noise and scale differences in multi-source heterogeneous data from the genome, transcriptome, and proteome, extracted corresponding basic features using bioinformatics tools, and integrated them to construct a multimodal original feature library. A dual-branch attention mechanism is introduced to strengthen intra-modal features and establish bidirectional inter-modal correlations, and weighted fusion is used to generate a cross-modal feature set; The pathogen invasion process is abstracted into a temporal interaction network. A dynamic framework is built by embedding gated recurrent units into a temporal graph convolutional network, and the edge weights are dynamically adjusted to characterize the network evolution. The expression time-series curves of pathogenic bacteria virulence genes and dynamic data of host immune response are extracted and added to the temporal interaction network. The contribution of time nodes is highlighted through temporal attention, and the node status is dynamically updated. Based on dynamic fusion characteristics and temporal network evolution results, a full-chain correlation model is constructed. It analyzes from three levels: characteristics, process, and differences, quantifies the differences in pathogen invasion mechanisms, and outputs an explanatory report.

2. The method for comparative analysis of multi-pathogenic bacterial invasion mechanisms based on illumination-adaptive planar representation according to claim 1, characterized in that: The method employs a combined normalization approach to eliminate noise and scale differences in multi-source heterogeneous data from the genome, transcriptome, and proteome. Corresponding fundamental features are extracted using bioinformatics tools, and a multimodal raw feature library is constructed, including: For raw genome sequencing data, quality control measures such as adapter removal and low-quality base splicing are used. For transcriptome RNA-seq data, adapter sequence and ploy-N tail removal and short sequence filtering are performed by detecting base quality distribution, GC content and sequence repetition rate. For proteome mass spectrometry data, peptide identification and quantification are performed, and low-reliability protein entries are removed by database matching and false positive result filtering. For genomic data, sequence alignment correction was used to align the quality-controlled sequences to the corresponding pathogenic bacteria reference genome, remove unaligned, multiple-aligned, and PCR repetitive sequences, and correct base quality values. For transcriptome data, a combination of expression level normalization and batch effect correction was used to achieve data scale uniformity. For proteome data, a combination of protein structure homology filtering and quantitative standardization was used to remove homologous redundant proteins, correct quantitative differences between samples, and fill in missing values. Genome feature extraction involves homology comparison between the corrected sequence and the virulence gene database to extract virulence gene-related sequences and functional domain features. Transcriptome feature extraction uses differential expression analysis combined with functional enrichment analysis to screen transcripts and extract their related features. Proteome feature extraction involves screening differentially expressed proteins, predicting their functional domains, and extracting domain-related features and quantitative and functional annotation information. The basic features of each modality are standardized in format, and the features of different modalities are associated and mapped through sample identifiers to construct a sample and feature data table. A screening strategy of filtering low-variation features and removing highly correlated features is adopted to retain effective features.

3. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 1, characterized in that: The introduction of a dual-branch attention mechanism, which respectively strengthens intra-modal features and establishes bidirectional inter-modal correlations, and weightedly fuses them to generate a cross-modal feature set, includes: The dimensionality unification and normalization processes were performed on the genomic, transcriptomic, and proteomic features in the multimodal original feature library, respectively. The feature mapping algorithm was used to transform the heterogeneous features of each modality to the same high-dimensional feature space. The numerical scale difference of each modality feature was eliminated by the standardization method to obtain the standardized feature matrix corresponding to each modality. For each modality's standardized feature matrix, an intramodal self-attention calculation module is constructed. The query, key, and value matrices are generated through linear transformation. The intramodal attention weight matrix is ​​obtained based on similarity calculation. After normalization, it is weighted and summed with the value matrix to enhance functional features and suppress redundant features, outputting the intramodal enhanced feature matrix for each modality. A cross-modal cross-attention module is constructed. The enhancement feature matrices of any two modalities are taken as input, and mutual query matrix and mutual key matrix are generated by linear transformation respectively. The mutual information similarity between features of different modalities is calculated to generate a bidirectional cross-attention weight matrix. Based on the weight matrix, the enhancement features of another modality are weighted and filtered to obtain the inter-modal association feature matrix, realizing the bidirectional association mapping between genome and transcriptome, genome and proteome, and transcriptome and proteome. An adaptive weight adjustment mechanism based on gating units is introduced. According to the information contribution of each modality feature and the degree of correlation between modalities, the fusion ratio of intramodal self-attention weights and intermodal cross-attention weights is dynamically learned. The adjusted weight matrix is ​​then normalized twice to ensure the rationality of weight allocation and the balance of feature contributions. The intra-modal enhancement feature matrix and the corresponding inter-modal correlation feature matrix of each modality are summed element-wise according to the adaptively adjusted weights to obtain the preliminary fusion feature matrix of each modality. The preliminary fusion feature matrices of all modalities are concatenated to form a high-dimensional cross-modal candidate feature set. The risk of overfitting is reduced by regularization, and the cross-modal fusion feature set is finally output.

4. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 2, characterized in that: For each modality, a standardized feature matrix is ​​used to construct an intra-modal self-attention calculation module. This module generates query, key, and value matrices through linear transformation, calculates intra-modal attention weight matrices based on similarity, and after normalization, performs a weighted sum with the value matrix to enhance functional features and suppress redundant features. The resulting intra-modal enhanced feature matrix for each modality is output, including: For a single-modality normalized feature matrix, a three-layer linear transformation layer with shared parameters is constructed to generate query, key, and value matrices, respectively. The parameters of the transformation layer are initialized using a parameter initialization method to ensure convergence. A scaling factor based on the feature dimension is set to eliminate the scale bias in similarity calculation in high-dimensional feature space. The standardized feature matrices are input into the three linear transformation layers respectively, and a query matrix, key matrix and value matrix with uniform dimensions are generated by matrix multiplication. The dimensions of the three types of matrices are validated, and feature rows with abnormal dimensions are removed to ensure that the matrix dimensions match and the data is complete. A scaling dot product attention mechanism is used to perform matrix multiplication on the transpose of the query matrix and the key matrix to obtain the original feature similarity matrix; the original similarity matrix is ​​divided by the preset scaling factor to complete the scale normalization, generating a scaled similarity matrix, which quantifies the association strength of the corresponding feature pairs. The scaled similarity matrix is ​​input into the Softmax activation function, and the elements of each row of the matrix are normalized to generate an intramodal attention weight matrix, thereby realizing the differentiated weight allocation between functional features and redundant features. The attention weight matrix and the value matrix are multiplied to obtain the initial enhanced feature matrix. A residual connection mechanism is introduced to add the original standardized feature matrix to the initial enhanced feature matrix element by element after linear transformation. Then, the numerical offset is eliminated by layer normalization to output the intramodal enhanced feature matrix. Calculate the weight ratio and variance contribution rate of functional features in the enhanced feature matrix. If the feature weight increase does not reach the preset threshold, fine-tune the linear transformation layer parameters or scaling factor, and repeat the above steps to screen high contribution dimensions in the enhanced feature matrix and remove low contribution redundant dimensions, and finally determine the intramodal enhanced feature matrix of each modality.

5. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 2, characterized in that: The construction of the cross-modal cross-attention module takes the enhancement feature matrices of any two modalities as input, generates a mutual query matrix and a mutual key matrix through linear transformation, calculates the mutual information similarity between features of different modalities to generate a bidirectional cross-attention weight matrix, and performs weighted filtering on the enhancement features of another modality based on this weight matrix to obtain the inter-modal association feature matrix, realizing bidirectional association mapping between genome and transcriptome, genome and proteome, and transcriptome and proteome, including: Two independent linear transformation layers corresponding to the bimodality are constructed, and the transformation layer parameters are initialized to ensure mapping stability. A mutual information calculation window is set to provide support for similarity measurement. A dual-modal enhanced feature matrix is ​​selected as the input pair, and four sets of mutual query and mutual key matrices are generated through linear transformation. The dimension consistency check ensures that the query and the corresponding key matrix column dimension match. Based on the generated mutual query and mutual key matrix, a two-way mutual information calculation process is constructed. The original association matrix and mutual information value are integrated to generate a two-way mutual information similarity matrix, which quantifies the association strength between the two modalities. The bidirectional mutual information similarity matrix is ​​input into the Softmax activation function for normalization, generating a bidirectional cross-attention weight matrix that reflects the contribution of the source mode to the target mode. By performing matrix multiplication operations between the bidirectional attention weight matrix and the other modality enhancement feature matrix, strong correlation features are filtered to obtain two sets of correlation feature matrices, which together constitute the intermodal correlation feature matrix of the bimodal pair. The genome and transcriptome, genome and proteome, and transcriptome and proteome were selected sequentially as bimodal input pairs. The aforementioned steps were repeated to generate bidirectional association feature matrices for each group, thus completing the bidirectional association mapping of the full cross-modal combination.

6. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 1, characterized in that: The process of pathogen invasion is abstracted into a temporal interaction network. A dynamic framework is built by embedding gated recurrent units into a temporal graph convolutional network, and the edge weights are dynamically adjusted to characterize the network evolution, including: Divide the time steps according to the pathogen invasion stage, define the types of network nodes and edges, collect interaction data of each time step to construct a static network snapshot, arrange them in time sequence to form an initial temporal interaction network and label the relevant attributes; Nodes are initialized with pathogenic and host features fused across modalities, and initial weights are assigned to edges based on experimental data and literature information. A time-step correlation matrix is ​​constructed to quantify the correlation degree of node states. A graph convolution module is built to extract spatial features of network snapshots, and a gated recurrent unit is embedded to integrate temporal dependency features, forming a spatiotemporal fusion structure and designing a forward propagation path. An edge weight update function is constructed based on the temporal characteristics of nodes. A differentiated adjustment threshold is set in combination with the biological interaction rules. The edge weight is adaptively updated as the invasion process is achieved by adjusting the coefficient. Input network snapshots into the fusion framework at time steps, output dynamic edge weights and node states, record evolutionary trajectories to generate evolutionary graphs, and highlight key evolutionary nodes and stages through dimensionality reduction visualization; Experimental data was collected and divided into training and validation sets. Framework parameters were optimized using a specific loss function, and the effectiveness of the core modules was verified through ablation experiments.

7. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 6, characterized in that: The constructed graph convolutional module extracts spatial features from the network snapshot, embeds gated recurrent units to integrate temporal dependency features, forms a spatiotemporal fusion structure, and designs the forward propagation path, including: First-order approximate graph convolution is selected as the operator. The weight matrix of the graph convolution layer is defined based on the node feature dimension of the temporal interaction network snapshot and the preset spatial feature output dimension. The weight matrix is ​​initialized using a parameter initialization method. A self-loop correction adjacency matrix processing function is constructed to add an identity matrix to the original network adjacency matrix and perform normalization processing to eliminate feature extraction bias caused by node degree differences. The activation function and dropout rate of the graph convolution module are set. The network snapshot at a single time step is decomposed into a node feature matrix and a normalized adjacency matrix. The node feature matrix is ​​input into the graph convolutional layer to perform matrix operations to extract spatial features and output the spatial feature matrix at that time step. The above operation is performed on the network snapshots at all time steps in sequence to obtain the spatial feature sequence, and the dimension of the feature sequence is verified. The GRU hidden layer dimension is set based on the spatial feature output dimension. The reset gate, update gate, weight matrix and bias term corresponding to the candidate hidden state of GRU are initialized. The GRU input adaptation layer is designed to convert the spatial feature matrix of each time step into a format that matches the input requirements of the GRU time series sequence, thus forming the GRU input sequence. The output of the graph convolution module is connected to the input of the GRU module to construct a serial fusion structure of graph convolution and GRU. A residual connection mechanism is introduced to linearly transform the spatial features output by the graph convolution and add them element-wise to the hidden state of the current time step output by the GRU, and then perform layer normalization. A fully connected layer is added after the fusion structure to map the dimension of the hidden state output by the GRU back to the dimension that matches the node features, forming a complete spatiotemporal fusion feature extraction module. The complete data flow path of forward propagation is clearly defined, which is as follows: the input layer receives network snapshots at each time step, the graph convolutional layer extracts spatial features, the adaptation layer converts the input format, the GRU layer calculates and outputs the hidden state in the order of time steps, and the fusion output layer outputs the spatiotemporal fusion features; the computational priority and data dependency of each layer are defined. A small number of temporal network snapshots are selected and input into the fusion module to verify the feature dimension matching step by step. If dimension incompatibility issues occur, the output dimension of the graph convolution or the dimension of the GRU hidden layer is adjusted and the weight parameters are reinitialized.

8. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 6, characterized in that: The method of constructing an edge weight update function based on the temporal characteristics of nodes, setting a differentiated adjustment threshold in conjunction with biological interaction patterns, and achieving adaptive updating of edge weights as the invasion process is achieved through adjusting coefficients, including: The spatiotemporal fusion feature vectors of the same node at adjacent time steps are extracted, and standardized node temporal feature difference values ​​are obtained through feature difference degree calculation and normalization. An edge weight update function is constructed with node feature difference degree and GRU temporal association weight as dual inputs. It is designed in two layers: basic adjustment coefficient calculation and biological rule correction, to ensure that the coefficients conform to the laws of biological interaction dynamics. The network edges are classified according to the type of biological event. Based on experimental data and pathway models, differentiated state change trigger thresholds and upper and lower limits for weight adjustment are set for different types of edges. Traverse the network edges at the current time step, calculate the basic adjustment coefficient and determine whether it exceeds the trigger threshold. Obtain the final adjustment coefficient through biological rule correction and truncation. The initial edge weights and the final adjustment coefficients are used to calculate the dynamic edge weights. After outlier removal and smoothing filtering, the dynamic edge weights are synchronously updated to the network snapshot at the corresponding time step. The accuracy of the weight update is evaluated by comparing with actual measured data, and the relevant parameters are adjusted in reverse. Through multiple rounds of iterative optimization, a stable adaptive update mechanism for edge weights is formed.

9. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 1, characterized in that: The extraction of pathogenic bacterial virulence gene expression time-series curves and host immune response dynamic data is supplemented into a temporal interaction network. Temporal attention is used to highlight the contribution of time nodes, enabling dynamic updates of node states, including: Collect time-series data on the expression of pathogenic bacteria virulence genes and dynamic data on the host immune response; perform noise reduction and normalization on the raw data, align the two types of data according to the preset time step of pathogen invasion, fill in the data of missing time steps, and generate a standardized time-series data sequence that corresponds one-to-one with the time step of the time-series interaction network. A node-dynamic data association dictionary is constructed to map the time-series data of virulence gene expression to the pathogenic virulence factor nodes in the time-series interaction network, and to map the dynamic data of host immune response to the corresponding host-related nodes; a unique time-series dynamic data sequence is bound to each node to generate a three-dimensional data table; the uniqueness and completeness of the mapping relationship are verified, and redundant dynamic data is removed. Design a temporal attention layer based on the contribution of time nodes, define the dimensions for calculating attention scores and output dimensions, initialize the parameters of the weight matrix and bias terms of the attention layer, define the logical framework for attention calculation, clarify the input and output content, and set the dropout rate of the attention module. Relevant features of dynamic data at each time step node are extracted to quantify the contribution of a single time step node to the overall state of the node. Attention scores are calculated using an attention mechanism, and the original attention scores are obtained by processing the contribution features of time steps with the spatiotemporal fusion features of the nodes. The original scores of all time steps are normalized to generate temporal attention weights for each time step, highlighting the role of high-contribution time steps. The temporal attention weights at each time step are weighted with the standardized dynamic data of the corresponding nodes to obtain the attention-weighted node dynamic feature sequence. This sequence is then concatenated with the original spatiotemporal fusion feature vector of the node and input into the fusion layer to generate an intermediate feature vector. A residual connection mechanism is introduced to fuse the intermediate feature vector with the original node feature vector to complete the dynamic update of the node state and output the updated node feature vector containing multiple features. The updated node states are compared with known pathogen and host interaction experimental data to assess the biological rationale for the node state update. For nodes that do not meet the preset requirements, the relevant parameters and contribution quantification dimensions of the temporal attention module are adjusted in reverse. The attention weight calculation and node state update process is repeated until the updated node state can accurately reflect the dynamic change pattern, forming a stable dynamic node state update mechanism.

10. The method for comparative analysis of invasion mechanisms of multiple pathogens based on illumination-adaptive planar representation according to claim 1, characterized in that: The model, constructed based on dynamic fusion features and temporal network evolution results, analyzes the entire chain of associations from three levels: features, processes, and differences. It quantifies the differences in pathogen invasion mechanisms and outputs an explanatory report, including: Collect three types of data and align them according to sample identifiers and time steps, unify the format, divide feature subsets and filter input variables, and initialize model parameters and a three-level analysis index system; To uncover the differential characteristics of pathogenic bacteria, quantify the contribution ratio of characteristics, analyze the correlation and functional classification of characteristics, and compare the differences in characteristic expression and activation patterns of different pathogenic bacteria. The evolutionary data is broken down by invasion stage, the synchronicity and differences in the invasion process are quantified, an evolutionary correlation map is constructed, and key turning points and stage difference patterns in the process are explored. A multi-dimensional difference assessment system was constructed, and a comprehensive difference score was obtained by integrating the results. The significance of the differences was verified and false differences were eliminated. A difference matrix of pathogen invasion mechanisms was constructed. Construct causal relationship paths for features, processes, and differences, quantify the strength of associations at each stage, adjust parameters through cross-validation, and optimize the model's fit to the laws of biological interactions. The biological functions and effects of the differences were analyzed, the quantitative results were transformed into biologically meaningful expressions, and an explanatory report with visualizations was written according to a standardized structure, with the reliability of the conclusions indicated. The model's accuracy was evaluated using an independent validation set, and points of contention were supplemented for verification. The model's algorithm and parameters were adjusted, and the logic for feature selection and process decomposition was optimized to improve the accuracy and rationality of the interpretation.

Images