A method for inferring specific regulatory modulation of glioblastoma cells

Abstract

This invention discloses a method for inferring cell-specific regulation in glioblastoma. It includes: a multi-source data acquisition and preprocessing module, a dual-channel gene feature embedding module, a multi-view graph attention feature update module, a channel feature adaptive fusion module, and a gene regulatory relationship prediction output module. This invention employs a dual-channel architecture combining Transformer channels and autoencoder channels, enabling effective modeling of long-range dependencies, nonlinear interaction patterns, and contextual associations across cell states from single-cell expression data. This overcomes the limitations of traditional methods in characterizing complex regulatory logic and is more suitable for analyzing the highly dynamic and multi-level transcriptional regulatory processes in glioblastoma. The model of this invention exhibits superior robustness and generalization ability, making it more valuable for biological research and clinical applications in glioblastoma.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of bioinformatics, computational biology and artificial intelligence. Specifically, it relates to a method for inferring glioblastoma cell type-specific gene regulatory networks based on single-cell transcriptome data using deep learning technology, particularly by combining convolutional neural networks, autoencoders and graph attention mechanisms. Background Technology

[0002] Glioblastoma (GBM) is one of the most common and most malignant primary brain tumors of the central nervous system, characterized by high invasiveness, diffuse infiltrative growth, rapid progression, and poor prognosis. Its pathogenesis is extremely complex, involving abnormal regulation at multiple levels, including the genome, transcriptome, epigenome, and tumor microenvironment. Gene Regulatory Networks (GRNs), as topological models describing the interactions between intracellular transcription factors (TFs) and target genes (TGs), are core tools for elucidating the molecular mechanisms of GBM development, understanding tumor cell heterogeneity, revealing tumor invasion and treatment resistance characteristics, and discovering potential therapeutic targets.

[0003] Traditional gene regulatory network inference relies primarily on bulk RNA-seq data. However, bulk data only reflects average gene expression signals at the tissue level, masking the heterogeneity among different cell types (such as tumor cells, immune cells, and stromal cells) in the tumor microenvironment and failing to reveal cell-specific regulatory landscapes. With the rapid development of single-cell RNA sequencing (scRNA-seq) technology, researchers can acquire massive amounts of gene expression data at single-cell resolution, providing an unprecedented data foundation for constructing cell-type-specific GRNs.

[0004] Nevertheless, directly inferring GRNs from scRNA-seq data still faces significant technical challenges: data sparsity and noise interference. scRNA-seq data generally suffers from severe "zero-inflation" or "dropout" phenomena, and has a low signal-to-noise ratio. Traditional inference methods based on correlation or mutual information are highly susceptible to this noise, leading to high false positive rates and difficulty in capturing true regulatory signals.

[0005] Insufficient utilization of prior knowledge: Many existing inference methods (such as GENIE3, GRNBoost2, etc.) are mainly data-driven, ignoring the large amount of prior regulatory knowledge accumulated in the field of biology (such as TF binding information in the ChIP-seq database). Although some methods (such as SCENIC+) attempt to integrate multi-omics data (such as scATAC-seq), the high cost and technical difficulty of acquiring paired multi-omics data limit their widespread application.

[0006] Lack of heterogeneity modeling capability: Most existing deep learning models are designed for general benchmark datasets and lack targeted modeling of the specific tumor microenvironment and cell lineage specificity of glioblastoma. They often fail to effectively capture indirect regulation (co-regulation) or higher-order regulatory patterns, resulting in insufficient inference accuracy in complex disease contexts.

[0007] Therefore, there is an urgent need to develop a high-precision gene regulatory network inference method that can effectively integrate single-cell expression data with multi-source prior knowledge, extract complementary features through deep learning, and adapt to the high heterogeneity of glioblastoma. (Invention Content)

[0008] The purpose of this invention is to provide a method for inferring cell-specific gene regulatory networks in glioblastoma cells by integrating variational autoencoders and graph Transformers. This method integrates single-cell transcriptome data with information from public transcription factor databases, designs a unique dual-channel feature extraction architecture and a multi-view graph attention update mechanism, and achieves accurate identification of cell-type-specific regulatory relationships in glioblastoma.

[0009] This invention includes the following modules: a multi-source data acquisition and preprocessing module, a dual-channel gene feature embedding module, a multi-view graph attention feature update module, a channel feature adaptive fusion module, and a gene regulatory relationship prediction output module.

[0010] The multi-source data acquisition and preprocessing module is used to acquire single-cell sequencing data and prior transcription factor regulatory information of glioblastoma, filter low-quality cells through multi-level quality control, identify tumor cells, immune cells, glial cells and other cell subpopulations by combining dimensionality reduction clustering technology, and screen for hypervariable genes specific to each cell type to construct a prior regulatory network framework specific to glioblastoma cell type.

[0011] The dual-channel gene feature embedding module employs a parallel Transformer channel and an autoencoder (AE) channel. The Transformer channel utilizes a self-attention mechanism to model long-range dependencies and complex regulatory patterns along the gene-cell relationship axis, effectively capturing the contextual features of different genes in multicellular states. The AE channel learns the global co-expression structure by reconstructing the cell-gene expression profile and converts the decoder weights into a global embedding representation of the gene.

[0012] The multi-view graph attention feature update module constructs three complementary graph structures based on prior regulatory information, gene expression association patterns, and functional similarity relationships: a prior regulatory graph, regulatory association graph, and functional similarity graph. A multi-head graph attention mechanism is used to propagate and weighted aggregate node features on these three graph structures, capturing explicit regulatory relationships, potential association patterns, and functional module information from different perspectives. The feature representations from each perspective are then fused to obtain more comprehensive and biologically interpretable node embedding results.

[0013] The aforementioned channel feature adaptive fusion module introduces a channel attention mechanism, dynamically calculates dual-channel weights based on the importance of different feature channels to the task, and adaptively fuses global structural features and local contextual features to generate the final node embedding representation.

[0014] The gene regulation relationship prediction output module, based on the fused node embedding representation, predicts the probability of connections between TF human transcription factors and target genes through a deep neural network, and constructs the final gene regulation network.

[0015] Compared with the prior art, the present invention has the following advantages:

[0016] Enhanced ability to characterize complex regulatory relationships: This invention employs a dual-channel architecture combining Transformer and autoencoder channels, which can effectively model long-range dependencies, nonlinear interaction patterns, and contextual associations across cell states from single-cell expression data. This breaks through the limitations of traditional methods in characterizing complex regulatory logic and is more suitable for analyzing highly dynamic and multi-level transcriptional regulation processes in glioblastoma.

[0017] Superior Model Robustness and Generalization Ability: This invention effectively reduces the interference of high noise, high dropout rate and batch differences in single-cell sequencing data on inference results through the collaborative design of multi-source prior knowledge constraints, dual-channel complementary representation learning and multi-adjacency graph attention propagation mechanism, thereby improving the stability and generalization performance of the model in different datasets, different cell subpopulations and complex tumor sample backgrounds.

[0018] The invention has significant value for biological research and clinical applications in glioblastoma: It can construct cell type-specific regulatory networks for different cell populations in glioblastoma, such as tumor cells, immune cells, astrocytes, and oligodendrocytes. This helps to identify key regulatory factors related to tumor invasion, stemness maintenance, treatment tolerance, and immune escape, providing strong support for molecular subtyping of glioblastoma, screening of potential drug targets, and formulation of personalized treatment strategies. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall framework of the method of the present invention. Detailed Implementation

[0020] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. This embodiment uses a publicly available glioblastoma single-cell dataset (such as GSE131928) as an example to illustrate the specific implementation steps of the present invention. Figure 1 As shown, the method specifically includes:

[0021] I. Multi-source data acquisition and preprocessing module. This module is the cornerstone of the entire method and aims to provide high-quality, low-noise input data.

[0022] (1) Data collection: Collect single-cell RNA sequencing (scRNA-seq) data of tumor tissue and adjacent normal tissue of glioblastoma patients; download the list of human transcription factors (TF) and their target gene (TG) regulatory relationship data from authoritative databases such as hTFtarget.

[0023] (2) Single-cell data quality control (QC):

[0024] To eliminate interference from low-quality cells (such as dead cells, double cells, or empty droplets) in subsequent analyses, this invention establishes a stringent four-stage filtration standard:

[0025] Gene detection number filtering: Count the number of genes detected in each cell (nFeature_RNA), and retain cells with a number of 200 to 2000. Cells with fewer than 200 may be empty droplets, and cells with more than 2000 may be double cells.

[0026] Sequencing depth filtering: Count the total number of UMIs (nCount_RNA) for each cell, retain cells with less than 5000, and further exclude duplicate cells.

[0027] Mitochondrial gene proportion filtering: Calculate the proportion of mitochondrial gene expression to total expression (percent.mt). Cells with a proportion greater than or equal to 15% are removed. A high mitochondrial proportion usually indicates that the cells are in a state of apoptosis or rupture.

[0028] Ribosome gene ratio filtering: Retain cells with a ribosome gene ratio (percent.ribo) of 10% to 30% to ensure that the cells are in a normal metabolically active state.

[0029] (3) Data standardization and screening for hypervariable genes:

[0030] Data standardization: The filtered expression matrix was normalized using the LogNormalize method. The normalized and logarithmically transformed gene expression values ​​were then used. , x represents the original expression count of a gene in a cell, Total represents the total expression count of all genes in the cell, i.e. the total UMI or sequencing depth of the cell; 10000 represents the scaling factor, used to standardize the total expression level of different cells to the same order of magnitude; through the above processing, the systematic bias caused by the difference in sequencing depth between different cells can be effectively reduced.

[0031] Cell annotation: Dimensionality reduction (PCA+UMAP) and clustering (Louvain algorithm) were performed using the Seurat workflow. Cells were annotated by combining typical marker genes (such as SOX2, OLIG2, EGFR, AIF1, CD68, C1QA, CD3D, NKG7, MBP, PLP1, GFAP, AQP4, CLDN5, PECAM1, etc.) and the cells were divided into six categories: malignant tumor cells, myeloid cells, T / NK cells, oligodendrocytes, astrocytes and endothelial cells.

[0032] Specific feature selection: For each annotated cell type, the variance of gene expression was calculated independently, and the top 500 genes with the highest degree of variation (Hvgs) were selected. These 500 genes represent the most significant biological changes in that cell type and are defined as the candidate target gene set.

[0033] (4) Construct a cell-specific prior network and use the "intersection mapping" strategy to construct the prior network framework:

[0034] The intersection of TF-TG regulatory pairs from the public database with the 500 hypervariable genes selected above was taken. That is, only the regulatory edges of those target genes belonging to hypervariable genes of this cell type were retained. The resulting network contains both general regulatory knowledge and cell type specificity. Based on this, three adjacency matrices were further constructed for subsequent module input.

[0035] II. Dual-channel Gene Feature Embedding Module. Addressing the inherent high dimensionality, sparsity, and high noise characteristics of single-cell transcriptome data (scRNA-seq), this invention designs a parallel dual-channel feature extraction framework. It utilizes a Transformer encoder to capture long-range dependencies in the sequence, and a variational autoencoder (VAE) to learn a robust probabilistic global co-expression latent representation, thereby generating an informative initial gene embedding vector.

[0036] (1) Construct the Transformer encoder channel:

[0037] This channel aims to treat the expression value of a single gene across all cells as a sequence, using a self-attention mechanism to capture complex contextual dependencies and long-distance co-expression patterns across different cell populations.

[0038] Input: a single first One gene in all Expression profile vector in each cell e j,n This represents the expression value of the j-th gene in the n-th cell; the superscript T indicates the transpose operation, which converts the row vector into a column vector.

[0039] Positional encoding and embedding: First, the input scalar representation values ​​are mapped to a high-dimensional vector sequence, and then positional encoding is superimposed to preserve the sequence information of cell arrangement (if the data has been sorted in a specific way) or enhance the model's ability to distinguish different input positions.

[0040] Multi-Head Self-Attention: The core component. By calculating attention weights between different positions in the sequence, the model can focus on the subset of cells that are most informative of the current gene expression pattern, capturing global contextual information. The multi-head mechanism allows the model to learn different features in parallel across different representation subspaces.

[0041] Feed-Forward Network: After each attention operation, the feature vector at each location is independently transformed non-linearly, enhancing the expressive power of the model.

[0042] Residual connections and layer normalization (Add & Norm): Residual connections and layer normalization are introduced at the output of each sub-layer (attention layer and feedforward layer) to accelerate training convergence and alleviate the gradient vanishing problem, thus building a deep network structure.

[0043] Output: After After processing by the Transformer encoder block, the output feature sequence is aggregated (e.g., by global average pooling) to obtain a fixed-length feature vector for the gene. .

[0044] Calculation process: For the th Layer, input is The output is The calculation is as follows:

[0045]

[0046]

[0047] Final gene embedding

[0048] AttentionOutput represents the intermediate feature representation obtained after multi-head self-attention computation, residual connections, and layer normalization. It reflects the updated result of the input features after global dependency modeling and can be regarded as the output of the first sub-layer of the Transformer module. LayerNorm represents the layer normalization operation, which normalizes the feature dimension of each sample to make the feature distribution more stable, thereby accelerating model convergence and improving training stability. MultiHeadAttention represents the multi-head attention mechanism, which learns the dependencies between input features in parallel from multiple different representation subspaces, thereby capturing richer global contextual information. FeedForward represents the feedforward neural network.

[0049] (2) Constructing the Variational Autoencoder (VAE) channel:

[0050] This channel aims to learn the global co-expression potential structure of genes from noisy cellular expression data using generative modeling methods. Unlike traditional autoencoders, VAEs map the input to a probability distribution in the latent space, thereby learning smoother and more robust feature representations.

[0051] Network architecture: Construct a probabilistic encoder and a generative decoder based on cell expression vectors.

[0052] Probabilistic encoder: Receives high-dimensional cell expression vectors (Contains all genes in the cell) The expression value in the model is used as input. Through a multi-layer neural network, the parameters of the posterior distribution of the cell in the low-dimensional latent space are estimated, typically the mean vector of a Gaussian distribution. Sum of logarithmic variance vector , i represents the i-th cell, This represents the mean vector of the distribution corresponding to the i-th cell in the low-dimensional latent space. Let represent the variance vector of the i-th cell in the potential spatial distribution.

[0053] Reparameterization Trick: To enable gradient backpropagation through a random sampling process from the standard normal distribution. Sample a noise vector and then transform it. Generate latent variables , z i This represents the low-dimensional representation of the i-th cell, used for subsequent decoding or inference.

[0054] Generative decoder: samples latent variables Map back to the original high-dimensional gene space to reconstruct the gene expression vector of the cell. .

[0055] Training objective: Maximize the lower bound of evidence (ELBO). Its loss function consists of two parts: one is the reconstruction error (such as negative binomial likelihood loss or MSE), which is used to ensure the quality of generation; the other is the KL divergence regularization term, which constrains the learned latent posterior distribution to be close to the prior standard normal distribution, thus preventing the model from overfitting.

[0056] Gene feature extraction (key step): Similar to the design of autoencoder channels, although VAEs are trained at the cellular level, their weight matrices encode the relationships between genes. This invention extracts the weight matrix of the first fully connected layer of the encoder. . For the real number field, N genes d represents the number of genes. hidden This matrix represents the number of neurons in the first hidden layer of the encoder. Each row of this matrix corresponds to the contribution weight of a gene to the first hidden unit. Therefore, this weight matrix directly constitutes the gene feature embedding matrix based on the global probability model. .

[0057] Brief description of the calculation formula:

[0058] Encoding process: EncoderNN is an encoder neural network, a non-linear mapping function responsible for mapping high-dimensional input genes to a low-dimensional latent space.

[0059] Sampling process: Follow the mean Covariance is The multivariate distribution By reparameterization accomplish; This indicates element-wise multiplication. This represents the normally distributed noise from random sampling.

[0060] Decoding process: DecoderNN is a decoder neural network that takes vectors inferred from the latent space and reconstructs them into the original representation data.

[0061] Gene embedding acquisition: ,in The first layer weights of the encoder.

[0062] III. Multi-view graph attention feature update module:

[0063] To more effectively utilize prior regulatory information and data-driven relationship patterns to guide node feature propagation, this invention constructs three complementary graph structures and employs a multi-head graph attention mechanism to iteratively update node representations, thereby enhancing the ability to characterize the complex regulatory relationships in glioblastoma.

[0064] (1) Construct three types of complementary graph structures:

[0065] ① Prior Regulation Graph: Based on a public transcription factor regulation database and cell type-specific expression screening results, a prior regulation graph was constructed, with transcription factors and target genes as nodes and known regulatory relationships as edges. This graph is used to characterize high-confidence explicit regulatory relationships, providing biological constraints for subsequent feature propagation.

[0066] ② Regulatory Relationship Graph: Based on the correlation, covariation, and local topological neighborhood relationships between genes in single-cell expression data, a regulatory relationship graph reflecting potential synergistic regulation, cascade regulation, and implicit association patterns is constructed. This graph is not limited to the direct connections already recorded in the database and can supplement potential indirect regulatory information.

[0067] ③ Functional similarity graph: Based on the similarity of transcription factors or target genes in regulatory patterns, functional annotation, or expression behavior, a functional similarity graph is constructed to connect gene nodes with similar biological functions or in similar regulatory programs. This graph can enhance the model's ability to identify modular functional structures and common regulatory rules.

[0068] (2) Feature update based on Graph Transformer with three types of graph structures:

[0069] To more effectively fuse multi-source prior information and capture complex nonlinear feature interactions, this invention abandons the traditional Graph Attention Network (GAT) and instead adopts a more powerful Graph Transformer architecture. This architecture explicitly integrates the prior graph structure into the Transformer's self-attention mechanism and achieves deep feature propagation through multi-layer stacking.

[0070] This module receives the gene feature matrix from the previous stage. As input (initial input is) and (the splicing), and through The Graph Transformer layer handles the processing. The computation flow for each layer is as follows:

[0071] ① Multi-Head Graph Self-Attention with Incorporating Prior Structure

[0072] The standard Transformer self-attention mechanism computes attention across all node pairs, forming a fully connected graph. To leverage prior knowledge, this invention modulates the calculation of attention scores through graph structural bias.

[0073] For the h-th attention head, calculate the attention score between node i and node j:

[0074]

[0075] in, This represents the influence weight of gene i on gene j under the h-th attention head. and This represents the feature vectors of gene i and gene j in layer l-1. and These are learnable query and key projection matrices, which are used to project genetic features into different subspaces. It is the dimension of the key vector. The key innovation lies in... The term is a structural mask used to restrict which genes can generate "attention," based on nodes i and j in four prior networks. Define the connection relationships in the table:

[0076]

[0077] This means that as long as at least one of the three networks G1, G2, and G3 contains the edge (i,j), Let G represent the k-th prior graph.k All prior edges in the set, i=j represents a self-loop connection of a gene.

[0078] This design is equivalent to introducing a hard attention mask, which forces the model to focus only on node pairs (or itself) that are connected in at least one prior network, thereby hardcoding the prior topology into the attention mechanism.

[0079] Subsequently, normalized attention weights are calculated using the Softmax function, and the value vector is aggregated:

[0080]

[0081]

[0082] This represents the regulatory weight of gene i on gene j under the h-th attention head. represents the original matching score between gene i and gene j, and Softmax is the normalization function. This indicates that the h-th attention head calculates the new feature vector of gene i. Gene j transmits the specific biological information of gene i.

[0083] Finally, the outputs of all K heads are concatenated and mapped through a linear layer to obtain the output of multi-head attention.

[0084] ② Residual Connectivity and Layer Normalization (Add & Norm)

[0085] To alleviate the vanishing gradient problem and accelerate convergence during deep network training, residual connections and layer normalization are introduced at the output of the attention sublayer:

[0086] , This represents the intermediate output features of this sublayer. This represents the gene feature matrix output by the previous layer of the network. This represents the attention gene in a multi-head graph, responsible for capturing the complex relationships between genes.

[0087] ③ Feed-Forward Network

[0088] Subsequently, the features are subjected to a pointwise nonlinear transformation through a feedforward network containing two fully connected layers and a ReLU activation function to enhance the model's feature extraction capability:

[0089] FFN represents a feedforward neural network, ReLU represents a non-linear activation function, and W1, b1 and W2, b2 are the weight matrices and bias vectors of the first and second fully connected layers, respectively.

[0090] Similarly, this sub-layer is followed by a residual connection and layer normalization operation to obtain the final output of the current layer:

[0091] , This represents the final gene feature representation after processing by the l-th transformer layer. This indicates that feature extraction is further performed on the output of the attention layer.

[0092] After iterative updates of L layers of Graph Transformer blocks, the output of the last layer... This refers to the high-level gene feature representation learned by this module, which incorporates multiple prior topological information.

[0093] IV. Channel Feature Adaptive Fusion Module:

[0094] This module aims to adaptively fuse two types of high-level gene feature representations output by the dual-channel gene feature embedding module and the multi-view graph attention feature update module. One type of feature is a feature representation obtained based on the variational autoencoder channel and updated using a graph structure. Another type of feature is the feature representation obtained based on Transformer channels and updated by graph structure. Since the two focus on characterizing the global co-expression probability distribution and complex context-dependent features in gene expression data, respectively, this module introduces an attention mechanism to dynamically model the feature contributions of the two channels and generate a unified fusion feature representation in order to fully leverage their complementary advantages. This serves as the input for the subsequent gene regulatory relationship prediction output module.

[0095] Let the updated features of the k-th gene node in the two channels be respectively and First, calculate the attention scores for each of the two channels:

[0096]

[0097]

[0098] in, and For learnable weight matrices, and Here, 'a' is the bias term, and 'a' is the shared attention vector. Represents a nonlinear activation function. and represents the importance score of the k-th gene in the variational autoencoder channel and the Transformer channel, respectively.

[0099] Subsequently, the scores of the two channels are normalized using the Softmax function to obtain the corresponding adaptive fusion weights:

[0100]

[0101]

[0102] in, and Let represent the normalized weights of the k-th gene in the two channels, and satisfy . .

[0103] The final fusion feature representation of the k-th gene is defined as follows:

[0104]

[0105] After fusing all gene nodes, the final node feature matrix can be obtained:

[0106]

[0107] Where N represents the number of gene nodes, This is the unified feature representation matrix after fusion.

[0108] V. Gene Regulation Relationship Prediction Output Module:

[0109] This module aims to utilize the high-level gene feature representation generated by the aforementioned channel feature adaptive fusion module to uniformly model the multi-source information learned by the dual-channel gene feature embedding module and the multi-adjacency graph Transformer feature update module. The high-level gene feature representation integrates sequence context features, global co-expression structure probability distribution features, and high-order graph structure features after prior topological propagation, and is used to quantitatively predict whether a regulatory relationship exists between any pair of genes (specific transcription factor i and target gene j).

[0110] (1) Scoring mechanism based on bilinear decoder:

[0111] The input is the final gene feature matrix output by the channel feature adaptive fusion module in the previous stage. For the potential regulatory relationships to be predicted (i.e., directed edges from transcription factor node i to target gene node j), their corresponding feature vectors are extracted. .

[0112] To effectively capture the complex and asymmetric interactions between nodes, this invention abandons simple inner product or general multilayer perceptron (MLP) methods, and instead employs a bilinear decoder specifically designed for link prediction tasks. This decoder introduces a learnable interaction weight matrix, which can explicitly model the directed interactions between features of two nodes.

[0113] The formula for calculating the raw score of transcription factor i regulating target gene j is as follows;

[0114]

[0115] in, It is a learnable bilinear weight matrix used to capture the interaction relationships between feature dimensions; This represents the transpose of the transcription factor eigenvector.

[0116] (2) Probability Mapping and Final Prediction

[0117] To transform the raw scores into probability values ​​representing the possibility of a regulatory relationship, this invention maps the scores to a sigmoid activation function. Interval:

[0118] , This represents the predicted probability that gene i has a regulatory relationship with gene j. It is the raw score of the correlation between gene i and gene j predicted by the model.

[0119] Final output This represents the probability predicted by the model that transcription factor i regulates target gene j. In practical applications, a threshold (such as 0.9) can be set to classify gene pairs with a regulatory relationship as having a relationship, thereby constructing the final gene regulatory network.

[0120] (3) Model training objective (loss function):

[0121] Given the high sparsity of real gene regulatory networks—meaning the number of real regulatory edges is far less than the number of non-existent edges—the proportion of positive and negative samples in the training data becomes extremely unbalanced. If a traditional binary cross-entropy loss function is used, the model tends to predict all samples as negative to obtain a lower overall loss.

[0122] To alleviate this problem and improve the model's ability to discover a small number of positive samples (true regulatory relationships), this invention employs a Focal Loss function during end-to-end training. Focal Loss reduces the weight of easily classified samples (the vast majority of negative samples), prompting the model to focus more on those difficult-to-classify samples (usually sparse positive samples) during training.

[0123] The total loss function is defined as the sum of Focal Losses on all possible gene pairs. The model jointly optimizes all learnable parameters in the Transformer encoder, VAE, Graph Transformer, and bilinear decoder through backpropagation.

[0124] Results evaluation and biological validation: This invention constructs a comprehensive quantitative and qualitative evaluation system to verify the accuracy and biological significance of the inferred results.

[0125] (1) Quantitative indicator evaluation:

[0126] The following metrics were used to evaluate the model's predictive performance on the test set (including retained prior edges and negative samples):

[0127] AUROC (Area Under ROC Curve): The graph uses the false positive rate (FPR) as the x-axis and the true positive rate (TPR) as the y-axis. It measures the model's global classification ability.

[0128] ;

[0129] AUPR (Area Under PR Curve): This metric uses recall as the x-axis and precision as the y-axis. In GRN inference tasks with extreme imbalance between positive and negative samples, AUPR is a more critical metric than AUROC.

[0130] ;

[0131] F1-Score: The harmonic mean of precision and recall.

[0132] ;

[0133] Accuracy:

[0134] .

[0135] (2) Robustness and Ablation Experiment:

[0136] KL divergence analysis: Calculate the KL divergence between the predicted score distributions of positive and negative samples. The larger the value, the stronger the model's discriminative ability.

[0137] Ablation experiments: Transformer channels, AE channels, or parts of the adjacency matrix were removed respectively, and changes in metrics such as AUPR were observed to verify the necessity of each module. The results showed that the combination of dual-channel and multi-adjacency matrix models significantly outperformed single-channel or single-graph models.

[0138] (3) Validation of biological function (GSEA and Hub gene analysis):

[0139] Hub gene identification: In the inferred glioblastoma cell type-specific gene regulatory network, topological indices such as node degree, betweenness centrality, or PageRank are calculated to identify hub transcription factors located at the core of the network. For GBM data, key regulatory factors closely related to tumor stemness maintenance, proliferative activity, and cell state transitions, such as STAT3, SOX2, OLIG2, and ASCL1, should be given special attention.

[0140] Enrichment analysis: Downstream target genes of Hub transcription factors were extracted and KEGG pathway enrichment analysis was performed. Taking STAT3-related target genes as an example, their downstream genes were significantly enriched in the JAK / STAT signaling pathway and the PI3K / AKT-related pathway, suggesting their important role in GBM cell proliferation, invasion, and the shaping of the immunosuppressive microenvironment. Taking SOX2 or OLIG2-related target genes as examples, their downstream genes were enriched in cell cycle regulation, neural stem maintenance, and MAPK-related signaling pathways, thus validating the biological rationale of the model identification results.

[0141] Cell-specific validation: The composition of Hub transcription factors and their downstream regulatory patterns were compared in networks of different cell types. For example, specific transcription factors such as SOX2, OLIG2, or ASCL1, which are related to the maintenance of tumor cell lineage and stemness characteristics, could be identified in malignant tumor cell networks; transcription factors such as SPI1 and IRF8, which are related to immune regulation and inflammatory responses, could be identified in myeloid cell networks; and transcription factors such as ERG and SOX17, which are related to angiogenesis and the maintenance of endothelial cell function, could be identified in vascular-associated cell networks. This validates that the present invention can effectively capture the cell type-specific regulatory characteristics of GBM tumor cells and their microenvironment.

[0142] This method, by systematically integrating single-cell multi-omics data and prior knowledge and utilizing an innovative deep learning architecture, successfully overcomes the limitations of existing technologies, providing strong technical support for the molecular mechanism analysis of glioblastoma, drug target discovery, and precision medicine.

Claims

1. A method for inferring cell-specific regulation in glioblastoma cells, characterized in that, include: The module includes a multi-source data acquisition and preprocessing module, a dual-channel gene feature embedding module, a multi-view graph attention feature update module, a channel feature adaptive fusion module, and a gene regulation relationship prediction output module. The multi-source data acquisition and preprocessing module is used to acquire single-cell sequencing data and prior transcription factor regulatory information of glioblastoma, filter low-quality cells through multi-level quality control, identify tumor cells, immune cells, glial cells and other cell subpopulations by combining dimensionality reduction clustering technology, screen hypervariable genes specific to each cell type, and construct a prior regulatory network framework specific to glioblastoma cell type. The dual-channel gene feature embedding module employs a parallel Transformer channel and an autoencoder AE channel. The Transformer channel utilizes a self-attention mechanism to model long-range dependencies and complex regulatory patterns along the gene-cell relationship axis, capturing the contextual association features of different genes in multicellular states. The autoencoder channel learns the global co-expression structure by reconstructing the cell-gene expression profile and converts the decoder weights into a global embedding representation of the gene; The multi-view graph attention feature update module constructs three complementary graph structures based on prior regulatory information, gene expression association patterns, and functional similarity relationships: a prior regulatory graph, regulatory association graph, and functional similarity graph. It then utilizes a multi-head graph attention mechanism to perform node feature propagation and weighted aggregation on these three graph structures, capturing explicit regulatory relationships, potential association patterns, and functional module information from different perspectives. Finally, it fuses the feature representations from each perspective to obtain biologically interpretable node embedding results. The aforementioned channel feature adaptive fusion module introduces a channel attention mechanism, dynamically calculates dual-channel weights based on the importance of different feature channels to the task, adaptively fuses global structural features and local context features, and generates the final node embedding representation. The gene regulation relationship prediction output module, based on the fused node embedding representation, predicts the probability of connections between TF human transcription factors and target genes through a deep neural network, and constructs the final gene regulation network.

2. The method for inferring cell-specific regulation in glioblastoma as described in claim 1, characterized in that, The multi-source data acquisition and refined preprocessing module specifically includes: (1) Data acquisition: Single-cell RNA sequencing data of tumor tissue and adjacent normal tissue of glioblastoma patients were collected, and the list of human transcription factor TF and its target gene TG regulatory relationship data were downloaded from the database. (2) Single-cell data quality control, including four filtering standards: gene detection number filtering, sequencing depth filtering, mitochondrial gene ratio filtering and ribosomal gene ratio filtering; (3) Data standardization and screening for hypervariable genes, including: Data standardization: The LogNormalize method is used to normalize the filtered expression matrix; Cell annotation: Dimensionality reduction and clustering were performed using the Seurat workflow, and cells were annotated using typical marker genes. Cells were classified into six categories: malignant tumor cells, myeloid cells, T / NK cells, oligodendrocytes, astrocytes, and endothelial cells. Specific feature selection: For each annotated cell type, the variance of gene expression is calculated independently, and the top few genes with the highest degree of variation, i.e. hypervariable genes, are selected as the candidate target gene set. (4) Constructing a cell-specific prior network: The intersection of the list of human transcription factors (TFs) and their target gene TG regulatory relationships with multiple screened hypervariable genes is taken, and the regulatory edges of the target genes belonging to the hypervariable genes of the cell type are retained. Three adjacency matrices are constructed for subsequent module input.

3. The method for inferring cell-specific regulation of glioblastoma cells as described in claim 2, characterized in that, The dual-channel gene feature embedding module employs a parallel dual-channel feature extraction framework, utilizing a Transformer encoder to capture serialized long-range dependency information and a variational autoencoder (VAE) to learn a robust probabilistic global co-expression latent representation, generating an initial gene embedding vector; including: (1) Construct a Transformer encoder channel, which is designed to treat the expression value of a single gene in all cells as a sequence and use the self-attention mechanism to capture complex contextual dependencies and long-distance co-expression patterns across different cell populations; (2) Constructing the variational autoencoder (VAE) channel: This channel aims to learn the global co-expression potential structure of genes from high-noise cell expression data through generative modeling methods, and obtain smoother and more robust feature representations.

4. The method for inferring cell-specific regulation of glioblastoma cells as described in claim 3, characterized in that, The multi-view graph attention feature update module includes: (1) Construct three types of complementary graph structures, including prior regulation graphs, regulation association graphs, and functional similarity graphs; (2) Based on the three types of graph structures, the Graph Transformer feature update explicitly integrates the prior graph structure into the Transformer's self-attention mechanism and achieves deep feature propagation through multi-layer stacking.

5. The method for inferring cell-specific regulation of glioblastoma cells as described in claim 4, characterized in that, The aforementioned channel feature adaptive fusion module adaptively fuses two types of high-level gene feature representations output by the dual-channel gene feature embedding module and the multi-view graph attention feature update module. It dynamically models the feature contribution of the two channels to generate a unified fused feature representation, which serves as the input for subsequent modules. The high-level gene feature representation fuses sequence context features, global co-expression structure probability distribution features, and high-order graph structure features after prior topology propagation.

6. The method for inferring cell-specific regulation of glioblastoma as described in claim 5, characterized in that, The gene regulation relationship prediction output module is a high-level gene feature representation generated by the channel feature adaptive fusion module, which quantitatively predicts whether there is a regulatory relationship between any pair of genes. Specifically: (1) Establish a scoring mechanism based on a bilinear decoder to calculate the original scores of transcription factor-regulated target genes; (2) The original scores are mapped to the (0,1) interval through the Sigmoid activation function. Gene pairs with scores greater than the set threshold are identified as having a regulatory relationship, thereby constructing the final gene regulatory network. (3) Model training: Focal Loss is adopted. The total loss function is defined as the sum of Focal Loss on all possible gene pairs. The model jointly optimizes all learnable parameters in Transformer encoder, VAE, GraphTransformer and bilinear decoder through backpropagation algorithm.

Images