Methods, devices and media for enhancing gene expression interactions in single-cell RNA sequencing data

Abstract

The present invention discloses a method, device and medium for enhancing gene expression interaction in single-cell RNA sequencing data, and belongs to the field of data processing technology. The method first performs principal component analysis on the cell-gene expression spectrum matrix to obtain characteristic genes that can indicate cell differences, calculates cell distance based on characteristic genes, obtains cell similarity based on cell distance, further obtains Markov transition probability based on cell similarity, and performs multiple interpolation on the cell-gene expression spectrum matrix. The method of the present invention can eliminate expression noise in the cell-gene expression spectrum matrix and fill in the missing expression, and finally can effectively enhance gene expression interaction, further can improve the clustering effect of cells, and effectively perform cell typing.

Description

[0001] Related patents

[0002] This application is a divisional application of the Chinese invention patent application with application number 2023107267258, application date June 19, 2023, and invention name “Methods, devices and media for enhancing gene expression interactions in scRNA-seq data”. Technical Field

[0003] The present invention belongs to the field of data processing technology, and in particular, relates to a method, device and medium for enhancing gene expression interactions in single-cell RNA sequencing data. Background Art

[0004] Single cell RNA sequencing (scRNA-seq), also known as single cell transcriptome sequencing, is a high-throughput experimental technology that uses RNA sequencing to quantify the gene expression profile of a specific cell population at the single cell level. It has become a popular technology in recent years. For multicellular organisms, there are differences between cells (cell heterogeneity), that is, cell heterogeneity. This cell heterogeneity can be reflected in different genetic backgrounds, different differentiation states, different physical characteristics, different gene mutation spectra and transcriptome and proteome expression profiles.

[0005] For a single cell, due to insufficient sampling of mRNA molecules, not all mRNA molecules can be captured, and due to the relatively shallow sequencing depth, generally only 10% to 50% of transcripts can be detected in each cell, which results in many gene counts in the cell being 0, resulting in the sparsity of the cell gene expression matrix in the scRNA-seq sequencing results. This sparsity of the cell gene expression matrix increases the computational difficulty of subsequent analysis work and may also seriously obscure the relationship between important genes.

[0006] In order to overcome the sparsity of the cell gene expression matrix in single-cell RNA sequencing results, most current methods cluster thousands of cells into a small number of clusters by clustering; or merge genes by other methods (such as principal component analysis [PCA]) to create "metagene". Although these methods deal with sparsity to a certain extent, they lose the resolution of single cells or single genes. Summary of the invention

[0007] In order to solve the above technical problems, the inventors aim to provide a new method for processing the sparsity in single-cell sequencing data to improve the cell clustering effect. To this end, the technical solution adopted by the present invention is as follows:

[0008] The first aspect of the present invention provides a method for enhancing gene expression interactions in scRNA-seq data, comprising the following steps:

[0009] S1, obtain the cell-gene expression profile matrix A of single-cell transcriptome sequencing data;

[0010] S2, based on the cell-gene expression profile matrix A, using a principal component analysis method to screen N principal components, where N=20-50, to obtain a PCA matrix;

[0011] S3, calculating the distance d(p, q) between any two cells based on the values ​​of the N principal components to obtain a cell distance matrix D;

[0012] S4, calculating the similarity C(p, q) between any two cells using a kernel function based on the cell distance matrix D to obtain a cell similarity matrix C;

[0013] S5, based on the cell similarity matrix C, according to the following algorithm, calculate the transition probability M(p, q) between any two cells to obtain the cell transition probability matrix M:

[0014] The two cells include a first cell and a second cell, and the transition probability between the first cell and the second cell is a ratio of the similarity between the first cell and the second cell to the sum of the similarities between the first cell and all cells;

[0015] S6, based on the cell transition probability matrix M, interpolate the cell-gene expression profile matrix A according to the following algorithm:

[0016] The expression level of the ith gene in the pth cell is converted into the sum of the product of the expression levels of the ith genes in the remaining cells and the transition probability between the pth cell and the corresponding cell.

[0017] Wherein, p=1-m and q=1-m, m represents the number of cells in the cell-gene expression spectrum matrix A, i=1-n, n represents the number of genes in the cell-gene expression spectrum matrix A.

[0018] In some embodiments of the present invention, before performing principal component analysis in step S2, a step of normalizing the cell-gene expression profile matrix A according to the library size is further included, and the normalization method is:

[0019]

[0020] Among them, A norm (p, i) represents the expression level of the i-th gene in the p-th cell after normalization, and A(p, i) represents the expression level of the i-th gene in the p-th cell before normalization; represents the sum of the expression levels of all genes in the pth cell, j = 1 to n; mean(LS) represents the average of the sum of the expression levels of all genes in all cells.

[0021] Principal Component Analysis (PCA) is the most widely used data dimensionality reduction algorithm. Dimensionality reduction is a method of preprocessing high-dimensional feature data. Dimensionality reduction is to retain the most important features of high-dimensional data and remove noise and unimportant features, thereby achieving the purpose of improving data processing speed.

[0022] In some embodiments of the present invention, the first N principal components (PCs) that can reflect the differences of the original cells are selected as data for subsequent analysis. When N=20-50, it can reflect more than about 80% of the information. In some preferred embodiments of the present invention, N=30 is selected, that is, 20 principal components are selected for analysis, and the results are sufficiently robust.

[0023] In some embodiments of the present invention, in step S3, the distance is the Euclidean distance, and the calculation formula of d(p, q) is as follows:

[0024]

[0025] Among them, PC i (p) represents the value of the i-th PC in the p-th cell, PC i (q) represents the value of the i-th PC in the q-th cell; i=1~N.

[0026] In some embodiments of the present invention, in step S4, before calculating the similarity, the PCA matrix is ​​further subjected to nonlinear dimensionality reduction using the UMAP (Uniform Manifold Approximation and Projection) algorithm, and the result after dimensionality reduction is two-dimensional coordinate information, and cells are clustered into different regions based on phenotypic similarity. In the present invention, cells in different regions have an insignificant interpolation effect on the sparse matrix, and therefore, for cells in different regions, the similarity between each pair is set to 0, and for cells in the same region, the similarity between any two cells is calculated using a kernel function,

[0027] The kernel function is a Gaussian kernel function, and the calculation formula of C(p, q) is as follows:

[0028]

[0029] Wherein, σ is the bandwidth, which is used to control the radial action range, σ = 1 to 30.

[0030] Different cells have different gene expressions, which reflect the possible phenotypes of cells. Phenotypes can be reflected in cells as cell types, such as the state of developmental period, etc. Since the phenotypes of different cells are not consistent, different σ values ​​are selected for different cell types. If σ is too small (less than 0.01), the results will be unstable and the accuracy will be reduced, that is, cells with the same phenotype will be difficult to classify as the same type of cells; the larger σ is, the larger the local influence range of the Gaussian kernel function will be. If it is too large (greater than 100), overfitting will occur, that is, cells with different phenotypes and far distances will be averaged together, and the data resolution will be lost.

[0031] In some embodiments of the present invention, the value of σ is determined according to the following method:

[0032] First determine the density of the area where cells p and q are located. When the density is less than 0.3, σ is 20 to 30. When the density is greater than or equal to 0.3 and less than or equal to 0.6, σ is 5 to 20. When the density is greater than 0.6, σ is 1 to 5.

[0033] In some embodiments of the present invention, after obtaining the cell transition probability matrix M, the following processing is further performed:

[0034] (1) For any cell in the same region, the 15 cells with the smallest sum of distances to the two cells are determined as the neighbor cells of the two cells. If there are less than 15 cells, all the remaining cells are neighbor cells;

[0035] (2) For each cell, the cell transition probability between it and non-neighboring cells is set to 0.

[0036] In some embodiments of the present invention, in step S6, the cell-gene expression profile matrix A is multiply interpolated by exponentiating the transition probability matrix M:

[0037] A imputed， t=M t ×A

[0038] Among them, A imputed，t represents the cell-gene expression spectrum matrix after the t-th interpolation; t represents the number of iterations, starting from 1 and increasing by 1 for multiple interpolation; when the t-th interpolation is performed, for each gene in each cell, the expression amounts used for conversion are the expression amounts of each gene in the cell-gene expression spectrum matrix A.

[0039] In the present invention, since the eigenvalues ​​of the Markov transition probability matrix M are all between [0, 1], the eigenvalues ​​will gradually decrease through exponential operation, and their range is also between [0, 1]. As the Markov transition probability matrix M is repeatedly exponentiated, the size of all eigenvalues ​​except 1 continues to decrease, which can reduce the importance of noise and its explanatory power is close to zero. As t increases, cells learn missing values ​​from their neighbor cells and quickly obtain relationships between biologically very similar cells.

[0040] In some embodiments of the present invention, the number of iterations is determined in the following manner:

[0041] After the tth interpolation, first calculate A imputed，t The relative expression of the ith gene in the pth cell:

[0042]

[0043] Then calculate the average gene expression after interpolation of the pth cell:

[0044]

[0045] Then, based on the residual sum of squares of the relative expression of all genes after and after the t-th interpolation and the deviation sum of squares of the relative expression of all genes after interpolation, the determination coefficient of cell p was calculated. The formula is as follows:

[0046]

[0047]

[0048]

[0049] If the determination coefficient of all cells is less than the preset threshold P1, the iteration is stopped and the cell-gene expression profile matrix A after multiple interpolation is obtained. imputed ,

[0050] Among them, j=1~n, P1=3%~10%.

[0051] In some embodiments of the present invention, when performing the tth interpolation, if M t (p, q) is less than the preset threshold P2, then M is set first t (p, q) = 0, and then interpolation is performed, where P2 = 0.00005 ~ 0.001.

[0052] A second aspect of the present invention provides a computer device, comprising:

[0053] Memory for storing computer programs;

[0054] A processor is used to implement the steps of a method for enhancing gene expression interactions in scRNA-seq data as described in any one of the first aspects of the present invention when executing the computer program.

[0055] A third aspect of the present invention provides a computer-readable storage medium,

[0056] The computer-readable storage medium stores a computer program, and when the computer program is executed by the processor, the steps of a method for enhancing gene expression interactions in scRNA-seq data as described in any one of the first aspects of the present invention are implemented.

[0057] Beneficial Effects of the Invention

[0058] Compared with the prior art, the technical effects of the present invention are as follows:

[0059] The method of the present invention performs multiple interpolation on single-cell transcriptome sequencing data based on a Markov stochastic model, thereby being able to eliminate expression noise in a cell-gene expression spectrum matrix and fill in missing transcript counts.

[0060] By using the method of the present invention to perform multiple interpolations on the cell-gene expression spectrum matrix obtained from single-cell transcriptome sequencing data, gene expression interactions can be effectively enhanced, and the clustering effect of cells can be further improved, thereby effectively performing cell typing. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 The cell-gene expression profile matrix (partial) obtained by single-cell transcriptome sequencing is shown.

[0062] Figure 2 The cluster diagram is shown after the PCA matrix information is subjected to nonlinear dimensionality reduction again using the UMAP algorithm in Example 1 of the present invention.

[0063] Figure 3 The results of CD3D gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 1 of the present invention are shown.

[0064] Figure 4 The results of MS4A1 gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 1 of the present invention are shown.

[0065] Figure 5 The results of the relationship analysis between the CD14 gene and the ITGAM gene before and after weighted multiple interpolation of the human PBMC single-cell transcriptome sequencing data in Example 1 of the present invention are shown.

[0066] Figure 6The results of CD3D gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 2 of the present invention are shown.

[0067] Figure 7 The results of MS4A1 gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 2 of the present invention are shown.

[0068] Figure 8 The results of the relationship analysis between the CD14 gene and the ITGAM gene before and after weighted multiple interpolation of the human PBMC single-cell transcriptome sequencing data in Example 2 of the present invention are shown.

[0069] Fig. 9 The results of CD3D gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 3 of the present invention are shown.

[0070] Fig.10 The results of MS4A1 gene expression analysis before and after weighted multiple interpolation of human PBMC single-cell transcriptome sequencing data in Example 3 of the present invention are shown.

[0071] Fig.11 The results of the relationship analysis between the CD14 gene and the ITGAM gene before and after weighted multiple interpolation of the human PBMC single-cell transcriptome sequencing data in Example 3 of the present invention are shown. DETAILED DESCRIPTION

[0072] Unless otherwise indicated, implied from the context, or customary in the prior art, all parts and percentages in this application are based on weight, and the tests and characterization methods used are all synchronized with the filing date of this application. Where applicable, the contents of any patent, patent application or disclosure involved in this application are fully incorporated herein by reference, and their equivalent patent families are also introduced as references, especially the definitions of relevant terms in the art disclosed in these documents. If the definition of a specific term disclosed in the prior art is inconsistent with any definition provided in this application, the definition of the term provided in this application shall prevail.

[0073] In order to make the technical problems, technical solutions and beneficial effects solved by the present invention more clearly understood, the present invention is further described in detail below in conjunction with embodiments.

[0074] Example

[0075] The following examples are used to demonstrate preferred embodiments of the present invention. It will be appreciated by those skilled in the art that the techniques disclosed in the following examples represent techniques discovered by the inventors that can be used to implement the present invention and therefore can be considered as preferred embodiments of the present invention. However, it will be appreciated by those skilled in the art based on this specification that many modifications may be made to the specific embodiments disclosed herein and still achieve the same or similar results without departing from the spirit or scope of the present invention.

[0076] Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs, and the disclosure and materials cited therein are hereby incorporated by reference.

[0077] Those skilled in the art will recognize, or be able to ascertain using no more than routine experimentation, many technical equivalents to the specific embodiments of the invention described herein. Such equivalents are intended to be encompassed by the claims.

[0078] The experimental methods in the following examples are conventional methods unless otherwise specified. The instruments and equipment used in the following examples are conventional laboratory instruments and equipment unless otherwise specified; the experimental materials used in the following examples are purchased from conventional biochemical reagent stores unless otherwise specified.

[0079] Example 1: Single-cell sequencing gene matrix sparsity processing

[0080] 1. Obtaining cell-gene expression profile matrix

[0081] Human peripheral blood mononuclear cells (PBMC) were collected and single-cell sequencing was performed to obtain a cell-gene expression matrix A, part of which is shown in Figure 1 As shown in Table 1, Figure 1 In the graph, the columns are different cells (represented by Barcode, which is composed of 16 random nucleotides), and the rows represent the expression levels of genes in different cells (transcript counts). · " represents the expression level is 0, that is, the gene transcript count is 0. It can be seen that the expression levels of many genes in the cell-gene expression matrix A are 0.

[0082] In this example, the cell-gene expression profile includes the expression levels of 36,477 genes in 2,887 cells.

[0083] Table 1 Cell-gene expression profile matrix

[0084]

[0085] 2. Normalization and principal component analysis of cell-gene expression matrix A

[0086] In order to ensure that the distance between cells reflects the real biological difference and is not affected by human factors, the original expression matrix is ​​processed in the following two steps:

[0087] (1) Normalize the cell-gene expression matrix A according to the library size

[0088] The normalization method is:

[0089]

[0090] Among them, A norm (p, i) represents the expression level (transcript count) of the i-th gene in the p-th cell after normalization, p = 1 to m, m represents the number of cells in the cell-gene expression matrix A, and A(p, i) represents the expression level of the i-th gene in the p-th cell before normalization; represents the sum of the expression levels of all genes in the p-th cell (the sum of the transcript counts), i=1~n and j=1~n, n represents the number of genes in the cell-gene expression matrix A; LS represents the cell library size set, that is, the set composed of the sum of the expression levels of all genes in each cell (the sum of the transcript counts), mean(LS) represents the average value of the sum of the transcript counts of all cells.

[0091] (2) Principal component analysis

[0092] Principal Component Analysis (PCA) is the most widely used data dimensionality reduction algorithm. Dimensionality reduction is a method of preprocessing high-dimensional feature data. Dimensionality reduction is to retain the most important features of high-dimensional data and remove noise and unimportant features, thereby achieving the purpose of improving data processing speed.

[0093] By performing principal component analysis on the normalized data, n genes, i.e., n-dimensional data, are reduced to 50 principal components (PCs), and the first 30 principal components (PCs) that can reflect the differences in the original cells are further selected as data for subsequent analysis. The inventors have confirmed that 30 principal components can reflect more than about 80% of the information, so 30 principal components are selected for analysis, and the results are sufficiently robust.

[0094] For Cell 1 and Cell 2, the 30 PC values ​​are shown in Table 2.

[0095] Table 2 PCA matrix of cell 1 and cell 2

[0096]

[0097]

[0098] The top 10 PC values ​​of some cells are shown in Table 3.

[0099] Table 3 PCA matrix of some cell parts

[0100]

[0101]

[0102] 3. Calculate the cell-cell distance matrix D

[0103] The Euclidean distance between cells was calculated using the PC values ​​based on the cells. The calculation formula is as follows:

[0104]

[0105] Where d(p,q) represents the Euclidean distance between the pth cell and the qth cell, p = 1~m and q = 1~m, m represents the number of cells, PC i (p) represents the value of the i-th PC in the p-th cell, PC i (q) represents the value of the i-th PC in the q-th cell; i = 1 to N, N represents the number of principal components used to calculate the Euclidean distance of cells, N = 30. When p = q, that is, for the same cell, the Euclidean distance value is 0. Calculation shows that the Euclidean distance between cell 1 and cell 2 is 22.95, and the Euclidean distance between any two cells of the 2887 cells can be calculated.

[0106] The Euclidean distance matrix D between some cells is shown in Table 4:

[0107] Table 4 Cell Euclidean distance matrix

[0108]

[0109] 4. Use adaptive Gaussian kernel function to convert the distance matrix into similarity matrix C

[0110] The UMAP algorithm is used to perform nonlinear dimensionality reduction on the PCA matrix information again. The result after dimensionality reduction is two-dimensional coordinate information, which can present the specific location of cells on low-dimensional coordinates. Cells are clustered into different areas based on phenotypic similarities, such as Figure 2 For cells in different regions, the similarity between them is set to 0.

[0111] For the correlation between cells clustered in the same cell region, the spearman / person algorithm was generally used in the past to calculate the correlation coefficient between two cells by comparing the normalized gene expression in the cells. However, due to the strong sparsity of single-cell transcriptome sequencing data, that is, many genes cannot be detected with transcripts, the calculated correlations are all lower than 0.3, making it difficult to distinguish the strength of similarity between cells.

[0112] Therefore, in this embodiment, the inventors introduce an adaptive Gaussian kernel function to represent the correlation between two cells clustered in the same area.

[0113] Gaussian kernel function, also known as radial basis function (RBF) function, is a commonly used kernel function that can map finite-dimensional data to high-dimensional space to achieve the purpose of distinguishing two vectors.

[0114] The Gaussian kernel function formula is as follows:

[0115]

[0116] Among them, C(p,q) represents the similarity between cells p and q; σ is the bandwidth, which controls the radial range of action.

[0117] It can be seen that the Gaussian kernel function is a monotonic function of the Euclidean distance between two cells, and the correlation between two cells shows a monotonically decreasing relationship as the Euclidean distance between the two cells increases.

[0118] Different cells have different gene expressions, which reflect the possible phenotypes of cells. Phenotypes can be reflected in cells as cell types, such as the state of developmental period, etc. Since the phenotypes of different cells are not consistent, different σ values ​​are selected for different cell types. If σ is too small (less than 0.01), the results will be unstable and the accuracy will be reduced, that is, cells with the same phenotype will be difficult to classify as the same type of cells; the larger σ is, the larger the local influence range of the Gaussian kernel function will be. If it is too large (greater than 100), overfitting will occur, that is, cells with different phenotypes and far distances will be averaged together, and the data resolution will be lost.

[0119] Depend on Figure 2It can be seen that the density of cells (different cell types) in different regions is different. For example, the number of mature cells is much larger than the number of stem cells, so their corresponding densities are also inconsistent. The range of density is [0, 1]. The neighbors of cells in denser areas will be larger than those in smaller areas. In order to ensure that the number of neighbors is as consistent as possible, smaller σ values ​​are selected for denser areas and larger σ values ​​are selected for smaller areas. When the density is less than 0.3, σ is 25, when the density is greater than or equal to 0.3 and the density is less than or equal to 0.6, σ is 10, and when the density is greater than 0.6, σ is 3. In this way, cells in different local areas on low-dimensional coordinates can be determined to have closer neighbor cells.

[0120] The above cells 1 and 2 are clustered in the same area, the density of which is 0.4, and σ is set to 15. The Gaussian kernel function for obtaining cells 1 and 2 is:

[0121] Thus, the similarity between any two cells can be obtained, and the similarity matrix C can be obtained. The similarity of the same cells is 1. The similarity matrix C between some cells is shown in Table 5:

[0122] Table 5 Similarity matrix C between cells

[0123]

[0124] 5. Calculate the Markov transition probability matrix

[0125] During the single-cell transcriptome sequencing process, since the capture of mRNA is random, different cells have the characteristics of random mRNA loss, and this process can be characterized as a Markov process. The Markov process has the following characteristics: under the condition of knowing the current state (now), its future evolution (future) does not depend on its previous evolution (past), and the transition of each state depends only on the previous state. For single-cell transcriptome sequencing data, the probability of different cells converting to other cells is different.

[0126] Based on the similarity matrix C, the Markov transition probability matrix M can be obtained. The probability that cell p transforms into cell q (Markov transition probability) is the ratio of the similarity between cell p and cell q to the sum of the similarities between cell p and each other cell:

[0127]

[0128] Wherein, M(p, q) represents the Markov transition probability between cell p and cell q, p = 1 ~ m and q = 1 ~ m, M(p, k) represents the Markov transition probability between cell p and cell k, k represents the kth cell (k = 1 ~ m, k ≠ p), and m represents the number of cells; Represents the sum of similarities between cell p and every other cell.

[0129] It can be seen that for each column or row of the Markov transition probability matrix, the sum of its transition probabilities is 1. The Markov transition probability matrix M between some cells is shown in Table 6:

[0130] Table 6 Markov transition probability matrix between cells

[0131] Cell 1 Cell 2 Cell 3 Cell 4 Cell 5 Cell 6 Cell 7 Cell 8 Cell 9 …… Cell 2887 Cell 1 1 3.91E-03 7.92E-03 3.90E-03 1.35E-03 1.60E-03 0 1.62E-03 1.62E-03 …… 1.62E-03 Cell 2 3.91E-03 1 4.60E-03 8.48E-03 1.67E-02 8.96E-03 0 2.47E-02 2.52E-02 …… 4.81E-03 Cell 3 7.92E-03 4.60E-03 1 1.52E-02 4.01E-03 4.23E-03 0 1.07E-03 6.96E-03 …… 7.11E-03 Cell 4 3.9E-03 8.48E-03 1.52E-02 1 9.31E-03 9.31E-03 9.31E-03 9.31E-03 9.31E-03 …… 1.39E-02 Cell 5 1.35E-03 1.67E-02 4.01E-03 9.31E-03 1 3.26E-03 0 1.14E-02 6.18E-02 …… 8.48E-03 Cell 6 1.60E-03 8.96E-03 4.23E-03 1.30E-02 3.26E-3 1 0 9.42E-03 9.42E-03 …… 1.12E-02 Cell 7 0 0 0 0 0 0 1 0 0 …… 0 Cell 8 1.62E-03 2.47E-02 1.07E-03 9.14E-03 1.14E-02 9.42E-03 0 1 7.13E-03 …… 1.52E-03 Cell 9 1.62E-03 2.52E-02 6.96E-03 1.43E-02 6.18E-02 9.42E-03 0 7.13E-03 1 …… 8.11E-03 …… …… …… …… …… …… …… …… …… …… 1 9.45E-03 Cell 2887 1.62E-03 4.81E-03 7.11E-03 1.39E-02 8.48E-03 1.12E-02 0 1.52E-03 8.11E-03 9.45E-03 1

[0132] 6. Weighted multiple imputation of gene-cell expression matrix using Markov transition probability matrix

[0133] For single-cell transcriptome sequencing data, the main basis for distinguishing cell phenotypes is highly variable genes (also called eigenvectors), and the real structural characteristics of the data are mainly reflected by the top eigenvectors, and the remaining eigenvectors may be noise. The eigenvalues ​​in the Markov transition probability matrix also reflect this feature information. Therefore, for the above steps, the top eigenvector is first retained based on the eigenvalue size of the Markov transition probability matrix M. The top eigenvector mainly comes from the low-frequency eigenvalues.

[0134] Since the noise in the data generally presents a higher frequency in the data, the Markov transition probability matrix M is exponentiated, and the eigenvalues ​​(transition probabilities) in the Markov transition probability matrix M can be low-pass filtered, that is, the low frequencies are allowed to pass through and the high frequencies are filtered out or attenuated, thereby filtering out the noise contained in the high frequencies.

[0135] Furthermore, the cell-gene expression profile matrix A is weighted multiple interpolation based on the exponentiation Markov transition probability matrix:

[0136] A imputed，t =M t ×A

[0137] Among them, A imputed，t represents the cell-gene expression spectrum matrix after the t-th interpolation; t represents the number of iterations, which is also a power, starting from 1 and increasing by 1 for multiple interpolation; when the t-th interpolation is performed, for each gene in each cell, the expression amounts used for the conversion are the expression amounts of each gene in the cell-gene expression spectrum matrix A.

[0138] After the tth interpolation, the cell-gene expression matrix A imputed，t The expression level of gene A in the pth cell imputed，t (p, i) will become the sum of the expression levels of the gene in each neighboring cell of the p-th cell in the cell-gene expression matrix A multiplied by the probability (power t) that the p-th cell is converted into each neighboring cell:

[0139]

[0140] Since the eigenvalues ​​of the Markov transition probability matrix M are all between [0, 1], the eigenvalues ​​will gradually decrease through exponential operations, and their range is also between [0, 1]. As the Markov transition probability matrix M is exponentially multiplied, the size of all eigenvalues ​​except 1 continues to decrease, which can reduce the importance of noise and its explanatory power is close to zero. As t increases, cells learn missing values ​​from their neighbor cells and quickly acquire relationships between biologically very similar cells.

[0141] When the noise removal is transformed into the removal of the real biological information signal, t reaches the optimal number of iterations. Since the noise usually has different frequencies from the real signal itself (i.e., high frequency and low frequency, respectively), as the high-frequency information is removed, the data will change rapidly, and then change slowly to reach convergence. When the data change rate is lower than the first lower threshold (set to 5% in this embodiment), it can be considered that convergence is achieved and t reaches the optimal number of iterations. Finally, the cell-gene weighted expression interpolation matrix A is obtained imputed，t , to reduce data noise and efficiently calculate missing transcript expression without overfitting the data.

[0142] The rate of change of cell p data is determined by the coefficient of determination To quantify:

[0143] For the number of iterations t, the relative expression of gene i after interpolation is calculated first (gene expression divided by the sum of all gene expression).

[0144]

[0145] Then calculate the average gene expression after cell p interpolation:

[0146]

[0147] Then, based on the residual sum of squares of the relative expression of all genes before and after the tth interpolation (SSE t (p)) and the sum of squares of the relative expression of all genes after interpolation (SST t (p)) Calculate the coefficient of determination of cell p The formula is as follows:

[0148]

[0149]

[0150]

[0151] In this implementation, when t=3, the coefficient of determination of all cells is less than 5%. Thus, the cell-gene expression profile matrix A after weighted interpolation is obtained: imputed , effectively reducing data noise and filling in missing transcript expression levels without overfitting the data.

[0152] Cell-gene expression profile matrix A before and after multi-weighted interpolation imputed After analysis, it was found that the expression of the same gene in adjacent cells was more consistent and more realistic after weighted interpolation. For example, for the T cell marker gene CD3D and the B cell marker gene MS4A1, the gene expression comparison is as follows: Figure 3 and Figure 4 As shown. It can be seen that the cell-gene expression profile matrix A after multiple weighted interpolation imputed In the multi-weighted interpolation, the gene expression levels are more clearly expressed between adjacent cells, and the characteristics of the cell type can be intuitively judged and the specific cell type can be determined. On the contrary, before multi-weighted interpolation, the gene expression levels are less consistent even among cells in the same cluster, and it is difficult to distinguish different cell populations based on the expression levels of cell type-specific marker genes.

[0153] In addition, as the expression of each gene in the cell increases with the real signal, the sparseness of the gene in each cell decreases, the number of cells with zero expression decreases, and the correlation between genes increases accordingly. For example, CD14 gene and ITGAM gene are both marker genes of myeloid cells, and the expression correlation of these two genes is strong in normal myeloid cells. However, it is difficult to find this correlation using the cell-gene expression matrix A before multiple weighted interpolation. imputed After analysis, we can see that the correlation between the two is obvious. Figure 5 Therefore, the method of the present invention also plays a very positive role in the identification of cell types.

[0154] Example 2 Optimization of single-cell sequencing gene matrix sparsity processing

[0155] This embodiment further optimizes the solution of embodiment 1, specifically:

[0156] When calculating the Markov transition probability matrix, the following processing is performed:

[0157] For any two cells in the same region (cluster, type), the 15 cells with the smallest sum of distances to the two cells are determined as the neighbor cells of the two cells. If there are less than 15 cells, all the remaining cells are neighbor cells.

[0158] For each cell, the Markov transition probability between it and non-neighboring cells is set to 0.

[0159] The Markov transition probability matrix optimized in this embodiment is used to perform weighted multiple interpolation on the gene-cell expression profile matrix for further analysis.

[0160] Similarly, for the T cell marker gene CD3D and the B cell marker gene MS4A1, the gene expression comparison is as follows Figure 6 and Figure 7 As shown. It can be seen that after the optimization method, the gene expression between adjacent cells is more clearly expressed, and the characteristics of the cell type can be intuitively judged and the specific cell type can be determined. For the interaction between cells, such as the CD14 gene and the ITGAM gene, the correlation between the two is further strengthened, such as Figure 8 shown.

[0161] Example 3 Further optimization of single-cell sequencing gene matrix sparsity processing

[0162] This embodiment further optimizes the solution of embodiment 2, specifically:

[0163] In the step of weighted multiple interpolation of the gene-cell expression spectrum matrix using the Markov transition probability matrix, when the Markov transition probability matrix M is exponentiated, when the power is, if a certain eigenvalue is lower than 0.001, the eigenvalue is assigned 0 and then interpolation is performed.

[0164] The optimized method of this embodiment was used to perform weighted multiple interpolation on the gene-cell expression profile matrix and then analyzed. The results showed that the optimized method can remove most of the noise in the data.

[0165] Similarly, for the T cell marker gene CD3D and the B cell marker gene MS4A1, the gene expression comparison is as follows Fig. 9 and Fig.10 As shown. It can be seen that after the optimization method, the gene expression between adjacent cells is more clearly expressed, and the characteristics of the cell type can be intuitively judged and the specific cell type can be determined. For the interaction between cells, such as the CD14 gene and the ITGAM gene, the correlation between the two is further strengthened, such as Fig.11 shown.

[0166] All documents mentioned in the present invention are cited as references in this application, just as each document is cited as reference individually. In addition, it should be understood that after reading the above teachings of the present invention, those skilled in the art can make various changes or modifications to the present invention, and these equivalent forms also fall within the scope defined by the claims attached to this application.

Claims

1. A method to enhance gene expression interactions in single-cell RNA sequencing data, It is characterized in that The following steps are involved: S1, Obtaining the cell-gene expression profile matrix of single-cell transcriptome sequencing data A , and the cell-gene expression profile matrix A Normalize according to library size; S2, based on the cell-gene expression profile matrix A , use the principal component analysis method to screen N principal components, where N=30, and get the PCA matrix; S3, calculating the distance between any two cells based on the values ​​of the N principal components , and obtain the cell distance matrix D ; S4, based on the cell distance matrix D Use kernel function to calculate the similarity between any two cells , and obtain the cell similarity matrix C ; S5, based on the cell similarity matrix C , according to the following algorithm, calculate the transition probability between any two cells , and obtain the cell transition probability matrix : The two cells include a first cell and a second cell, and the transition probability between the first cell and the second cell is the ratio of the similarity between the first cell and the second cell to the sum of the similarities between the first cell and all cells, and the cell transition probability matrix is ​​obtained. After that, further processing is performed as follows: (1) For any cell in the same region, determine the 15 cells with the smallest sum of distances to the two cells as the neighbor cells of the two cells. If there are less than 15 cells, all the remaining cells are neighbor cells; (2) For each cell, the cell transition probability between it and non-neighboring cells is set to 0; S6, based on the cell transfer probability matrix , according to the following algorithm, the cell-gene expression matrix A To do interpolation: The first Cell The expression level of each gene is converted to the The expression level of each gene is related to the The sum of the products of cells and the transition probabilities of the corresponding cells, in, =1~ and =1~ , Representing the cell-gene expression profile matrix A The number of cells in =1~ , Representing the cell-gene expression profile matrix A The number of genes in By transforming the transition probability matrix Exponentiation, the cell-gene expression profile matrix A Perform multiple imputation: in, Representative The cell-gene expression profile matrix after sub-interpolation; Represents the number of iterations, starting from 1 and increasing by 1 for multiple interpolation; when the During the interpolation, for each gene in each cell, the expression values ​​used for conversion are the cell-gene expression profile matrix A The expression level of each gene in The number of iterations is determined as follows: Conduct the After the first interpolation, if Cell If the expression level of a gene is lower than 0.001, the expression level of the gene is assigned a value of 0 and then calculated. Middle Cell The relative expression of genes: Calculate the Average gene expression after cell interpolation: Based on the The residual sum of squares of the relative expression of all genes after and after interpolation and the deviation sum of squares of the relative expression of all genes after interpolation were calculated. The coefficient of determination , the formula is as follows: If the determination coefficient of all cells is less than the preset threshold P1, the iteration is stopped to obtain the cell-gene expression profile matrix after multiple interpolation , in, =1~ , P1=3%~10%.

2. A method for enhancing gene expression interactions in single-cell RNA sequencing data according to claim 1, It is characterized in that In step S3, the distance is the Euclidean distance, The calculation formula is as follows: in, Representative i PC in the The value in the cell, Representative i PC in the The value in each cell; =1~N.

3. A method for enhancing gene expression interactions in single-cell RNA sequencing data according to claim 1 or 2, It is characterized in that In step S4, before calculating the similarity, the PCA matrix is ​​further subjected to nonlinear dimensionality reduction using the UMAP algorithm. The result after dimensionality reduction is two-dimensional coordinate information. Cells are clustered into different regions based on phenotypic similarity. For cells in different regions, the similarity between each other is set to 0. For cells in the same region, the kernel function is used to calculate the similarity between any two cells. The kernel function is a Gaussian kernel function. The calculation formula is as follows in, is the bandwidth, used to control the radial range of action, =1~30.

4. A method for enhancing gene expression interactions in single-cell RNA sequencing data according to claim 3, It is characterized in that Determined by Value: First determine the cells and cells The density of the area. When the density is less than 0.3, The value is 20~30, when the density is greater than or equal to 0.3 and the density is less than or equal to 0.6, The value ranges from 5 to 20. When the density is greater than 0.6, The value range is 1~5.

5. A computer device, It is characterized in that include: Memory for storing computer programs; A processor, configured to implement the steps of a method for enhancing gene expression interactions in single-cell RNA sequencing data as described in any one of claims 1 to 4 when executing the computer program.

6. A computer-readable storage medium, It is characterized in that The computer-readable storage medium stores a computer program, which, when executed by a processor, implements the steps of a method for enhancing gene expression interactions in single-cell RNA sequencing data as described in any one of claims 1 to 4.

Images