A gene expression gap filling method and system based on deep learning

Abstract

The present invention belongs to the field of electronic digital data processing, and specifically relates to a deep learning-based gene expression gap filling method and system. The method comprises: dividing a detection point into multiple tiles and spots; performing preliminary interpolation of tiles within a spot that lack expression of an expressed gene using a KNN algorithm; selecting tiles from the spot to construct positive and negative samples; training a tile feature embedding module and a spot feature embedding module using a comparative learning module based on the population feature embedding representations of the positive and negative samples; selecting the first one tiles that have the most similar embedding representation of the expressed gene to the tile to be interpolated as similar tiles to the tile to be interpolated; and re-interpolating the expression value of the expressed gene for the tile to be interpolated using the similar tiles to the tile to be interpolated. The present invention improves the local consistency, accuracy, and spatial continuity of the interpolation results.

Description

Technical Field

[0001] The present invention belongs to the field of electrical digital data processing, and in particular relates to a method and system for filling gene expression gaps based on deep learning. Background Art

[0002] Spatial transcriptomics, with its unique advantage of preserving tissue spatial location information while measuring gene expression levels, has become a key tool in biomedical research for exploring tissue microenvironments and cellular heterogeneity. However, due to technical limitations and experimental conditions, spatial transcriptomics data are often subject to gene expression deficits, which not only undermines data integrity but also directly affects the accuracy of subsequent analyses.

[0003] Currently, methods for imputing missing values in spatial transcriptome gene expression mainly fall into three categories: interpolation methods based on statistical models, prediction methods based on machine learning, and generative models based on deep learning. While statistical model-based interpolation methods, such as mean interpolation and K-nearest neighbor interpolation, are simple to implement, they ignore the complex relationships between gene expression and spatial location information, resulting in poor imputation accuracy. Prediction methods based on machine learning, using models such as random forests and support vector machines, have improved imputation accuracy to a certain extent, but have failed to fully exploit the spatial information and local regional interactions of spatial transcriptome data. Generative models based on deep learning, while capable of capturing nonlinear relationships in gene expression, lack consideration of local regional spatial dependencies when processing high-resolution data, resulting in poor spatial continuity in imputation results.

[0004] Existing methods have three significant flaws: first, insufficient utilization of spatial information leads to insufficient spatial accuracy of the filling results; second, they fail to effectively capture the interactions between local regions, resulting in a lack of consistency in local filling results; third, although some deep learning methods have improved the filling accuracy, the models are complex and the computational efficiency is low, making it difficult to adapt to large-scale data sets. Summary of the Invention

[0005] In order to solve the problems existing in the background technology, one aspect of the present invention provides a method for filling gene expression gaps based on deep learning, comprising:

[0006] S1: Preprocessing the raw data of the detection points on the tissue section to obtain high-resolution spatial transcriptome data, wherein the high-resolution spatial transcriptome data includes the gene expression vector of the detection point and the position information of the detection point;

[0007] S2: Divide the detection points into multiple tiles and spots, where each spot contains multiple tiles and each tile contains multiple detection points;

[0008] S3: Genes whose expression frequency in the spot is higher than the set threshold are counted as expressed genes to construct an expression gene set; tiles that lack expression of expressed genes in the spot are preliminarily interpolated using the KNN algorithm;

[0009] S4: Select tiles from the spot to construct positive and negative samples, and calculate the group feature embedding of the positive and negative samples through the tile feature embedding module and the spot feature embedding module; wherein the tile feature embedding module is used to obtain the embedded representation of the tile; and the spot feature embedding module is used to obtain the group feature embedded representation of the positive and negative samples based on the embedded representation of the tile;

[0010] S5: Based on the group feature embedding representation of positive and negative samples, the tile feature embedding module and the spot feature embedding module are trained through the contrastive learning module; the tile embedding representation is extracted through the trained tile feature embedding module;

[0011] S6: The tile that lacks expression of the expressed gene is used as the tile to be interpolated, and the first I tiles that have expression of the expressed gene and are most similar to the embedded representation of the tile to be interpolated are selected as similar tiles to the tile to be interpolated; the expression value of the expressed gene of the tile to be interpolated is re-interpolated using the similar tiles of the tile to be interpolated.

[0012] Another aspect of the present invention provides a gene expression gap filling system based on deep learning, the system comprising a memory and a processor; the memory is used to store an application; the processor is used to run the application and execute the gene expression gap filling method based on deep learning.

[0013] The present invention has at least the following beneficial effects:

[0014] The present invention divides the detection point into multiple tiles and spots. By counting the genes expressed in the spot and processing the tiles, it can effectively capture the interaction between local areas, so that the filling results are more consistent in the local area, and improve the local consistency of the filling results. The method proposed in the present invention avoids the problem of overly complex models to a certain extent through a strategy based on graph neural networks and contrastive learning. While ensuring the accuracy of filling, it can process data more efficiently, achieve a balance between model complexity and efficiency, and is more suitable for large-scale spatial transcriptome data filling. The present invention can more accurately extract the embedded representation of the tile through the calculation of the population feature embedding of positive and negative samples and contrastive learning training, thereby selecting more appropriate similar tiles to re-interpolate the tiles to be interpolated, thereby improving the accuracy of the filling results. At the same time, full consideration of spatial information and the interaction between local areas also improves the spatial continuity of the filling results, solving the shortcomings of existing methods in terms of spatial continuity. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 Schematic diagram of the method of the present invention. DETAILED DESCRIPTION

[0016] The following describes the embodiments of the present invention by means of specific examples, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the illustrations provided in the following embodiments are only schematic illustrations of the basic concept of the present invention, and the following embodiments and features in the embodiments can be combined with each other without conflict.

[0017] See also Figure 1 One aspect of the present invention provides a method for filling gene expression gaps based on deep learning, comprising:

[0018] S1: Preprocessing the original data of the detection points on the tissue section to obtain high-resolution spatial transcriptome data, wherein the high-resolution spatial transcriptome data includes the gene expression vector of the detection point and the position information of the detection point; the original data includes: the original position of the detection point, the position of the tissue section, and the reads Set R;

[0019] Preferably, the preprocessing of the raw data of the detection points on the tissue section includes:

[0020] S11: Use fastqc software to check the data collected at the detection point readsSet R, filtered set R The adapter sequences and low-quality reads in the test point are used to obtain the read set collected at the test point In this example, reads refer to short sequence fragments generated during high-throughput sequencing. These short sequences are the result of sequencing RNA samples. Adapter sequences are non-target sequences introduced during high-throughput sequencing. Adapter sequences are used to fix RNA fragments on sequencing chips. Fastqc software is used in this invention to detect the presence of adapter sequences in reads and filter them. Low-quality reads are determined based on the Phred quality score calculated by Fastqc. The Phred score is used to quantify the probability of sequencing errors and is calculated as: Q = −10 × log 10 (P), where Q is the Phred score, P is the probability of base sequencing errors, and a Phred score of 30 means the base error rate is 0.1%, or 10 -3 In this method, the average Phred quality score of a single read ≥ 30 is defined as a high-quality read, while reads with a score < 30 are defined as low-quality reads and need to be filtered. The judgment criteria of the Phred score can be adjusted according to the actual situation. In this embodiment, only 30 is used as an example.

[0021] S12: Use STAR software to group reads Each read in Aligned to the reference genome Get the comparison result set ;

[0022] S13: Use HTseq software to quantify the target gene set based on the alignment result set M Each gene Count the reads and generate the gene expression vector of the detection point;

[0023] S14: matching the spatial information of the tissue slice with the position coordinates of all detection points, and correcting the position information of each detection point to obtain the final position information of the detection point;

[0024] Preferably, the step S14 includes:

[0025] S141: Based on the spatial information of the tissue slice, a unified spatial coordinate system is determined. For example, a two-dimensional spatial coordinate system is established, with the upper left corner of the tissue slice as the coordinate origin (0, 0), with the horizontal direction defined as the x-axis and the vertical direction as the y-axis. If the tissue slice data is three-dimensional, the z-axis also needs to be defined.

[0026] S142: Mapping the position coordinates of the detection points into the unified spatial coordinate system; for example, converting relative coordinates based on sequencing technology into absolute coordinates based on the tissue section image.

[0027] S143: Using image processing techniques and algorithms, the locations of the detection points are matched to the spatial information of the tissue slice. For example, an image registration algorithm can be used to match the location information of the detection points to specific features on the tissue slice image. Specifically, by identifying landmarks in the tissue slice image (such as cell structures, specific stained areas, etc.), the detection points can be associated with the locations of these landmarks.

[0028] S144: For each detection point, find the corresponding position in the spatial information of the tissue slice;

[0029] S145: Based on the position matching results, the position information of each detection point is corrected. If there is a deviation in the original coordinates of the detection point, the coordinate value of the detection point is adjusted by matching it with the spatial information of the tissue slice. For example, if the actual position of the detection point in the image differs from the original coordinates, the coordinates of the detection point can be corrected according to the coordinate system of the image. Taking into account factors such as possible deformation of the tissue slice during preparation and detection, the position information of the detection point is further optimized and corrected. For example, an elastic registration algorithm can be used to dynamically adjust the position of the detection point to adapt to the deformation of the tissue slice.

[0030] S146: Verify the corrected detection point position information to ensure its accuracy and reliability. This can be done by comparing it with known spatial information on the tissue slice to check whether the detection point positions are appropriate. The corrected position information is evaluated using specific metrics, such as position error and matching accuracy. If the evaluation results are unsatisfactory, return to the previous steps to adjust and optimize the position matching and correction process.

[0031] S15: Use the TPM method to normalize the gene expression vector of the detection point to obtain the gene expression vector of the detection point after normalization.

[0032] Preferably, the step S15 includes:

[0033] S151: For target gene sets Each gene , obtain its gene length information; gene length can be obtained by querying gene databases or reference genome annotation files. Generally, gene length refers to the nucleotide sequence length from the gene's start codon to the stop codon, measured in base pairs (bp).

[0034] S152: Target gene set is generated using HTseq software according to step S13 Each gene The result of counting the reads is to get the number of reads for each gene Here Indicates that at this detection point, the gene The number of times it was sequenced.

[0035] S153: For each gene , calculate the number of reads per kilobase , the calculation formula is:

[0036]

[0037] in, It's genes Length (in kilobases); The value reflects the gene The expression level per unit length eliminates the effect of gene length on read counts.

[0038] S154: Add up the RPK values of all genes to get the total RPK value of the detection point ,in, is the total number of genes in the target gene set Genes.

[0039] S155: For each gene , calculate its TPM value:

[0040]

[0041] in, Indicates gene TPM value;

[0042] The TPM value represents the number of transcripts per million for a gene and is a standardized gene expression indicator that can be used to compare gene expression levels at different detection points or between different samples.

[0043] Through the above steps, the TPM method can be used to standardize the gene expression vectors of the detection points to obtain the gene expression vectors of the detection points after standardization, so that the expression levels of different genes are comparable, providing a unified standard for subsequent data analysis and comparison.

[0044] In this embodiment, in step S11, the raw data is comprehensively tested using fastqc software, low-quality reads and adapter sequences are accurately filtered, and interference factors such as sequencing errors and contamination are effectively eliminated to obtain a pure reads set R1, thereby building a reliable data foundation for subsequent analysis; in step S12, the reads in R1 are accurately aligned with the reference genome Gref with the help of STAR software, and the correspondence between the sequencing data and the genome is established, and the position of the reads in the genome is clarified, which not only provides a basis for subsequent gene expression calculations, but also can mine gene structure and variation information, thereby enhancing the biological significance of the data. In step S13, the reads of the target gene set Genes are counted using HTseq software, the original sequencing data is converted into an intuitive gene expression vector, the expression level of each gene at the detection point is quantified, and core data is provided for gene expression pattern analysis; in step S14, the spatial information of the tissue section is matched with the detection point coordinates to correct the position information, so that the gene expression data is accurately associated with the spatial structure of the tissue, laying the foundation for studying the spatial distribution of gene expression and analyzing the tissue microenvironment; in step S15, the TPM method is used to standardize the gene expression vector to eliminate the deviation caused by the difference in sequencing depth of different samples, ensure the comparability of gene expression data between different detection points, and enhance the analytical value of the data.

[0045] After this series of preprocessing, the resulting high-resolution spatial transcriptome data contains both precise gene expression information and reliable location information, providing high-quality data support for subsequent operations such as tile and spot division, gene expression statistics, construction of positive and negative samples, and comparative learning training. This not only significantly improves the accuracy of missing gene expression values but also enhances the spatial continuity of the imputed results. It also effectively balances model complexity and computational efficiency, creating favorable conditions for the in-depth analysis and application of large-scale spatial transcriptome data.

[0046] S2: Divide the detection points into multiple tiles and spots, where each spot contains multiple tiles and each tile contains multiple detection points. A tile is a square detection area containing N×N detection points, and a spot is a square area containing M×M tile areas, where N≥10 and M≥10.

[0047] Preferably, the dividing the detection points into multiple tiles and spots includes: the detection points are distributed in an array on the tissue slice to form an M×M two-dimensional detection point grid area, each grid representing a detection point; the M×M two-dimensional grid area is divided into N×N tiles of size A×A, where N=M÷A; each tile is used as a grid to form an N×N two-dimensional tile grid area, and the N×N two-dimensional tile grid area is divided into K×K spots of size B×B; where K=N÷B, M, N, A, B and K are all positive integers; the gene expression vectors of all detection points in each tile are added to obtain the gene expression vector of the tile; the position of each tile is the center value of all detection point positions in the tile.

[0048] In this example, the M×M two-dimensional grid area is first divided into N×N tiles of size A×A, where N=M÷A. This division is equivalent to a preliminary "gridding subdivision" of the entire detection area, clustering the previously scattered detection points into groups of A×A points. Each tile becomes a basic unit containing multiple detection points. For example, if M=100 and A=10, then 10×10 tiles are formed, each containing 10×10=100 detection points. This division facilitates the subsequent integration and analysis of gene expression features in the local area. Next, each tile is treated as a new "grid," constructing an N×N two-dimensional tile grid area. This is then further divided into K×K spots of size B×B, where K=N÷B. This step further aggregates the tile units, forming a larger set of local areas. Continuing with the previous example, if B=2, then N=10 and K=5, the resulting number is 5×5 spots, each containing 2×2 tiles. This hierarchical partitioning strategy, from detection points to tiles and then to spots, constructs a spatial region system from small to large, from microscopic to macroscopic, achieving a structured decomposition and integration of spatial information in tissue sections. The gene expression vector for each tile is obtained by summing the gene expression vectors of all detection points within the tile. This operation effectively aggregates the gene expression information of each detection point within the tile. On the one hand, the summation operation highlights the overall trend of gene expression within the local region, reducing the impact of individual detection point data fluctuations on the overall characteristics. On the other hand, compressing the gene expression vectors of multiple detection points into a single tile gene expression vector significantly reduces data dimensionality, reduces the complexity of subsequent calculations, and improves model processing efficiency. Furthermore, the position of each tile is the center value of all detection points within the tile. This setting assigns clear and representative spatial coordinates to the tile, allowing accurate consideration of the spatial relationship between tiles in subsequent analysis and providing a precise positional reference for imputing missing gene expression values using spatial information. This partitioning approach offers multiple advantages.First, by dividing the detection points into tiles and spots, the spatial location information of the spatial transcriptome data is fully utilized, so that in the subsequent process of filling missing gene expression values, analysis and calculation can be performed based on the spatial structural characteristics of the local area, thereby improving the spatial accuracy of the filling results; second, the constructed tile and spot hierarchy system strengthens the interaction between local areas. When performing gene expression analysis and missing value interpolation, the collaborative relationship within and between local areas can be better considered to ensure the consistency of the filling results in local areas; third, by aggregating gene expression vectors and regularizing position information, the model complexity and efficiency are balanced to a certain extent, which can not only retain key spatial and gene expression features but also reduce data redundancy, so that models based on graph neural networks and contrastive learning can process data more efficiently and are more suitable for the analysis and application of large-scale spatial transcriptome datasets.

[0049] S3: Genes whose expression frequency in the spot is higher than the set threshold are counted as expressed genes to construct an expression gene set; tiles that lack expression of expressed genes in the spot are preliminarily interpolated using the KNN algorithm;

[0050] Preferably, the KNN algorithm is used to preliminarily interpolate tiles within a spot that lack expression of the expressed gene, including: selecting the tile lacking expression of the expressed gene as the tile to be interpolated; selecting, based on the location of the tile to be interpolated, the K0 tiles closest to the tile to be interpolated that express the expressed gene; and using the mean expression value of the expressed gene in the K0 tiles closest to the tile to be interpolated that express the expressed gene as the expression value of the expressed gene in the tile to be interpolated, to preliminarily interpolate the tile to be interpolated. The value of K0 ranges from 0.05×M×M to 0.1×M×M, where M×M is the number of tiles contained in the spot.

[0051] In the present embodiment, in each spot, the gene expression frequency is counted and the genes above the set threshold are screened to construct an expression gene set. This operation focuses on identifying genes with active expression characteristics in local areas. The purpose of setting the threshold is to filter out low expression, gene signals that may be caused by experimental noise or accidental factors, and retain genes that really function in the area. For example, if the threshold is set to 0.8, it means that in a spot, only when a gene is expressed in 80% and above tiles, it will be included in the expression gene set. This screening mechanism effectively reduces data redundancy, highlights key gene expression information, and enables subsequent analysis to be carried out around core genes, avoiding the impact of secondary or interfering gene expression data on the filling results. At the same time, constructing the expression gene set clarifies the target gene range for subsequent interpolation operations, making the interpolation process more targeted and improving computational efficiency.

[0052] After marking the tiles within the spot that lack expression of the target gene as tiles to be interpolated, the KNN algorithm works based on spatial location information. Based on the location of the tile to be interpolated, the algorithm searches for the K0 tiles closest to it within the same spot that have expression of the corresponding gene. The "distance" here is usually calculated using metrics such as Euclidean distance to quantify the spatial proximity between tiles. For example, if K0 is set to 5, the five tiles closest to the tile to be interpolated that have expression of the target gene are selected. Subsequently, the mean expression value of the target gene of these K0 tiles is used as the expression value of the tile to be interpolated for preliminary interpolation. This method is based on the spatial proximity assumption and assumes that tiles with similar spatial locations have similar gene expression patterns. Using the expression information of surrounding tiles to provide reference values for the tile to be interpolated can preliminarily fill in the missing gene expression data.

[0053] From the perspective of the overall technical process, the dual operation of step S3 brings significant results. At the data feature extraction level, the construction of the expression gene set realizes the feature screening and dimensionality reduction of spatial transcriptome data, focusing on core gene expression information, and providing purer and more representative data input for subsequent comparative learning and model training. At the missing value processing level, the preliminary interpolation of the KNN algorithm is based on the spatial position relationship, effectively utilizing the spatial dependence and expression similarity between tiles in the local area, and to a certain extent alleviates the problem of missing gene expression. Although the preliminary interpolation results may be insufficient in accuracy, they lay the foundation for subsequent fine-tuning through comparative learning and deep learning models, reducing the difficulty and computational complexity of subsequent model training. At the same time, it also provides a more complete data foundation for subsequent complex analysis based on tile and spot features, improving the consistency of the filling results in local areas and the accuracy in space, which helps to achieve more accurate filling of missing gene expression values.

[0054] S4: Select tiles from the spot to construct positive and negative samples, and calculate the group feature embedding of the positive and negative samples through the tile feature embedding module and the spot feature embedding module; wherein the tile feature embedding module is used to obtain the embedded representation of the tile; and the spot feature embedding module is used to obtain the group feature embedded representation of the positive and negative samples based on the embedded representation of the tile;

[0055] Preferably, constructing positive and negative samples includes: taking all tiles of the spot as first positive samples; selecting tiles in an 80%-90% square area within the spot as second positive samples;

[0056] Select C×C tiles at the four corners of the spot, and randomly select P tiles in each C×C tile as negative samples, where , .

[0057] Preferably, the tile feature embedding module adopts a first transformer encoder, and obtaining the embedded representation of the tile includes: expanding the gene expression vector of the interpolated tile into a two-dimensional matrix to adapt to the input format of the first transformer encoder, and mapping the gene expression vector expanded into a two-dimensional matrix to a high-dimensional space; using Sinusoidal to encode the position coordinates of the tile, and splicing the position encoding features of the tile and the gene expression vector mapped to the high-dimensional space to obtain the input features of the first transformer encoder, and encoding the input features through the first transformer encoder to obtain the embedded representation of the tile.

[0058] Preferably, the spot feature embedding module includes: a graph structure construction module, a GAT encoding module and a deep feature extraction module; the graph structure construction module is used to construct a graph structure according to the embedded representation and position of the tile in the positive and negative samples; the GAT encoding module uses multiple GAT layers to aggregate and update the feature matrix of the graph structure according to the adjacency matrix of the graph structure; the deep feature extraction module includes: a second transformer encoder, an attention pooling layer, a global pooling layer and a fully connected layer; the input features of the second transformer encoder are the feature matrix after aggregation and update of the GAT encoding module; the input features of the second transformer encoder are respectively used as the input features of the attention pooling layer and the global pooling layer; the output features of the attention pooling layer and the global pooling layer are spliced as the input features of the fully connected layer; the output features of the fully connected layer are the group feature embedding representation of the positive and negative samples.

[0059] Preferably, constructing a graph structure according to the embedded representation and position of tiles in positive and negative samples includes:

[0060] S41: Use the tiles in the positive and negative samples as nodes and the embedded representation of the tiles as the features of the tile nodes to construct the initial feature matrix of the graph structure;

[0061] S42: Connect each tile to its K1 nearest neighbor tiles, and calculate the local weight of the edge based on the positions of the tile and its connected neighbor tiles; where K1 is an integer ranging from 1 / 4×M to 1 / 2×M; calculating the local weight of the edge includes calculating the Euclidean distance between the positions of the tile and its connected neighbor tiles, and performing normalization to obtain the local weight of the edge. The local weight calculation formula for the edge is: ,in, represents the local weight of the edge between tile i and tile j, is the position distance between two connected tiles i and j, and max_distance is the maximum distance used for normalization.

[0062] S43: Connect each tile to the K2 similar tiles with the most similar features, and calculate the global weight of the edge based on the features of the tile and the similar tiles it is connected to; where the value range of K2 is: 0.05×M×M~0.1×M×M, where M×M is the number of tiles contained in the spot; calculating the global weight of the edge includes: calculating the similarity between the features of the tile and the similar tiles it is connected to, and performing normalization processing to obtain the global weight of the edge. The global weight calculation formula of the edge is: ,in is the position distance between two connected tiles i and j, and is the gene expression vector (feature) between the two, max_distance is the maximum distance used for normalization, is the cosine similarity between tile i and tile j; represents the global weight of the edge between tile i and tile j;

[0063] S44: Remove the duplicated edges between tiles and retain the maximum value of the local weight and the global weight as the final weight of the edge, and then construct the adjacency matrix of the graph structure; the final weight calculation formula of the edge is: , is the final weight of the edge between tile i and tile j; construct the adjacency matrix A of the graph structure: if there is an edge between node tile i and tile j, then A ij =W ij , otherwise A ij =0;A ij Represents the element in row i and column j of the adjacency matrix A;

[0064] In this embodiment, the construction of positive and negative samples adopts a hierarchical and targeted selection method. All tiles of the spot are used as the first positive sample, and tiles of 80%-90% square area in the spot are selected as the second positive sample. This selection method ensures that the positive sample covers the main gene expression information and spatial structural characteristics in the spot, and can represent the common pattern of gene expression in the spot. Tiles of specific size areas are selected at the four corners of the spot and negative samples are randomly extracted. The spatial difference between the corner area and the central area is used to introduce samples with different expression patterns; it is ensured that the negative samples and positive samples are effectively distinguished in spatial position and expression characteristics. The reasonable construction of positive and negative samples provides rich and contrasting data samples for comparative learning, so that the model can learn the difference characteristics of different gene expression patterns and spatial structures, thereby enhancing the model's ability to discriminate the characteristics of spatial transcriptome data.

[0065] In this embodiment, the tile feature embedding module adopts the first transformer encoder, and its processing process fully integrates gene expression information and spatial position information. First, the gene expression vector of the interpolated tile is expanded to a two-dimensional matrix. This operation not only adapts to the input format of the transformer encoder, but also provides a structural basis for the subsequent mining of potential features of gene expression data. Then, by mapping the gene expression vector to a high-dimensional space, the complex nonlinear relationship between gene expression can be captured and richer expression features can be mined. The position coordinates of the tile are encoded using Sinusoidal, the position information is converted into a computable feature, and spliced with the gene expression feature, so that the model can simultaneously consider gene expression and spatial position information during the learning process, avoiding the defect of existing methods that ignore spatial information. Finally, after encoding by the first transformer encoder, the tile embedding representation obtained contains the dual features of gene expression and spatial position, providing a more comprehensive and accurate feature input for the calculation and comparative learning of the subsequent spot feature embedding module.

[0066] In this embodiment, the spot feature embedding module is composed of a graph structure construction module, a GAT encoding module, and a deep feature extraction module working together. The graph structure construction module constructs an initial feature matrix by using tiles as nodes and embedding representations as node features, connects nodes based on position and features, and calculates edge weights. The constructed graph structure not only reflects the spatial proximity relationship between tiles, but also reflects the similarity relationship of gene expression features. The GAT encoding module uses multiple GAT layers to aggregate and update the feature matrix based on the adjacency matrix of the graph structure, effectively capturing the interaction between tiles in the local area, making the feature representation more locally consistent. In the deep feature extraction module, the second transformer encoder further extracts features, combines the advantages of the attention pooling layer and the global pooling layer, screens and integrates the features, and finally outputs the group feature embedding representation of positive and negative samples through the fully connected layer. This multi-level feature extraction and integration mechanism fully exploits the characteristics of the tile groups in the positive and negative samples, retains the detailed features of the local area, and extracts the overall global features, providing high-quality feature representation for contrastive learning, helping the model to learn the essential differences between positive and negative samples, and improving the model's understanding and processing capabilities of spatial transcriptome data features.

[0067] Through the construction of positive and negative samples and the calculation of feature embeddings, a solid foundation is laid for contrastive learning at both the data sample and feature representation levels. At the data level, properly constructed positive and negative samples provide rich and discriminative training material for contrastive learning, guiding the model to learn the characteristic differences between different gene expression patterns and spatial structures. At the feature level, the tile and spot feature embedding modules deeply integrate spatial location and gene expression information. The extracted population feature embedding representation accurately reflects the spatial and expression characteristics of the sample, enabling the model to more accurately learn the associations and differences between features in subsequent contrastive learning training. This effectively addresses the problems of existing methods in insufficient utilization of spatial information and inadequate capture of local regional interactions, providing strong support for improving the accuracy and spatial continuity of gene expression missing value filling results. It also balances model complexity and efficiency, making it suitable for processing large-scale spatial transcriptome data.

[0068] S5: Based on the group feature embedding representation of positive and negative samples, the tile feature embedding module and the spot feature embedding module are trained through the contrastive learning module; the tile embedding representation is extracted through the trained tile feature embedding module;

[0069] Preferably, the contrastive learning module uses a noise contrast estimation function InfoNCE Loss:

[0070]

[0071] Where log represents the natural logarithm, exp represents the exponential function, L represents the loss function, q represents the population feature embedding representation of the first positive sample, k+ represents the population feature embedding representation of the second positive sample, k− represents the population feature embedding representation of the negative sample, τ is the temperature parameter; sim represents the cosine similarity function.

[0072] S6: The tile that lacks expression of the expressed gene is used as the tile to be interpolated, and the first I tiles that have expression of the expressed gene and are most similar to the embedded representation of the tile to be interpolated are selected as similar tiles to the tile to be interpolated; the expression value of the expressed gene is re-interpolated for the tile to be interpolated using the similar tiles of the tile to be interpolated. The value range of I is 0.05×M×M~0.1×M×M, where M×M is the number of tiles contained in the spot.

[0073] In the previous steps of this embodiment, the detection points on the tissue section have been divided into tiles and spots, and the expression gene set has been counted. In this step, the tiles that lack expression of the expression genes are marked as tiles to be interpolated. These tiles may not detect the expression of certain expression genes in the original data due to technical limitations, experimental errors, etc. In order to interpolate the tiles to be interpolated, it is necessary to find tiles similar to them. The embedding representation of the tile is used here to measure the similarity. In step S4, the embedding representation of each tile has been calculated by the tile feature embedding module and the spot feature embedding module. These embedding representations contain various features such as the gene expression information and spatial position information of the tile. From all tiles that have expression of the expression gene, the first I tiles that are most similar to the embedding representation of the tile to be interpolated are selected as similar tiles to the tile to be interpolated. Common distance metrics such as Euclidean distance and cosine similarity can be used to measure similarity. For example, if cosine similarity is used, the higher the similarity, the more similar the embedding representations of the two tiles are. The expression values of the genes expressed in the tile to be interpolated are re-interpolated using the expression values of the genes expressed in the tiles similar to the tile to be interpolated. A common interpolation method is to calculate the mean expression value of the genes expressed in these tiles and use this mean as the expression value of the genes expressed in the tile to be interpolated.

[0074] By selecting similar tiles for interpolation based on a tile-based embedding representation, we fully leverage gene expression and spatial location information. This embedding representation, meticulously calculated in the previous steps, comprehensively reflects the various characteristics of each tile. Compared to traditional interpolation methods based solely on gene expression values or spatial location, this approach more accurately identifies tiles that are similar to the tile to be interpolated in both gene expression patterns and spatial distribution, thereby improving interpolation accuracy. Because the embedding representation incorporates spatial location information, the selected similar tiles are typically spatially close to the tile to be interpolated. This results in more spatial continuity after interpolation, avoiding obvious discontinuities or abrupt changes. For example, in a tissue section, the interpolated gene expression values of adjacent tiles will better align with actual biological patterns, resulting in a more reasonable gene expression distribution across the entire tissue section. Step S6 builds on the previous steps, which include data preprocessing, tile and spot partitioning, gene expression set construction, and feature embedding calculation, all of which provide strong support for this interpolation step. Through this step-by-step approach, the quality of imputing missing gene expression values is gradually improved, ensuring that the final imputation results better reflect the true gene expression status of the tissue section. The number of similar tiles selected, I, can be adjusted based on the specific dataset and analysis requirements. When the data volume is large and the gene expression patterns are complex, the value of I can be appropriately increased to obtain more similar information. When the data volume is small or local features are of greater interest, the value of I can be reduced, so that the imputation results are more focused on the few tiles that are most similar to the tile to be imputed. This flexibility makes this method adaptable to different types of spatial transcriptome data.

[0075] Another aspect of the present invention provides a gene expression gap filling system based on deep learning, the system comprising a memory and a processor; the memory is used to store an application; the processor is used to run the application and execute the gene expression gap filling method based on deep learning.

[0076] Those skilled in the art will appreciate that all or part of the processes in the above-mentioned embodiment methods can be implemented by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to memory, storage, database or other media used in the embodiments provided in this application may include non-volatile and / or volatile memory. Non-volatile memory may include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory may include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in many forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous 8DRAM (SDRAM), double data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct 10RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0077] In summary, the present invention divides the detection point into multiple tiles and spots. By counting the genes expressed in the spot and processing the tiles, it can effectively capture the interaction between local areas, so that the filling results are more consistent in the local area, and improve the local consistency of the filling results. The method proposed in the present invention avoids the problem of overly complex models to a certain extent through a strategy based on graph neural networks and contrastive learning. While ensuring the accuracy of filling, it can process data more efficiently, achieve a balance between model complexity and efficiency, and is more suitable for large-scale spatial transcriptome data filling. The present invention can more accurately extract the embedded representation of the tile through the calculation of the population feature embedding and contrastive learning training of positive and negative samples, thereby selecting more appropriate similar tiles to re-interpolate the tiles to be interpolated, thereby improving the accuracy of the filling results. At the same time, full consideration of spatial information and the interaction between local areas also improves the spatial continuity of the filling results, solving the shortcomings of existing methods in terms of spatial continuity.

[0078] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not limiting. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solutions, which should all be included in the scope of the claims of the present invention.

Claims

1. A gene expression gap filling method based on deep learning, characterized in that: include: S1: Preprocessing the raw data of the detection points on the tissue section to obtain high-resolution spatial transcriptome data, wherein the high-resolution spatial transcriptome data includes the gene expression vector of the detection point and the position information of the detection point; S2: Divide the detection points into multiple tiles and spots, where each spot contains multiple tiles and each tile contains multiple detection points; The method of dividing the detection points into a plurality of tiles and spots includes: distributing the detection points in an array on the tissue slice to form an M×M two-dimensional detection point grid area, where each grid represents a detection point; dividing the M×M two-dimensional grid area into N×N tiles of size A×A, where N=M÷A; using each tile as a grid to form an N×N two-dimensional tile grid area, and dividing the N×N two-dimensional tile grid area into K×K spots of size B×B, where K=N÷B, and M, N, A, B, and K are all positive integers; summing the gene expression vectors of all detection points in each tile to obtain the gene expression vector of the tile; and taking the position of each tile as the center value of all detection point positions in the tile; S3: Genes whose expression frequency in the spot is higher than the set threshold are counted as expressed genes to construct an expression gene set; tiles that lack expression of expressed genes in the spot are preliminarily interpolated using the KNN algorithm; S4: Select tiles from the spot to construct positive and negative samples, and calculate the group feature embedding of the positive and negative samples through the tile feature embedding module and the spot feature embedding module; wherein the tile feature embedding module is used to obtain the embedded representation of the tile; and the spot feature embedding module is used to obtain the group feature embedded representation of the positive and negative samples based on the embedded representation of the tile; The constructing of positive and negative samples includes: taking all tiles of the spot as the first positive sample; selecting tiles of 80%-90% of the square area in the spot as the second positive sample; Select C×C tiles at the four corners of the spot, and randomly select P tiles in each C×C tile as negative samples, where , ; S5: Based on the group feature embedding representation of positive and negative samples, the tile feature embedding module and the spot feature embedding module are trained through the contrastive learning module; the tile embedding representation is extracted through the trained tile feature embedding module; S6: The tile that lacks expression of the expressed gene is used as the tile to be interpolated. The first I tiles that have expression of the expressed gene and are most similar to the embedded representation of the tile to be interpolated are selected as similar tiles to the tile to be interpolated. The expression value of the expressed gene of the tile to be interpolated is re-interpolated using the similar tiles of the tile to be interpolated. The value range of I is 0.05×M×M~0.1×M×M.

2. A method for filling gene expression gaps based on deep learning according to claim 1, characterized in that: The preprocessing of the raw data of the detection points on the tissue slice comprises: S11: Use fastqc software to check the read set R collected at the detection point, filter the adapter sequences and low-quality reads in set R to obtain the read set collected at the detection point ; S12: Use STAR software to group reads Each read in Aligned to the reference genome Get the comparison result set ; S13: Use HTseq software to quantify the target gene set based on the alignment result set M Each gene Count the reads and generate the gene expression vector of the detection point; S14: matching the spatial information of the tissue slice with the position coordinates of all detection points, and correcting the position information of each detection point to obtain the final position information of the detection point; S15: Use the TPM method to normalize the gene expression vector of the detection point to obtain the gene expression vector of the detection point after normalization.

3. The method for filling gene expression gaps based on deep learning according to claim 1, characterized in that: The contrastive learning module uses the noise contrast estimation function InfoNCE Loss: Among them, log represents the natural logarithm, exp represents the exponential function, and L represents the loss function. represents the group feature embedding representation of the first positive sample, k+ represents the group feature embedding representation of the second positive sample, k− represents the group feature embedding representation of the negative sample, τ is the temperature parameter; sim represents the cosine similarity function.

4. The method for filling gene expression gaps based on deep learning according to claim 3, characterized in that: The preliminary interpolation of tiles lacking expression of the expressed gene in the spot by the KNN algorithm includes: taking the tile lacking expression of the expressed gene as the tile to be interpolated; selecting K0 tiles with expression of the expressed gene that are closest to the tile to be interpolated according to the position of the tile to be interpolated; and using the average expression amount of the expressed gene of the K0 tiles with expression of the expressed gene that are closest to the tile to be interpolated as the expression value of the expressed gene of the tile to be interpolated, and performing preliminary interpolation on the tile to be interpolated; the range of K0 values is: 0.05×M×M~0.1×M×M.

5. The method for filling gene expression gaps based on deep learning according to claim 4, characterized in that: The tile feature embedding module adopts a first transformer encoder to obtain the embedded representation of the tile, including: expanding the gene expression vector of the interpolated tile into a two-dimensional matrix to adapt to the input format of the first transformer encoder, and mapping the gene expression vector expanded into a two-dimensional matrix to a high-dimensional space; using Sinusoidal to encode the position coordinates of the tile, and splicing the position encoding features of the tile and the gene expression vector mapped to the high-dimensional space to obtain the input features of the first transformer encoder, and encoding the input features through the first transformer encoder to obtain the embedded representation of the tile.

6. The method for filling gene expression gaps based on deep learning according to claim 1, characterized in that: The spot feature embedding module includes: a graph structure construction module, a GAT encoding module and a deep feature extraction module; the graph structure construction module is used to construct a graph structure according to the embedded representation and position of the tile in the positive and negative samples; the GAT encoding module uses multiple GAT layers to aggregate and update the feature matrix of the graph structure according to the adjacency matrix of the graph structure; the deep feature extraction module includes: a second transformer encoder, an attention pooling layer, a global pooling layer and a fully connected layer; the input features of the second transformer encoder are the feature matrix after aggregation and update of the GAT encoding module; the input features of the second transformer encoder are respectively used as the input features of the attention pooling layer and the global pooling layer; the output features of the attention pooling layer and the global pooling layer are spliced as the input features of the fully connected layer; the output features of the fully connected layer are the group feature embedding representation of the positive and negative samples.

7. The method for filling gene expression gaps based on deep learning according to claim 6, wherein constructing a graph structure based on the embedded representation and position of tiles in positive and negative samples comprises: S41: Use the tiles in the positive and negative samples as nodes and the embedded representation of the tiles as the features of the tile nodes to construct the initial feature matrix of the graph structure; S42: Connect each tile to its K1 nearest neighbor tiles and calculate the local weight of the edge based on the positions of the tile and its connected neighbor tiles; S43: Connect each tile to the K2 similar tiles with the most similar features, and calculate the global weight of the edge based on the features of the tile and the similar tiles it is connected to; S44: Remove the duplicated edges between tiles and retain the maximum value of the local weight and the global weight as the final weight of the edge, thereby constructing the adjacency matrix of the graph structure; The value range of K1 is: 1 / 4×M~1 / 2×M; The value range of K2 is: 0.05×M×M~0.1×M×M.

8. A deep learning-based gene expression gap filling system, the system comprising a memory and a processor; the memory is used to store an application; the processor is used to run the application and execute a deep learning-based gene expression gap filling method as described in any one of claims 1 to 7.

Images