A gene expression data clustering method based on deep autoencoder

Abstract

The invention belongs to the field of gene expression data analysis, and specifically relates to a gene expression data clustering method based on a deep autoencoder; the method comprises: acquiring gene expression data; preprocessing the gene expression data to obtain mRNA data; inputting the mRNA data into a trained deep autoencoder to obtain mRNA data after dimensionality reduction; using an mclust package to determine the number of clusters of the mRNA data after dimensionality reduction; clustering the mRNA data after dimensionality reduction using a k-means clustering algorithm according to the number of cluster clusters to obtain a clustering result; the invention adopts the mclust package to determine the number of cluster clusters, thereby avoiding the limitation of artificially setting the number of cluster clusters, and the clustering result is more accurate; combining the deep autoencoder with the clustering method, speeding up the clustering operation speed, and effectively improving the clustering effect, and has good application prospects.

Description

Technical Field

[0001] The present invention belongs to the field of gene expression data analysis, and in particular relates to a gene expression data clustering method based on a deep autoencoder. Background Art

[0002] As gene chip detection technology becomes more and more mature, a large amount of gene expression data has been generated, and people can obtain gene expression data from public data sets. Using data mining methods to mine important information from gene expression data, we can explore the mechanism and biomarkers of diseases at the molecular level. Cluster analysis is one of the important technologies of data mining and plays a very important role in the study of gene expression data. Cluster analysis of gene expression data can be divided into three categories: gene clustering, sample clustering, and bidirectional clustering. Among them, gene clustering takes genes as objects and sample expression values ​​as features. It can be used to judge the expression patterns of differential genes under different conditions and cluster genes with the same or similar expression patterns. These genes of the same type may have similar functions, participate in the same metabolic process or exist in the same cell pathway.

[0003] Gene clustering can group genes with similar expression patterns together to discover unknown gene functions and regulatory relationships between genes. In the prior art, there are classic clustering algorithms such as partition-based clustering algorithms, hierarchical clustering algorithms, density-based clustering algorithms, and graph-based clustering algorithms, all of which have their own advantages and disadvantages. At present, the gene clustering algorithms that are widely used for gene expression data include partition-based k-means, hierarchical agglomerative hierarchical and Birch, mixed Gaussian model (GMM), and graph-based spectral clustering. There are also some studies on clustering algorithms for gene expression data that improve the shortcomings of clustering algorithms, but almost all of them require the specification of the number of clustering clusters, so they all use simulated data sets or data sets with known clustering cluster numbers. K-means is a classic partition-based unsupervised static clustering method, which is faster and simpler than the agglomerative hierarchical clustering algorithm, and is applicable to a variety of data sets and has good clustering effects, but it also requires the given number of clustering clusters, and in reality, the number of gene clustering clusters in gene expression data sets is often unknown, which has become a difficulty in cluster analysis. The R language mclust package is specifically used to determine the number of clusters, and it contains 14 algorithms for selecting the number of clusters. However, when the R language mclust package is used for gene clustering of gene expression datasets, it often fails to obtain results due to the high dimensionality of the gene expression dataset and the large amount of computation.

[0004] In summary, there is an urgent need for a gene expression data clustering method that does not require a given number of clusters and has less computational complexity. Summary of the invention

[0005] In view of the shortcomings of the prior art, the present invention proposes a gene expression data clustering method based on a deep autoencoder, the method comprising:

[0006] S1: Obtain gene expression data;

[0007] S2: preprocess the gene expression data to obtain mRNA data;

[0008] S3: Input the mRNA data into the trained deep autoencoder to obtain the reduced-dimensional mRNA data;

[0009] S4: The mclust package was used to determine the number of clusters of mRNA data after dimensionality reduction;

[0010] S5: According to the number of clusters, the k-means clustering algorithm is used to cluster the mRNA data after dimensionality reduction to obtain the clustering results.

[0011] Preferably, preprocessing the acquired gene expression data includes: extracting mRNA data from the gene expression data; performing gene ID conversion, missing value processing and duplicate gene deletion on the mRNA data to obtain preliminarily processed mRNA data; removing genes with too low expression levels from the preliminarily processed mRNA data; and using variance stabilization transformation to standardize the data after removing genes with too low expression levels to obtain preprocessed mRNA data.

[0012] Preferably, training the deep autoencoder comprises: the deep autoencoder comprises an encoder and a decoder;

[0013] S31: using an encoder to perform feature learning on the input mRNA data to obtain mRNA data after dimensionality reduction;

[0014] S32: reconstructing the mRNA data after dimensionality reduction using a decoder to obtain reconstructed mRNA data;

[0015] S33: Calculate the loss function according to the reconstructed mRNA data and the original mRNA data, and return to step S31. When the loss function is minimized, a trained deep autoencoder is obtained.

[0016] Furthermore, the formula for the encoder of the deep autoencoder to process data is:

[0017] y=f θ (x) = f k (f k-1 (...f1(x)))=σ k (σ k-1 (...σ1(wx+b)))

[0018] Among them, x represents the input vector, y represents the mRNA data after dimensionality reduction, and f θ represents the encoder, f i represents the mapping function of the i-th layer in the deep network, σ k represents the activation function of the kth layer, w represents the weight matrix value, and b represents the offset matrix.

[0019] Furthermore, the formula for processing data by the decoder of the deep autoencoder is:

[0020] z=g θ′ (y) = g k (g k-1 (...g1(y)))=σ k (σ k-1 (...σ1(w′y+b′)))

[0021] Among them, y represents the input mRNA data after dimension reduction, z represents the reconstructed gene vector, and g θ′ ={g1,g2,...,g k} represents the decoder, σ k represents the activation function, w′ represents the weight matrix value, and b′ represents the offset matrix.

[0022] Furthermore, the loss function is:

[0023] L=-logp(x|z)=-logp(x|z=g θ′ (f θ (x)))

[0024] Among them, L represents the loss function, x represents the original mRNA data of the input, z represents a reconstructed gene vector of the input x, and f θ (x) represents the mRNA data after dimensionality reduction, g θ′ Represents a decoder.

[0025] Preferably, using the mclust package to determine the number of clusters of the mRNA data after dimensionality reduction includes: setting a range of the initial number of clusters; within the range of the number of clusters, all algorithms in the mclust package calculate BIC values ​​according to different numbers of clusters; and using the number of clusters when the BCI value is maximized as the number of clusters of the mRNA data.

[0026] Furthermore, the calculation formula of BIC value is:

[0027]

[0028] Among them, M represents the Mth algorithm in the mclust package, G represents the number of clusters, and x represents the input mRNA data after dimensionality reduction. represents maximum likelihood estimation, n represents the number of genes, and v represents the number of estimated parameters.

[0029] Preferably, clustering the mRNA data after dimensionality reduction using a k-means clustering algorithm includes:

[0030] S51: define initial cluster centers;

[0031] S52: clustering the mRNA data according to the cluster center to obtain clustered mRNA data;

[0032] S53: Calculate the objective function based on the distance between all mRNA data and the cluster center;

[0033] S54: Update the cluster center and return to step S52. When the objective function is minimized, stop clustering.

[0034] Furthermore, the objective function is:

[0035]

[0036] Among them, S i represents the i-th cluster, μ i Represents cluster S i The mean vector cluster center of is, k represents the number of clusters, and x represents the mRNA data.

[0037] The beneficial effects of the present invention are as follows: the present invention adopts the mclust package to determine the number of clusters, thereby avoiding the limitation of artificially setting the number of clusters, and the clustering result is more accurate; the deep autoencoder is combined with the clustering method, and the deep autoencoder reduces the data dimension through feature learning, accelerates the clustering operation speed, and effectively improves the clustering effect; the present invention does not need to give the number of clusters and has a small amount of calculation, and the clustering result has high accuracy. The high-precision gene clustering result is helpful to better study the unknown functions of genes and the regulatory relationship between genes, and has good application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 It is a flow chart of the gene expression data clustering method based on deep autoencoder in the present invention;

[0039] Figure 2 This is a schematic diagram of the effect of preprocessing a thyroid cancer gene expression dataset in a preferred embodiment of the present invention;

[0040] Figure 3 It is the network model of the deep automatic encoder in the present invention. DETAILED DESCRIPTION

[0041] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0042] The present invention proposes a gene expression data clustering method based on deep autoencoder, such as Figure 1 As shown, the method includes:

[0043] S1: Obtain gene expression data;

[0044] S2: preprocess the gene expression data to obtain mRNA data;

[0045] S3: Input the mRNA data into the trained deep autoencoder to obtain the reduced-dimensional mRNA data;

[0046] S4: The mclust package was used to determine the number of clusters of mRNA data after dimensionality reduction;

[0047] S5: According to the number of clusters, the k-means clustering algorithm is used to cluster the mRNA data after dimensionality reduction to obtain the clustering results.

[0048] A specific embodiment of a gene expression data clustering method based on a deep autoencoder is as follows:

[0049] The gene expression dataset of thyroid cancer in the Cancer Genome Atlas (TCGA) database was used to extract mRNA data. The extracted mRNA data were converted into gene IDs, missing values ​​were processed, and duplicate genes were deleted to obtain the initially processed mRNA data. At this time, the number of genes was reduced from 56,404 to 19,193.

[0050] The mRNA data after preliminary processing were further processed, including: using CPM normalization to remove genes with too low expression, retaining genes with CPM normalization expression greater than 0.5 in at least two samples, and using variance stabilization transformation (VST) to standardize the data. The final gene expression data matrix was 16971x568, with 16971 genes, 568 samples, and an expression data range of 6.29-23.38; the results after preprocessing are as follows Figure 2 shown.

[0051] like Figure 3As shown in the figure, the deep autoencoder network includes an input layer, a hidden layer and an output layer; the hidden layer includes 5 fully connected layers; the encoder of the deep autoencoder includes 3 fully connected layers of the input layer and the hidden layer, and the decoder of the deep autoencoder includes 2 fully connected layers of the output layer and the hidden layer. Both the encoder and the decoder use the ReLU function as the activation function; among them, the number of neurons in the first fully connected layer of the encoder is 128, the second fully connected layer is 64, and the third fully connected layer is 25. The number of fully connected layers in the third layer essentially compresses the input data to this dimension; the decoder fully connected layer can reconstruct the features extracted by the encoder and compare them with the original input data. The number of neurons in the first fully connected layer of the decoder is 64, and the second fully connected layer is 128; the number of neurons in the output layer is 568, which is consistent with the original data size, and is used to compare with the input gene expression data to iteratively optimize the model.

[0052] Training the deep autoencoder includes: the deep autoencoder includes an encoder and a decoder;

[0053] S31: using an encoder to perform feature learning on the input mRNA data to obtain mRNA data after dimensionality reduction;

[0054] S32: reconstructing the mRNA data after dimensionality reduction using a decoder to obtain reconstructed mRNA data;

[0055] S33: Calculate the loss function according to the reconstructed mRNA data and the original mRNA data, and return to step S31. When the loss function is minimized, a trained deep autoencoder is obtained.

[0056] The formula for the encoder of the deep autoencoder to process data is:

[0057] y=f θ (x) = f k (f k-1 (...f1(x)))=σ k (σ k-1 (...σ1(wx+b)))

[0058] Among them, x is the input vector, which represents the expression of a gene in each sample, y is the characteristic representation of this gene, here is the mRNA data after dimensionality reduction, and f θ ={f1,f2,...,f k} is the encoder, f i is the mapping function of the i-th layer in the deep network, σ krepresents the activation function of the kth layer, w is the weight matrix value, that is, the weight value of each feature extracted by the gene in each sample, and b is the offset matrix; in this neural network, the output of the first layer is used as the input of the second layer, and k single hidden layers are stacked as a deep neural network to fit the mapping function f θ .

[0059] The formula for the decoder of the deep autoencoder to process data is:

[0060] z=g θ′ (y) = g k (g k-1 (...g1(y)))=σ k (σ k-1 (...σ1(w′y+b′)))

[0061] Among them, y represents the input mRNA data after dimension reduction, z represents the reconstructed gene vector, and g θ′ ={g1,g2,...,g k} represents the decoder, σ k represents the activation function of the kth layer, w′ represents the weight matrix value, and b′ represents the offset matrix.

[0062] The loss function is:

[0063] L=-logp(x|z)=-logp(x|z=g θ′ (f θ (x)))

[0064] Among them, L represents the loss function, x represents the input original mRNA data, z represents a reconstructed gene vector of input x, and f θ (x) represents the mRNA data after dimensionality reduction, g θ′ Represents a decoder.

[0065] The trained deep autoencoder is used to extract features from the preprocessed mRNA data to obtain the reduced-dimensional mRNA data.

[0066] The mclust package was used to determine the number of clusters of the mRNA data after dimensionality reduction, including: setting the range of the initial number of clusters; within the range of the number of clusters, all algorithms in the mclust package calculated the BIC value according to different numbers of clusters; and the number of clusters that maximized the BCI value was used as the number of clusters of the mRNA data.

[0067] The BIC value is calculated as follows:

[0068]

[0069] Among them, M represents the Mth algorithm in the mclust package, G represents the number of clusters, and x represents the input mRNA data after dimensionality reduction. represents maximum likelihood estimation, n represents the number of genes, and v represents the number of estimated parameters.

[0070] For a given gene expression dataset D, D = (x1, x2, x3, ..., x n ), each gene data object x i is an m-dimensional vector, and the k-means clustering algorithm is used to cluster the reduced-dimensional mRNA data, including:

[0071] S51: define initial cluster centers;

[0072] S52: Cluster the mRNA data according to the cluster center, and divide the n gene data objects into k clusters, that is, S = {S1, S2, S3, ..., S k}, get clustered mRNA data;

[0073] S53: Calculate the objective function based on the distance between all mRNA data and the cluster center;

[0074] S54: Update the cluster center and return to step S52. When the objective function is minimized, stop clustering.

[0075] The formula for updating the cluster center is:

[0076]

[0077] Among them, S i represents the i-th cluster, μ i Represents cluster S i The mean vector of cluster centers, x represents the mRNA data.

[0078] The objective function is:

[0079]

[0080] Among them, E represents the sum of the squared differences between the data points in the cluster and the cluster center, S i represents the i-th cluster, μ i Represents cluster S i The mean vector cluster center of is, k represents the number of clusters, and x represents the mRNA data.

[0081] When the distance metric is Euclidean distance, the objective function E is minimized based on the error square function. When the objective function E converges, the best clustering result is obtained. Since the genes clustered together have similar expression patterns in all samples, this may indicate the co-regulation of these genes. In addition, genes from the same cluster may perform similar cellular functions, which is helpful for annotating newly discovered genes. Users can analyze the clustering results to explore the mechanisms and biomarkers of diseases at the molecular level, helping users obtain more biological information.

[0082] The Davidson-Borgin Index (DBI), Calinski-Harabasz Index (CH Index) and Silhouette Coefficient are used as evaluation indicators to evaluate the present invention. The Davidson-Borgin Index (DBI), CH Index and Silhouette Coefficient are internal evaluation methods of clustering algorithms. The DBI calculates the sum of the average distances within any two categories divided by the distance between the two cluster centers to obtain the maximum value, and then calculates the average; the CH Index measures the compactness within the class by calculating the square sum of the distances between each point in the class and the class center, and measures the separation of the data set by calculating the square sum of the distances between the center points of each class and the center point of the data set. The final CH Index is obtained by the ratio of separation to compactness; the Silhouette Coefficient quantifies the compactness and separation of the clusters according to the Euclidean distance.

[0083] DBI is represented by:

[0084]

[0085] Among them, n is the number of categories, c i is the center of the ith category, σ i is the average distance from all points in category i to the centroid of the cluster; d(c i ,c j ) is the distance between the center points of two categories. The smaller the DBI value is, the smaller the intra-class distance is and the larger the inter-class clustering is, so the clustering effect is better.

[0086] The CH index is expressed as:

[0087]

[0088] Among them, n is the number of clusters, k is the current class, is the trace of the inter-class distance difference matrix, is the trace of the intra-class distance difference matrix, z is the mean of the entire data set, and z j is the jth cluster c j The larger the CH(k) value is, the better the clustering effect is.

[0089] The silhouette coefficient is expressed as:

[0090]

[0091] Among them, s(i) represents the silhouette coefficient of the ith gene; a(i) is the dissimilarity within the first cluster, which represents the average distance between the ith gene and other different genes in the same cluster, reflecting the compactness; b(i) is the dissimilarity between the second clusters, which represents the minimum value of the average distance between the ith gene and all genes in other clusters, reflecting the separation. For a data set, its silhouette coefficient is the average silhouette coefficient of all data points (genes). The value of the silhouette coefficient is between [-1,1], and the larger the value, the better the clustering effect.

[0092] DBI, silhouette coefficient and CH index are used as comprehensive evaluation indicators of clustering methods; the DBI, silhouette coefficient and CH index of the common principal component analysis (PCA) combined with k-means are 1.26, 23059.76 and 0.24 respectively; the DBI, silhouette coefficient and CH index of the deep autoencoder combined with k-means of the present invention are 0.64, 71963.67 and 0.46 respectively; compared with the present invention, the clustering effect of the present invention is better.

[0093] The present invention adopts the mclust package to determine the number of clusters, avoiding the limitation of artificially setting the number of clusters, and the clustering result is more accurate; the deep autoencoder is combined with the clustering method, and the deep autoencoder reduces the data dimension through feature learning, speeds up the clustering operation speed, and effectively improves the clustering effect; the present invention does not need to give the number of clusters and has a small amount of calculation, and the clustering result has high accuracy. The high-precision gene clustering result is helpful to better study the unknown functions of genes and the regulatory relationship between genes, and has good application prospects.

[0094] The above embodiments further illustrate the purpose, technical solutions and advantages of the present invention in detail. It should be understood that the above embodiments are only preferred implementation modes of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made to the present invention within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A gene expression data clustering method based on deep autoencoder, characterized in that: include: S1: Obtain gene expression data; S2: preprocess the gene expression data to obtain mRNA data; S3: Input the mRNA data into the trained deep autoencoder to obtain the reduced-dimensional mRNA data; Training the deep autoencoder includes: the deep autoencoder includes an encoder and a decoder; S31: Use the encoder to perform feature learning on the input mRNA data to obtain the mRNA data after dimensionality reduction; the formula for the encoder to process the data is: y=f θ (x)=f k (f k-1 (...f1(x)))=σ k (s k-1 (...σ1(wx+b))) Among them, x represents the input vector, y represents the mRNA data after dimensionality reduction, and f θ represents the encoder, f i represents the mapping function of the i-th layer in the deep network, σ k represents the activation function of the kth layer, w represents the encoder weight matrix value, and b represents the encoder offset matrix; S32: Use a decoder to reconstruct the mRNA data after dimensionality reduction to obtain reconstructed mRNA data; the formula for the decoder to process data is: z=g θ′ (y)=g k (g k-1 (...g1(y)))=σ k (s k-1 (...σ1(w′y+b′))) Among them, z represents the reconstructed gene vector, g θ′ ={g1,g2,...,g k } represents the decoder, w′ represents the decoder weight matrix value, and b′ represents the decoder offset matrix; S33: Calculate the loss function according to the reconstructed mRNA data and the original mRNA data, and return to step S31. When the loss function is minimized, a trained deep autoencoder is obtained; the loss function is: L=-logp(x|z)=-logp(x|z=g θ′ (f θ (x))) Among them, L represents the loss function, f θ (x) represents the mRNA data after dimensionality reduction; S4: The mclust package was used to determine the number of clusters of mRNA data after dimensionality reduction; S5: According to the number of clusters, the k-means clustering algorithm is used to cluster the mRNA data after dimensionality reduction to obtain the clustering results.

2. A gene expression data clustering method based on deep autoencoder according to claim 1, characterized in that: The preprocessing of the acquired gene expression data includes: extracting mRNA data from the gene expression data; performing gene ID conversion, missing value processing and duplicate gene deletion on the mRNA data to obtain preliminarily processed mRNA data; removing genes with too low expression levels from the preliminarily processed mRNA data; and using variance stabilization transformation to standardize the data after removing genes with too low expression levels to obtain preprocessed mRNA data.

3. A gene expression data clustering method based on deep autoencoder according to claim 1, characterized in that: The mclust package is used to determine the number of clusters of mRNA data after dimensionality reduction, including: setting the range of the initial number of clusters; within the range of the number of clusters, all algorithms in the mclust package calculate BIC values ​​according to different numbers of clusters; and the number of clusters that maximizes the BCI value is used as the number of clusters of the mRNA data.

4. A gene expression data clustering method based on deep autoencoder according to claim 3, characterized in that: The BIC value is calculated as follows: Among them, M represents the Mth algorithm in the mclust package, G represents the number of clusters, and x represents the input mRNA data after dimensionality reduction. represents maximum likelihood estimation, n represents the number of genes, and v represents the number of estimated parameters.

5. A gene expression data clustering method based on deep autoencoder according to claim 1, characterized in that: The k-means clustering algorithm is used to cluster the mRNA data after dimensionality reduction, including: S51: define initial cluster centers; S52: clustering the mRNA data according to the cluster center to obtain clustered mRNA data; S53: Calculate the objective function based on the distance between all mRNA data and the cluster center; S54: Update the cluster center and return to step S52. When the objective function is minimized, stop clustering.

6. A gene expression data clustering method based on deep autoencoder according to claim 5, characterized in that: The objective function is: Among them, S i represents the i-th cluster, μ i Represents cluster S i The mean vector cluster center of is, k represents the number of clusters, and x represents the mRNA data.

Images