Single cell sequencing data interpolation method based on biclustering and graph convolutional neural network

Abstract

A single cell sequencing data interpolation method based on a biclustering and graph convolutional neural network comprises the following steps: 1) respectively downloading matrix-form gene expression data of human and mouse from a gene database, and screening out cells of which the total number of gene expressions is less than 3 and genes of which the total number of cell expressions is less than 200; 2) performing feature selection on the gene expression matrix, finding out cells with consistent growth trend, forming a new two-dimensional matrix, and constructing a cell connection diagram; 3) constructing a BiGCN single-cell interpolation model by using a graph convolutional neural network, learning a nonlinear relationship between cells, dividing a data set into a training set and a test set, training the model by using the training set, and reserving the trained BiGCN single-cell interpolation model; 4) inputting the test set into the trained BiGCN single cell interpolation model to obtain an interpolation result, and evaluating the interpolation result by using an evaluation index; the method has the characteristics that the feature extraction capability is improved, interference of redundant zero expression is avoided, data sparsity is effectively relieved, and the interpolation accuracy is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of single-cell sequencing data interpolation, and in particular relates to a single-cell sequencing data interpolation method based on biclustering and graph convolutional neural network. Background Art

[0002] Single-cell sequencing is a technology that amplifies and sequences genomes and transcriptomes at the level of individual cells. It plays a vital role in fields such as immunology, neuroscience, and developmental biology. Using single-cell sequencing technology, we can not only effectively identify cell subpopulations within a sample and identify different types of cells, but also explore intercellular heterogeneity and reveal the characteristics of tumor cells, which is of great help in exploring the pathogenesis of the human body. With the vigorous development of single-cell sequencing, current sequencing technology can detect a large number of samples, and the amount of data obtained is becoming increasingly large. However, in the process of sequencing library construction, it is easy for some data to be missing due to amplification failure, which is also called the missing phenomenon. In order to solve the problem of missing data, some interpolation methods have been proposed, but there are still some problems, such as false positives caused by over-interpolation, data smoothing caused by over-training, and the inability to fully learn the nonlinear relationship between cells.

[0003] Most of the current interpolation models based on graph neural networks construct cell connection graphs based on K nearest neighbors. This method simply selects the nearest K cells based on the Euclidean distance to form a cell graph. Due to the particularity of biological samples, cells of the same type may be too far apart when cells are obtained. As the data set becomes larger and larger, the distance between cells of the same type in the matrix will become more dispersed, making it difficult to obtain relevant connections between cells. Summary of the invention

[0004] In order to overcome the shortcomings of the above-mentioned prior art, the purpose of the present invention is to provide a single-cell sequencing data interpolation method based on biclustering and graph convolutional neural networks. Aiming at the problem of insufficient learning caused by the nonlinear relationship between cells in the interpolation process, a single-cell interpolation model BiGCN (full name Bi-clustering GraphConvolution Network, biclustering graph convolutional neural network) based on biclustering and graph convolutional neural networks is designed. The biclustering is used to simultaneously cluster the rows and columns of the expression matrix to obtain similar cell subpopulations, and the graph convolutional neural network is fused with the autoencoder to replace the linear layer with the graph convolution layer to further explore the relationship between cells and genes. It has the characteristics of improving feature extraction ability, avoiding the interference of redundant zero expression, effectively alleviating data sparsity, and improving interpolation accuracy.

[0005] To achieve the above object, the technical solution adopted by the present invention is: a single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network, comprising the following steps:

[0006] Step 1, downloading the gene expression data of humans and mice in matrix form from the gene database, and filtering out cells with a total gene expression of less than 3 and genes with a total cell expression of less than 200;

[0007] Step 2: Select features of the gene expression matrix and use biclustering to find cells with consistent growth trends to form a new two-dimensional matrix. Use cells as nodes and cell correlations as features to build a cell connection map.

[0008] Step 3: Use graph convolutional neural network to build BiGCN single-cell interpolation model, learn nonlinear relationships between cells, divide the data set into training set and test set, use the training set to train the model, and retain the trained BiGCN single-cell interpolation model;

[0009] Step 4: Input the test set into the trained BiGCN single-cell interpolation model, and evaluate the interpolation results using evaluation indicators.

[0010] The step 2, using double clustering to find cells with consistent gene expression level change trends, specifically includes the following steps:

[0011] Step 2-1, preprocessing stage

[0012] Divide the data matrix into single-row and single-column sub-biclusters, which together represent all rows and columns of the expression matrix, and each bicluster contains a row set and a column set, from which some columns and rows are randomly selected to form the index matrix of the starting bicluster;

[0013] Step 2-2, find associated cells

[0014] The constructed index matrix can reflect the changing trend of the expression value. The binary joint information entropy is used to determine the correlation between two cells. Formula (1) is selected

[0015]

[0016] In the formula, X represents the gene expression intensity, m represents the total number of cells, n represents the total number of genes, log represents the logarithmic function with e as the base, and x i Represents a cell in a row, y j represents a gene in a column, p(x i ,y j ) represents the probability density between cells and genes, H(X) represents the degree of association between the cell and the gene, and the larger the value, the higher the degree of association;

[0017] Step 2-3, iterative optimization

[0018] The binary joint information entropy of the double clustering is optimized in an iterative manner. In each iteration, a row is randomly added or deleted to improve the binary joint information entropy of the cell. The iterative process of seeking the maximum binary joint information entropy uses greedy search, trying to maximize the binary joint information entropy at each step. When the value no longer changes, the iteration is stopped.

[0019] The step 3, constructing the BiGCN single-cell interpolation model, specifically includes the following steps:

[0020] The BiGCN single-cell interpolation model is improved on the basis of the traditional autoencoder. In order to further capture the relationship between cells, the graph convolution layer is used instead of the linear layer to maximize the learning of the characteristics between cells. Before entering the graph convolutional neural network, the expression matrix is ​​converted into a cell connection graph, and then the ZINB loss function is used to compare the matrices before and after interpolation to reduce the error.

[0021] The autoencoder has a four-layer structure, the first and third layers are graph convolution layers, the second and fourth layers are Dropout layers. In order to avoid overfitting problems in the autoencoder, the number of convolution layers used cannot be too deep, otherwise it will cause data smoothing; the input layer and output layer of the autoencoder are both 64 layers; the decoder part has three fully connected layers, which respectively represent the mean, variance and sparsity of the model, representing the sparsity of zero values ​​in the current sample, and using the negative binomial distribution result as the loss function.

[0022] The BiGCN single-cell interpolation model described above has a large sparsity of single-cell sequencing data and cannot be fitted using a simple distribution. Since the zero value distribution in the data is not uniform, it is not suitable to use the Poisson distribution. The negative binomial distribution is currently the closest one. The present invention uses a zero-inflated negative binomial distribution based on a modified version of the negative binomial distribution. The results of using it are improved to a certain extent compared to the binomial distribution. The specific setting of the BiGCN single-cell interpolation model includes the following steps:

[0023] Step S-1: In the obtained double cluster expression matrix, the expression value of the i-th cell in the k-th cell subset is The connection between two cells is obtained through the cell connection graph. The similarity calculation formula between two cells is shown in formula (2). The three cases of k values ​​of 5, 10 and 15 are compared. After experimental comparison, k=5 is set:

[0024]

[0025] In the formula, i and j represent cells in the cell subpopulation, and represents the expression level of the gene in cells i and j, k represents the current cell subpopulation, and represent and A vector of , sim represents the similarity between two cells in the same subpopulation,

[0026] calculate Afterwards, a single cell in the single-cell expression matrix will form a cell connection graph together with the five nearest cells. There are two-way arrows between the target node and the first-order aggregated neighbors, and there are also two-way arrows between the first-order aggregated neighbors and the second-order aggregated neighbors, so that the target node can learn the neighbor node information.

[0027] Step S-2: Use GCN to extract features

[0028] After the expression matrix is ​​converted into a cell connection graph, it can enter the graph neural network and update the network nodes of the model through learning; each cell node will change according to the changes of neighboring nodes. The closer it is to the central cell node, the greater the impact on it.

[0029] The first-order graph convolution update is divided into three stages: message transmission stage, message reception stage and message transformation stage.

[0030] The single cell node is R1, and its first-order neighbor node is R2. The transmission function from the central cell to the first-order neighbor is shown in formulas (3) and (4):

[0031]

[0032] In the formula, F k (1) represents the first-order aggregation in subgroup k, and Represent the message sequences received by the target node and the first-order neighbor node, H k (1) represents the emission function of the first-order neighbors in subgroup k,

[0033]

[0034] In the formula, H 1 Represents the output of the first layer, H 0 represents the output result of the input layer, f represents the emission function, represents the weight matrix at the first layer, D represents the degree matrix of the cell node, and represents the different degrees of influence of the first-order neighbors on the core node. represents the weight matrix at the first layer, is the output of a single target node, σ is the ReLU activation function, which aggregates the features of surrounding cells and is used to update the nodes in the next layer. The second layer is the dropout layer. The purpose of setting up this layer is to reduce the occurrence of overfitting, as shown in formula (5):

[0035] H 2 = p(H 1 ) (5)

[0036] In the formula, H 2 Represents the output of the second layer, H 1 represents the output result of the first layer, and p represents the node missing rate. In the cell connection graph, some cell node values ​​will be randomly set to zero.

[0037] The third layer is also a graph convolution layer, which is used to deepen the connection between cells. The formula of this layer is shown in formula (6):

[0038]

[0039] In the formula, H 2 represents the output of the second layer, A represents the adjacency matrix, f represents the emission function, and D represents the degree matrix of the cell node. represents the weight matrix in the second layer, H 3 represents the third layer function, is the output of the target node in the first two layers, σ is the ReLU activation function, which aggregates the features of the surrounding cells and is used to update the nodes in the next layer;

[0040] Step S-3, interpolation using zero-inflated negative binomial distribution

[0041] After the encoder extracts the cell features, it enters the decoder part. The decoder used is different from the conventional decoder. The output layer is three linear layers, and the values ​​expressed are W μ ,W θ and W π , respectively represent the expectation, variance and missing probability of the relevant cells in each subpopulation. The calculation formulas are shown in formula (7), formula (8) and formula (9):

[0042]

[0043] Π=sigmoid(DW π ) (8)

[0044] Θ=exp(DW θ ) (9)

[0045] Where D represents the degree matrix of the cell node, represents the mean of the output in the cell subpopulation, Π represents the variance of the output in the cell subpopulation, Θ represents the sparse probability of the subpopulation output, exp represents the exponential function, sigmoid represents the activation function, and sigmoid is performed on D because the sparsity level is kept between 0 and 1. The loss function used by the model is shown in formula (10) and formula (11):

[0046]

[0047] Loss ZINB(x；π,μ,θ) =min(πδ0(x)+(1-π)NB(x;μ,θ)) (11)

[0048] Where x represents the gene expression value, μ represents the expectation, θ represents the variance, π represents the missing probability, θ / (θ+μ) represents the probability of successfully fitting the zero-value distribution, μ / (θ+μ) represents the probability of fitting failure, Γ(x+θ) represents the factorial of (x+θ-1), Γ(θ) represents the factorial of (θ-1), NB represents the negative binomial distribution, ZINB represents the zero-inflated negative binomial distribution, and min represents minimizing the function. The training parameters are updated iteratively until the algorithm converges.

[0049] In step 4, the specific method of using evaluation indicators for specific evaluation is:

[0050] The evaluation indicators used include Pearson correlation coefficient and mean square error. The performance of the single-cell interpolation model is measured according to the two indicators.

[0051] 1) Pearson correlation coefficient

[0052] The Pearson Correlation Coefficient (PCC) is used to measure the correlation between the true value of the sample and the interpolated value, and is defined as shown in Formula 12:

[0053]

[0054] In the formula, X(gene) and X method (gene) represents the true expression level of a single cell gene in a cell subpopulation and the expression level of a single cell gene after interpolation, Cov is the covariance matrix, Var is the variance, and Corr reflects the similarity between the true gene expression level and the gene expression level after interpolation. The closer its value is to 1, the greater the similarity is and the better the interpolation effect is.

[0055] 2) Mean square error

[0056] The mean square error (MSE) is used to measure the difference between the interpolated value and the true value, and can evaluate the quality of the single-cell interpolation results. It is defined as shown in Formula 13:

[0057]

[0058] In the formula, x i represents the interpolated gene expression value in the cell, y iIt represents the true gene expression value, N represents the total number of samples in the cell, and the smaller the MSE value, the better the interpolation model effect.

[0059] The beneficial effects of the present invention are:

[0060] The present invention solves the problem of insufficient learning of the nonlinear relationship between cells and genes during the interpolation process, and proposes a single-cell transcriptome interpolation method BiGCN based on biclustering and graph convolution. First, similar cells and similar genes are selected using biclustering to form a new biclustering matrix. Secondly, the cell connection graph is constructed using the matrix and input into BiGCN. The neighbor node features of similar cells are extracted, and finally the results are interpolated into the original matrix using the ZINB loss function. The model designed by the present invention compares five algorithms and four evaluation indicators. The results show that it is better than SAVER, MAGIC, scImpute, DCA and DeepImpute, indicating that the use of graph convolution can effectively extract the nonlinear relationship between cells and genes and extract deeper features. In the comparison experiment of different mask rates, the present invention is more stable and has less impact in a high noise environment.

[0061] In view of the problem of insufficient learning caused by nonlinear relationships between cells during the interpolation process, the present invention provides a single-cell transcriptome interpolation method BiGCN based on biclustering and graph convolutional neural network. When constructing a cell connection graph, this method uses biclustering instead of the original KMeans clustering, so that the newly constructed cell connection graph can aggregate more similar cells and genes, and uses a graph convolution layer instead of the original linear layer in the autoencoder to extend the cell relationship from Euclidean space to non-Euclidean space, solving the problem of insufficient learning of similar cells in neural networks. The Pearson correlation coefficient in the mouse data set reached 0.92, which is better than the interpolation method compared by the present invention. In addition, after adding masks in four proportions of 20%, 40%, 50% and 60% in the original sample, the interpolation of the model designed by the present invention is more stable, with fluctuations of no more than 6%.

[0062] The present invention adopts a single-cell interpolation model based on the fusion of graph convolution and zero-inflated negative binomial distribution, which distinguishes the nonlinear relationship between different cells according to the diffusion relationship between the nodes of the cell graph, thereby improving the feature extraction ability. Based on the proposed interpolation model, a cell screening strategy based on double clustering is adopted to aggregate the dispersed similar cells before constructing the cell connection graph, thereby avoiding the interference of redundant zero expression, effectively alleviating data sparsity, and further improving the interpolation accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 The figure is an overall flow chart of the BiGCN single-cell interpolation model of the method of the present invention.

[0064] Figure 2Flowchart of iterative optimization of biclustering for the present invention.

[0065] Figure 3 Schematic diagram of the framework of the BiGCN single-cell interpolation model of the present invention.

[0066] Figure 4 Schematic diagram of the encoder network architecture of BiGCN.

[0067] Figure 5 Diagram of the cell network constructed for the present invention.

[0068] Figure 6 This is a schematic diagram of first-order neighbor sampling of the present invention.

[0069] Figure 7 The data flow diagram is simulated for an embodiment of the present invention.

[0070] Figure 8 A bar chart comparing the MSE error and PCC value of BiGCN according to an embodiment of the present invention and the traditional interpolation method.

[0071] Fig. 9 This is a bar chart comparing the MSE error and PCC value of BiGCN and the deep learning interpolation method according to an embodiment of the present invention.

[0072] Fig.10 It is a histogram of the MSE error and MAE error after interpolation in the simulation data of the embodiment of the present invention.

[0073] Fig.11 It is a histogram of the interpolated PCC value and CS value in the simulation data of the embodiment of the present invention.

[0074] Fig.12 This is a graph of MSE error and PCC values ​​at different mask rates in GSE67602 according to an embodiment of the present invention.

[0075] Fig.13 This is a visualization result diagram of the matrix after DeepImpute interpolation according to an embodiment of the present invention.

[0076] Fig.14 This is a visualization result diagram of the matrix after BiGCN interpolation in an embodiment of the present invention. DETAILED DESCRIPTION

[0077] The present invention is further described in detail below with reference to the accompanying drawings and embodiments.

[0078] The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network includes the following steps:

[0079] See also Figure 1Step 1: Download the gene expression data of human and mouse in matrix form from 10X (https: / / www.10xgenomics.com / ) and GEO (Gene Expression Omnibus) gene databases, respectively, and filter out cells with a total gene expression of less than 3 and genes with a total cell expression of less than 200;

[0080] Step 2: Select features of the gene expression matrix and use biclustering to find cells with consistent growth trends to form a new two-dimensional matrix. Use cells as nodes and cell correlations as features to build a cell connection map.

[0081] Step 3: Use graph convolutional neural network to build BiGCN (full name Bi-clustering Graph Convolution Network) single-cell interpolation model to learn the nonlinear relationship between cells. Divide the data set into training set and test set in a ratio of 8 to 2, use the training set to train the model, and retain the trained BiGCN single-cell interpolation model;

[0082] Step 4: Input the test set into the trained BiGCN single-cell interpolation model, and evaluate the interpolation results using evaluation indicators.

[0083] The step 2, using double clustering to find cells with consistent gene expression level change trends, specifically includes the following steps:

[0084] Step 2-1, preprocessing stage

[0085] Divide the data matrix into single-row and single-column sub-biclusters, which together represent all rows and columns of the expression matrix, and each bicluster contains a row set and a column set, from which some columns and rows are randomly selected to form the index matrix of the starting bicluster;

[0086] Step 2-2, find associated cells

[0087] The constructed index matrix can reflect the changing trend of the expression value. The binary joint information entropy is used to determine the correlation between two cells. Formula (1) is selected

[0088]

[0089] In the formula, X represents the gene expression intensity, m represents the total number of cells, n represents the total number of genes, log represents the logarithmic function with e as the base, and x i represents the i-th cell in a row, y j represents the jth gene in a column, p(x i ,yj ) represents the probability density between cells and genes, H(X) represents the degree of association between the cell and the gene, and the larger the value, the higher the degree of association;

[0090] Step 2-3, iterative optimization

[0091] See also Figure 2 , the binary joint information entropy of the double clusters is optimized in an iterative manner. In each iteration, a row is randomly added or deleted to improve the binary joint information entropy of the cell. The iterative process of seeking the maximum binary joint information entropy uses greedy search, trying to maximize the binary joint information entropy at each step. When the value no longer changes, the iteration is stopped.

[0092] The step 3, constructing the BiGCN single-cell interpolation model, specifically includes the following steps:

[0093] See also Figure 3 BiGCN improves on the traditional autoencoder. In order to further capture the relationship between cells, the graph convolution layer is used instead of the linear layer to learn the characteristics between cells to the greatest extent. Before entering the graph convolutional neural network, the expression matrix is ​​converted into a cell connection graph, and then the ZINB (full name Zero-inflated NegativeBinomial) loss function is used to compare the matrices before and after interpolation to reduce the error.

[0094] See also Figure 4 The autoencoder has a four-layer structure. The first and third layers are graph convolution layers, and the second and fourth layers are Dropout layers. In order to avoid overfitting problems in the autoencoder, the number of convolution layers used cannot be too deep, otherwise it will cause data smoothing; the input layer and output layer of the autoencoder are both 64 layers; the decoder part has three fully connected layers, which respectively represent the mean, variance and sparsity of the model, representing the sparsity of zero values ​​in the current sample, and using the negative binomial distribution result as the loss function.

[0095] The BiGCN single-cell interpolation model described above has a large sparsity of single-cell sequencing data and cannot be fitted using a simple distribution. Since the zero value distribution in the data is not uniform, it is not suitable to use the Poisson distribution. The negative binomial distribution is currently the closest one. The present invention uses a zero-inflated negative binomial distribution based on a modified version of the negative binomial distribution. The results of using it are improved to a certain extent compared to the binomial distribution. The specific setting of the BiGCN single-cell interpolation model includes the following steps:

[0096] Step S-1: In the obtained double cluster expression matrix, the expression value of the i-th cell in the k-th cell subset is The connection between two cells is obtained through the cell connection graph. The similarity calculation formula between two cells is shown in formula (2). The three cases of k values ​​of 5, 10 and 15 are compared. After experimental comparison, k=5 is set:

[0097]

[0098] In the formula, i and j represent cells in the cell subsets, respectively, and k represents the current cell subset. and represent and A vector of , sim represents the similarity between two cells in the same subpopulation,

[0099] calculate After that, a single cell in the single-cell expression matrix will form a cell connection map together with the five nearest cells, see Figure 5 As shown, Figure 5 There are two-way arrows between the target node and the first-order aggregated neighbors, and there are two-way arrows between the first-order aggregated neighbors and the second-order aggregated neighbors, so that the target node can learn the neighbor node information;

[0100] Step S-2: Use GCN (Graph Convolution Network) to extract features

[0101] After the expression matrix is ​​converted into a cell connection graph, it can enter the graph neural network and update the network nodes of the model through learning; Figure 5 Each cell node in the network will change according to the changes of neighboring nodes. The closer it is to the central cell node, the greater the impact on it. The learning process of each layer is as follows: Figure 6 As shown, the five-pointed star part represents the target node, the circle part represents the first-order neighbor, and the pentagon part represents the second-order neighbor.

[0102] The first-order graph convolution update is divided into three stages: message transmission stage, message reception stage and message transformation stage.

[0103] The single cell node is R1, and its first-order neighbor node is R2. The emission function from the central cell to the first-order neighbor is shown in formulas (3) and (4):

[0104]

[0105] In the formula, F k (1) represents the first-order aggregation in subgroup k, and Represent the message sequences received by the target node and the first-order neighbor node, H k (1)represents the emission function of the first-order neighbors in subgroup k,

[0106]

[0107] In the formula, H 1 Represents the output of the first layer, H 0 represents the output result of the input layer, f represents the emission function, represents the weight matrix at the first layer, D represents the degree matrix of the cell node, and represents the different degrees of influence of the first-order neighbors on the core node. represents the weight matrix at the first layer, is the output of a single target node, σ is the ReLU activation function, which aggregates the features of surrounding cells and is used to update the nodes in the next layer. The second layer is the dropout layer. The purpose of setting up this layer is to reduce the occurrence of overfitting, as shown in formula (5):

[0108] H 2 = p(H 1 ) (5)

[0109] In the formula, H 2 Represents the output of the second layer, H 1 represents the output result of the first layer, and p represents the node missing rate. In the cell connection graph, some cell node values ​​will be randomly set to zero.

[0110] The third layer is also a graph convolution layer, which is used to deepen the connection between cells. The formula of this layer is shown in formula (6):

[0111]

[0112] In the formula, H 2 represents the output of the second layer, A represents the adjacency matrix, f represents the emission function, and D represents the degree matrix of the cell node. represents the weight matrix in the second layer, H 3 represents the third layer function, is the output of the target node in the first two layers, σ is the ReLU activation function, which aggregates the features of the surrounding cells and is used to update the nodes in the next layer;

[0113] Step S-3, interpolation using zero-inflated negative binomial distribution

[0114] After the encoder extracts the cell features, it enters the decoder part. The decoder used in this embodiment is different from the conventional decoder. The output layer is three linear layers, and the values ​​expressed are W μ ,W θ and W π, respectively represent the expectation, variance and missing probability of the relevant cells in each subpopulation. The calculation formulas are shown in formula (7), formula (8) and formula (9):

[0115]

[0116] Π=sigmoid(DW π ) (8)

[0117] Θ=exp(DW θ ) (9)

[0118] Where D represents the degree matrix of the cell node, represents the mean of the output in the cell subpopulation, Π represents the variance of the output in the cell subpopulation, Θ represents the sparse probability of the subpopulation output, exp represents the exponential function, sigmoid represents the activation function, and sigmoid is performed on D because the sparsity level is kept between 0 and 1. The loss function used by the model is shown in formula (10) and formula (11):

[0119]

[0120] Loss ZINB(x；π,μ,θ) =min(πδ0(x)+(1-π)NB(x;μ,θ)) (11)

[0121] Where x represents the gene expression value, μ represents the expectation, θ represents the variance, π represents the missing probability, θ / (θ+μ) represents the probability of successfully fitting the zero-value distribution, μ / (θ+μ) represents the probability of fitting failure, Γ(x+θ) represents the factorial of (x+θ-1), Γ(θ) represents the factorial of (θ-1), NB represents the negative binomial distribution, ZINB represents the zero-inflated negative binomial distribution, and min represents minimizing the function. The training parameters are updated iteratively until the algorithm converges.

[0122] In step 4, the specific method of using evaluation indicators for specific evaluation is:

[0123] The evaluation metrics used included Pearson correlation coefficient and mean square error, and the performance of the single-cell imputation model was measured according to the two metrics.

[0124] 1) Pearson correlation coefficient

[0125] The Pearson Correlation Coefficient (PCC) is used to measure the correlation between the true value of the sample and the interpolated value, and is defined as shown in Formula 12:

[0126]

[0127] In the formula, X(gene) and X method (gene) represents the true expression level of a single cell gene in a cell subpopulation and the expression level of a single cell gene after interpolation, Cov is the covariance matrix, Var is the variance, and Corr reflects the similarity between the true gene expression level and the gene expression level after interpolation. The closer its value is to 1, the greater the similarity is and the better the interpolation effect is.

[0128] 2) Mean square error

[0129] The mean square error (MSE) is used to measure the difference between the interpolated value and the true value, and can evaluate the quality of the single-cell interpolation results. It is defined as shown in Formula 13:

[0130]

[0131] In the formula, x i represents the interpolated gene expression value in the cell, y i It represents the true gene expression value, N represents the total number of samples in the cell, and the smaller the MSE value, the better the interpolation model effect.

[0132] The experiment of this embodiment is implemented based on the Windows operating system, and the detailed information is shown in Table 1.

[0133] Table 1 Experimental environment settings

[0134]

[0135] BiGCN interpolation model training parameters

[0136] The training parameters of the BiGCN interpolation model are shown in Table 2. The number of rounds for each training is 500. The training results tend to be flat after 300 rounds.

[0137] Table 2 Specific values ​​or methods of hyperparameters

[0138]

[0139]

[0140] Model data source

[0141] Simulation data plays a very important role in the development of single-cell sequencing technology. Reasonable simulation data will directly affect the subsequent algorithm process. On the one hand, single-cell sequencing simulation data can help researchers develop and optimize algorithms for single-cell data analysis. By simulating the generated data, the performance of different algorithms in different situations can be evaluated, thereby improving the accuracy and efficiency of the results. On the other hand, by generating known simulation data, the accuracy and reliability of the data processing process can be verified. See the simulation data process for details. Figure 7 shown.

[0142] At present, single-cell sequencing data can be used to derive simulated data based on a variety of models. There are eight common types: Poisson distribution, zero-inflated Poisson distribution, mixed distribution of zero-inflated distribution and log-normal Poisson distribution, Dirichlet-multinomial distribution, gamma-multinomial distribution, negative binomial-binomial distribution, gamma-Poisson / negative binomial distribution, and zero-inflated negative binomial distribution. Among them, gamma-Poisson is the most widely used, and many toolkits are also based on this method.

[0143] There are some relatively mature data simulation packages, such as Splatter, SCONE, SymSim, and SinQC, etc. Use the Splatter package in R language to generate simulation data.

[0144] Users can generate corresponding single-cell RNA sequencing data as needed, including different types of cells and different biological conditions. They can specify parameters such as the number of cells, the number of genes, the distribution of cell types, and the gene expression pattern to generate simulated data. Splatter can also introduce simulated noise based on the noise in real single-cell sequencing data, adding labels before data loss, which helps to more realistically reflect the actual situation. The gamma-Poisson model formulas used are shown in formulas (14) and (15).

[0145] λ i,j ~Gamma(α i ,β i ) (14)

[0146] Y i,j ~Poisson(α j λ j,j ) (15)

[0147] In the formula, α i and β i represents the shape parameter and inverse scale parameter of the i-th row gamma distribution, λ i,j Y represents the expression value of the jth gene in the i-th cell. i,j Represents a set of real data, using Poisson distribution to select data. Gamma is the gamma function, and Poisson is the Poisson function.

[0148] Comparison with traditional single-cell interpolation methods

[0149] This example compares BiGCN with three other model-based single-cell interpolation algorithms, MAGIC, SAVER, and scImpute, on four real data sets, with the aim of comparing the differences between traditional algorithms and deep learning algorithms. Figure 8 Table 3 shows the performance comparison.

[0150] Compared with the traditional single-cell interpolation model, the model designed by the present invention has the best effect and is better than other algorithms. Although it is not as good as MixNet in Chapter 3, the PCC index is higher than MixNet. The second is scImpute. Although its PCC effect is good, the interpolation accuracy is not very good, which is 16% lower than the model designed by the present invention. The reason may be that it over-distinguishes true zeros and false zeros on a small data set, resulting in multiple cells overlapping, making the interpolation value random. Next is MAGIC. Due to the use of the Markov affinity matrix, the similarity information between cells is taken into account, but due to the particularity of biological information, it is not completely close to the expression distribution of single cells. It also loses some performance in interpolation accuracy. SAVER has the worst effect, only reaching 0.46 on the mouse data set, which is 30% lower than BiGCN. Because it considers the correlation between genes, but not the relationship between cells and genes.

[0151] Table 3 Comparison of MSE and PCC between BiGCN and different traditional single-cell interpolation methods

[0152]

[0153]

[0154] Note: Bold represents the best effect

[0155] In summary, model-based algorithms have their common disadvantages. Certain false positive values ​​will be generated during the interpolation process. Although the accuracy appears to be improved, the degree of fitting is not high.

[0156] Comparison with Deep Learning Single Cell Interpolation

[0157] BiGCN is compared with three other deep learning single-cell interpolation algorithms DCA, DeepImpute and GNNImpute on four real data sets, in order to compare the difference between the model designed by the present invention and the currently popular deep learning algorithms and reveal the advantages and disadvantages of the model designed by the present invention. Fig. 9 A performance comparison is shown.

[0158] Table 4 Comparison of MSE and PCC between BiGCN and different deep learning single-cell interpolation methods

[0159]

[0160]

[0161] Note: Bold indicates the best result

[0162] On the human data set, the MSE of DeepImpute achieved the best result, followed by the model designed by the present invention, and DCA was the worst. DeepImpute is 3% higher than the model designed by the present invention. The reason is that there are more effective genes in human samples, and the model designed by the present invention will first distinguish true zeros and false zeros based on the model, resulting in more repeated genes, which affects the accuracy. In terms of the PCC indicator, the model designed by the present invention achieved the best result, 4% higher than DeepImpute. The reason is that the model designed by the present invention uses a graph convolution layer, while DeepImpute uses a linear layer, and its ability to extract cell features is slightly worse than that of the model designed by the present invention. DCA uses a combination of mathematical models and autoencoders, taking into account the differences between similar cells, but does not take into account the influence of aggregated neighbor nodes, so the effect is not optimal. Nonlinear relationships cannot be effectively extracted.

[0163] Interpolation results under simulated data set

[0164] The interpolation performance of the model is evaluated using simulated data, because the real data can be known in advance. In the comparative experiment, two evaluation indicators are added to measure the interpolation accuracy, namely the mean absolute error and cosine similarity, as shown in formula (16) and formula (17). The dataset used is generated by the Splatter toolkit, and the interpolation results are shown in Fig.10 , Fig.11 And as shown in Table 5.

[0165] 1) Mean absolute error

[0166] Mean Absolute Error (MAE) is used to measure the average difference between the imputed value and the true value.

[0167]

[0168] Where n is the total number of samples. i represents the true value, Represents the predicted value

[0169] 2) Cosine similarity

[0170] Cosine Similarity (CS), used to evaluate the similarity between two cells

[0171]

[0172] In the formula, i represents the current i-th cell, n represents the total number of cells, A represents the gene expression value of a single cell after interpolation, and B represents the true gene expression value of a single cell. The closer the cosine similarity is to 1, the greater the degree of recovery of the cell.

[0173] According to the experimental results, it can be seen that the single-cell interpolation model proposed by the present invention is the best in most indicators. It can be seen from the bar graph that the results of DCA are closest to the model designed by the present invention, and the MSE is higher. The reason is that DCA adopts negative binomial distribution, which is closer to the simulated data set, while the MSE of SAVER and scImpute are very large. The reason is excessive interpolation, adding a lot of useless zeros, resulting in distorted results. The PCC values ​​of the models designed by the present invention are all the best, because the models designed by the present invention use semi-supervised learning, with some real information, and the use of graph convolution models can effectively learn cell-gene features. The experiment verifies that in a relatively sparse Poisson-gamma distribution environment, the model designed by the present invention still has a good recovery effect.

[0174] Table 5 Comparison of simulation data interpolation results

[0175]

[0176] Note: Bold indicates the best result

[0177] Robustness comparison under different masking rates

[0178] This example tests the interpolation capabilities of different algorithms on different datasets when facing different degrees of masking. Taking the mouse dataset GSE67602 as an example, the original dataset is masked at a ratio of 20%, 40%, 50% and 60%. The results are compared. Fig.12 shown.

[0179] After restoring the original values ​​through the model, the data before and after masking are recorded, the data before masking is used as label data, and the data after masking is restored as prediction data. Three comparative data sets are used to compare MSE and PCC to judge the interpolation effect of the model designed by the present invention. Under different masking rates, the model designed by the present invention is the least affected one in several scenarios. When the masking rate reaches 60%, the PCC can still reach 0.92, and the MSE is 10% lower than that of scImpute. DCA is also almost unaffected, but the expression effect is second only to the model designed by the present invention. Other algorithms are obviously affected by masks to varying degrees, especially scImpute. In scenarios where the masking rate is higher than 40%, the performance loss is very serious.

[0180] Other models will drop significantly when the mask changes, which is obviously affected by the proportion of zero values ​​in the initial matrix. This shows that compared with other models, the model designed by the present invention has certain stability under different masking rates, and has a better recovery effect when the initial value is small.

[0181] Visual analysis

[0182] In this example, the mouse dataset GSE67602 was selected for visualization experiments. The Umap function under the Scanpy package was used for dimensionality reduction. The specific resolutions were set to 0.25, 0.5, and 1, and the comparison algorithms were BiGCN and DeepImpute. The visualization results after clustering are shown in Figure 2. Fig.13 and 14 As shown:

[0183] Figure 13-14 A color represents a cell subpopulation. From the visualization results, it can be seen that after the original matrix is ​​aggregated, the general outline can be distinguished, but the specific boundaries are not clear. When the resolution is 0.25 and 0.5, the aggregation effect of the model designed by the present invention is the best. Compared with the original matrix, this algorithm can divide obvious boundaries, but the aggregation effect of DeepImpute in small subpopulations is more obvious than that of the other two algorithms.

Claims

1. A single-cell sequencing data interpolation method based on biclustering and graph convolutional neural network, characterized in that: The following steps are involved: Step 1, downloading the gene expression data of humans and mice in matrix form from the gene database, and filtering out cells with a total gene expression of less than 3 and genes with a total cell expression of less than 200; Step 2: Select features of the gene expression matrix and use biclustering to find cells with consistent growth trends to form a new two-dimensional matrix. Use cells as nodes and cell correlations as features to build a cell connection map. Step 3: Use graph convolutional neural network to build BiGCN single-cell interpolation model, learn nonlinear relationships between cells, divide the data set into training set and test set, use the training set to train the model, and retain the trained BiGCN single-cell interpolation model; Step 4: Input the test set into the trained BiGCN single-cell interpolation model, and evaluate the interpolation results using evaluation indicators.

2. The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network according to claim 1, characterized in that: The step 2, using double clustering to find cells with consistent gene expression level change trends, specifically includes the following steps: Step 2-1, preprocessing stage Divide the data matrix into single-row and single-column sub-biclusters, which together represent all rows and columns of the expression matrix, and each bicluster contains a row set and a column set, from which some columns and rows are randomly selected to form the index matrix of the starting bicluster; Step 2-2, find associated cells The constructed index matrix can reflect the changing trend of the expression value. The binary joint information entropy is used to determine the correlation between two cells. Formula (1) is selected In the formula, X represents the gene expression intensity, m represents the total number of cells, n represents the total number of genes, log represents the logarithmic function with e as the base, and x i represents the i-th cell in a row, y j represents the jth gene in a column, p(x i ,y j ) represents the probability density between cells and genes, H(X) represents the degree of association between the cell and the gene, and the larger the value, the higher the degree of association; Step 2-3, iterative optimization The binary joint information entropy of the biclustering is optimized in an iterative manner. In each iteration, a row is randomly added or deleted to improve the binary joint information entropy of the cell. The iterative process of seeking the maximum binary joint information entropy uses greedy search, trying to maximize the binary joint information entropy at each step. When the binary joint information entropy no longer changes, the iteration is stopped.

3. The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network according to claim 1, characterized in that: The step 3, constructing the BiGCN single-cell interpolation model, specifically includes the following steps: The BiGCN single-cell interpolation model is improved on the basis of the traditional autoencoder. In order to further capture the relationship between cells, the graph convolution layer is used instead of the linear layer to maximize the learning of the characteristics between cells. Before entering the graph convolutional neural network, the expression matrix is ​​converted into a cell connection graph, and then the ZINB loss function is used to compare the matrices before and after interpolation to reduce the error.

4. The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network according to claim 1, characterized in that: The autoencoder has a four-layer structure, the first and third layers are graph convolution layers, the second and fourth layers are Dropout layers. In order to avoid overfitting problems in the autoencoder, the number of convolution layers used cannot be too deep, otherwise it will cause data smoothing; the input layer and output layer of the autoencoder are both 64 layers; the decoder part has three fully connected layers, which respectively represent the mean, variance and sparsity of the model, representing the sparsity of zero values ​​in the current sample, and using the negative binomial distribution result as the loss function.

5. The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network according to claim 1, characterized in that: The specific setting of the BiGCN single-cell interpolation model includes the following steps: Step S-1: In the obtained double cluster expression matrix, the expression value of the i-th cell in the k-th cell subset is The connection between two cells is obtained through the cell connection graph. The similarity calculation formula between two cells is shown in formula (2). The three cases of k values ​​of 5, 10 and 15 are compared. After experimental comparison, k=5 is set: In the formula, i and j represent cells in the cell subpopulation, and represents the expression level of the gene in cells i and j, k represents the current cell subpopulation, and represent and A vector of , Sim represents the similarity between two cells in the same subpopulation, calculate Afterwards, a single cell in the single-cell expression matrix will form a cell connection graph together with the five nearest cells. There are bidirectional arrows between the target node and the first-order aggregated neighbors, and bidirectional arrows between the first-order aggregated neighbors and the second-order aggregated neighbors, so that the target node can learn the neighbor node information. Step S-2: After using GCN to realize feature extraction and expression matrix conversion into a cell connection graph, it can enter the graph neural network and update the network nodes of the model through learning; The first-order graph convolution update is divided into three stages: message transmission stage, message reception stage and message transformation stage. The single cell node is R1, and its first-order neighbor node is R2. The emission function from the central cell to the first-order neighbor is shown in formulas (3) and (4): In the formula, F k (1) represents the first-order aggregation in subgroup k, and Represent the message sequences received by the target node and the first-order neighbor node, H k (1) represents the emission function of the first-order neighbors in subgroup k, In the formula, H 1 Represents the output of the first layer, H 0 represents the output result of the input layer, A represents the adjacency matrix, f represents the emission function, represents the weight matrix at the first layer, D represents the degree matrix of the cell node, and represents the different degrees of influence of the first-order neighbors on the core node. is the output of a single target node, σ is the ReLU activation function, which aggregates the features of surrounding cells and is used to update the nodes in the next layer. The second layer is the Dropout layer. The purpose of setting up this layer is to reduce the occurrence of overfitting, as shown in formula (5): H 2 =p(H 1 ) (5) In the formula, H 2 Represents the output of the second layer, H 1 represents the output result of the first layer, and p represents the node missing rate. In the cell connection graph, some cell node values ​​will be randomly set to zero. The third layer is also a graph convolution layer, which is used to deepen the connection between cells. The formula of this layer is shown in formula (6): In the formula, H 2 represents the output of the second layer, A represents the adjacency matrix, f represents the emission function, and D represents the degree matrix of the cell node. represents the weight matrix in the second layer, H 3 represents the third layer function, is the output of the target node in the first two layers, σ is the ReLU activation function, which aggregates the features of the surrounding cells and is used to update the nodes in the next layer; Step S-3, interpolation using zero-inflated negative binomial distribution After the encoder extracts the cell features, it enters the decoder part. The decoder used is different from the conventional decoder. The output layer is three linear layers, and the values ​​expressed are W μ , W θ and W π , respectively represent the expectation, variance and missing probability of the relevant cells in each subpopulation. The calculation formulas are shown in formula (7), formula (8) and formula (9): Π=sigmoid(DW π ) (8) Θ=exp(DW θ ) (9) Where D represents the degree matrix of the cell node, represents the mean of the output in the cell subpopulation, Π represents the variance of the output in the cell subpopulation, Θ represents the sparse probability of the subpopulation output, exp represents the exponential function, sigmoid represents the activation function, and sigmoid is performed on D because the sparsity level is kept between 0 and 1. The loss function used by the model is shown in formula (10) and formula (11): Loss ZINB(x；π，μ，θ) =min(πδ0(x)+(1-π)NB(x;μ,θ)) (11) In the formula, x represents the gene expression value, μ represents the expectation, θ represents the variance, π represents the missing probability, θ / (θ+μ) represents the probability of successfully fitting the zero value distribution, and μ / (θ+μ) represents the probability of fitting failure. Γ(x+θ) represents the factorial of (x+θ-1), Γ(θ) represents the factorial of (θ-1), NB represents negative binomial distribution, ZINB represents zero-inflated negative binomial distribution, min represents minimizing the function, and the training parameters are updated iteratively until the algorithm converges.

6. The single-cell sequencing data interpolation method based on double clustering and graph convolutional neural network according to claim 1, characterized in that: The specific method of using evaluation indicators to conduct specific evaluation in step 4 is: The evaluation indicators used include Pearson correlation coefficient and mean square error. The performance of the single-cell interpolation model is measured according to the two indicators. 1) Pearson correlation coefficient The Pearson Correlation Coefficient (PCC) is used to measure the correlation between the true value of the sample and the interpolated value, and is defined as shown in Formula 12: In the formula, X(gene) and X method (gene) represents the true expression level of a single cell gene in a cell subpopulation and the expression level of a single cell gene after interpolation, Cov is the covariance matrix, Var is the variance, and Corr reflects the similarity between the true gene expression level and the gene expression level after interpolation. The closer its value is to 1, the greater the similarity is and the better the interpolation effect is. 2) Mean square error The mean square error (MSE) is used to measure the difference between the interpolated value and the true value, and can evaluate the quality of the single-cell interpolation results. The definition is shown in Formula 13: In the formula, x i represents the interpolated gene expression value in the cell, y i It represents the true gene expression value, N represents the total number of samples in the cell, and the smaller the MSE value, the better the interpolation model effect.

Images