A twin clustering method and system for interpreting single-cell RNA sequencing data

Abstract

The application discloses a twin clustering method and system for interpreting single-cell RNA sequencing data, and belongs to the technical field of bioinformatics. The method comprises the following steps: preprocessing single-cell RNA sequencing data to obtain a preprocessed gene expression matrix, and constructing a K-nearest neighbor graph between cells based on the preprocessed gene expression matrix; performing feature space enhancement on the preprocessed gene expression matrix to obtain a perturbation attribute matrix; performing structure space enhancement on the K-nearest neighbor graph between cells to obtain an enhanced graph; and performing information fusion and coding and decoding based on the perturbation attribute matrix and the enhanced graph to obtain a clustering result corresponding to the single-cell RNA sequencing data. The application can effectively solve the problems of insufficient exploration of information between cells, low data noise resistance and poor scalability for large data sets in the existing scRNA-seq data clustering method.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology and relates to a twin clustering method and system for interpreting single-cell RNA sequencing data. Background Technology

[0002] Single-cell RNA sequencing (scRNA-seq) is a groundbreaking advancement in bioinformatics. This technology enables the parallel capture of comprehensive genetic information within a single cell, providing higher-resolution data to reveal cellular heterogeneity and diversity. This is crucial for understanding complex biological systems and disease progression. Cluster analysis, as a primary step in scRNA-seq data analysis, is essential. Grouping individual cells by gene expression patterns reveals the complex characteristics of different cell populations, providing deeper insights into their functions and interactions. Furthermore, cluster analysis of gene expression patterns in large-scale scRNA-seq data facilitates wet experimental design in cell and gene-related research and promotes various downstream biological tasks, such as marker gene identification, cell type annotation, and developmental trajectory prediction.

[0003] For decades, clustering analysis of scRNA-seq data has been a vibrant research area. The inherent characteristics of scRNA-seq data—high dimensionality, high sparsity, high noise, nonlinearity, and frequent missing events—make the relationships between variables complex and nonlinear. In recent years, research has increasingly focused on applying deep learning frameworks to learn representations of scRNA-seq data and classify samples into different clusters. Self-supervised learning methods are a typical example, revealing efficient representations of scRNA-seq data by recreating the input data, such as DCA, scDeepCluster, and scNAME. However, these methods primarily focus on extracting features from individual cells, neglecting the crucial intercellular structural information needed to accurately describe cellular differences. Graph Neural Networks (GNNs), by fully utilizing graph information, have been applied to analyze scRNA-seq data, effectively extracting cell topological embeddings and achieving good performance, such as scGAE, scDSC, and scGNN.

[0004] Despite the encouraging results achieved by GNN-based cell clustering methods, limitations remain in addressing the following challenges: (1) insufficient exploration of intercellular information; (2) inadequate data robustness; and (3) low scalability to large datasets. These challenges have prompted researchers to explore a general scRNA-seq data clustering model for interpreting scRNA-seq data. By gathering fundamental knowledge from raw scRNA-seq data and extending it to various downstream tasks, breakthroughs are expected to be achieved in overcoming these challenges. Summary of the Invention

[0005] The purpose of this invention is to provide a twin clustering method and system (hereinafter referred to as scSCF) for interpreting single-cell RNA sequencing data. This method not only offers higher clustering accuracy but also stronger noise resistance and the ability to process large-scale data, providing a solid data foundation for various downstream biological tasks. This invention effectively addresses the problems of insufficient exploration of intercellular information, low data noise resistance, and poor scalability for large datasets faced by existing scRNA-seq data clustering methods.

[0006] To achieve the above-mentioned objectives, the technical solution of the present invention includes the following:

[0007] A twin clustering method for interpreting single-cell RNA sequencing data, the method comprising:

[0008] Preprocessing of single-cell RNA sequencing data yields a preprocessed gene expression matrix. A K-nearest neighbor graph between cells is constructed based on the preprocessed gene expression matrix;

[0009] The preprocessed gene expression matrix is ​​subjected to feature space enhancement to obtain a perturbation attribute matrix.

[0010] Structural spatial enhancement is performed on the K-nearest neighbor graph between the cells to obtain an enhanced graph; wherein, the enhanced graph includes: an enhanced graph after edge removal. And the enhanced image after diffusion operation A r This represents the enhanced image after edge removal. The adjacency matrix, A d This represents the enhanced graph after the diffusion operation. The adjacency matrix;

[0011] Information fusion and encoding / decoding are performed based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data.

[0012] Furthermore, the preprocessing of the single-cell RNA sequencing data to obtain the preprocessed gene expression matrix includes:

[0013] Filtering single-cell RNA sequencing data, wherein the filtering includes: filtering out genes expressed in fewer than one cell and cells with fewer than one gene expression;

[0014] The filtered data is standardized and transformed using logTPM to filter out any genes that are not counted in the cells.

[0015] Highly expressed genes were screened to construct a pre-processed gene expression matrix.

[0016] Furthermore, structural spatial enhancement is performed on the K-nearest neighbor graph between the cells to obtain an enhanced graph after edge removal. include:

[0017] After calculating the pairwise cosine similarity matrix of the samples based on the preprocessed gene expression matrix, a mask matrix M is generated.

[0018] Calculate the adjacency matrix after removing edges. Where D is the degree matrix of the K-nearest neighbor graph between the cells. Let I be the normalized adjacency matrix of the K-nearest neighbor graph between the cells, and let I represent the identity matrix.

[0019] Generate an enhanced map after edge removal.

[0020] Furthermore, the K-nearest neighbor graph between the cells is structurally enhanced to obtain an enhanced graph after diffusion operation. include:

[0021] The normalized adjacency matrix is ​​obtained using the Personalized PageRank algorithm. Convert to adjacency matrix A d ;

[0022] Generate an enhanced graph after diffusion operation

[0023] Furthermore, based on the perturbation attribute matrix and the enhancement graph, information fusion and encoding / decoding are performed to obtain the clustering results corresponding to the single-cell RNA sequencing data, including:

[0024] Construct a twin information fusion network autoencoder and a FAE decoder based on a denoising ZINB model; wherein each branch of the twin information fusion network autoencoder is equipped with a feature autoencoder and a symmetric graph autoencoder;

[0025] Perturbation attribute matrix As input to the feature autoencoder in each branch, the first feature representation matrix is ​​obtained respectively. Second feature representation matrix

[0026] Perturbation attribute matrix and adjacency matrix A r As input to the symmetric graph autodecoder in one branch, and the perturbation attribute matrix and adjacency matrix A d As input to the symmetric graph autodecoder in another branch, the first symmetric graph representation matrix is ​​obtained. First adjacency matrix and the second symmetric graph representation matrix Second adjacency matrix

[0027] Fusion first feature representation matrix Second feature representation matrix The feature representation matrix Z is obtained. FAE ;

[0028] Fusion First Symmetric Graph Representation Matrix Second symmetric graph representation matrix The symmetric graph representation matrix Z is obtained. SGAE ;

[0029] For the feature representation matrix Z FAE and the symmetric graph representation matrix Z SGAE The initial fusion representation matrix Z is obtained by performing a weighted summation calculation. I ;

[0030] Based on normalized adjacency matrix Merging local structural information between data to enhance the initial fusion representation matrix Z I To obtain the expression matrix Z L ; where, the normalized adjacency matrix This represents the normalized adjacency matrix of the K-nearest neighbor graph among the cells;

[0031] By employing an autocorrelation learning mechanism to leverage the nonlocal relationships between data in the fusion space based on initial information, the representation matrix Z is obtained. G The initial information fusion space is based on the perturbation attribute matrix. Build;

[0032] For the expression matrix Z L and expression matrix Z G A weighted summation is performed to obtain the fusion representation matrix Z;

[0033] Self-supervised clustering based on the optimal transmission strategy is performed on the fusion characterization matrix Z to obtain the clustering results corresponding to the single-cell RNA sequencing data;

[0034] Extract and visualize the distribution of the fusion representation matrix Z to obtain a visual representation of the clustering results.

[0035] Furthermore, the process of training the twin information fusion network autoencoder, the FAE decoder based on the denoised ZINB model, and the self-supervised clustering based on the optimal transmission strategy includes:

[0036] Based on the perturbation attribute matrix Reconstruct the weighted attribute matrix Reconstruct the weighted attribute matrix Normalized adjacency matrix and reconstructing the adjacency matrix Constructing a reconstruction loss function based on feature autoencoders and symmetric graph autoencoders Wherein, the reconstructed weighted attribute matrix It is the output obtained by the feature autoencoder based on the fusion representation matrix Z, the reconstructed weighted attribute matrix and the reconstructed adjacency matrix It is the output of the symmetric graph autoencoder based on the fusion representation matrix Z;

[0037] Based on a two-layer cross-correlation reduction mechanism, a cell-level cross-correlation reduction loss function is constructed. Reducing loss function based on gene-level mutual information Among them, the cell-level mutual information reduction loss function Includes: feature cross-correlation reduction, cell-level cross-correlation reduction, and loss function. Reducing graph cross-correlation reduces cell-level cross-correlation in channels and decreases the loss function. The gene-level mutual information reduction loss function This includes: feature-level cross-correlation reduction loss function to reduce gene-level cross-correlation in the channel. Reducing gene-level cross-correlation in channels and graph cross-correlation reduces the loss function.

[0038] Based on the fusion representation matrix Z and the normalized adjacency matrix Constructing the propagation regularization loss function

[0039] Integrating cell-level cross-correlation information to reduce loss function Gene-level cross-correlation information reduction loss function Reconstruct the loss function and propagation regularization loss function The loss function of the twin information fusion network autoencoder is obtained.

[0040] By assuming that single-cell RNA sequencing data follows a ZINB distribution, an overall negative likelihood loss for the FAE decoder based on a denoising ZINB model is constructed.

[0041] Construct a self-supervised clustering loss function based on the optimal transmission strategy.

[0042] Based on the loss function The overall negative likelihood loss and the clustering loss function Complete the training of the twin information fusion network autoencoder and the FAE decoder based on the denoising ZINB model.

[0043] Furthermore, based on the two-layer cross-correlation reduction mechanism, a cell-level cross-correlation reduction loss function is constructed in the feature cross-correlation reduction channel. Reducing loss function based on gene-level mutual information include:

[0044] By calculating the first feature representation matrix The i-th feature and the second feature representation matrix Cosine similarity between the j-th dimension features To construct matrix S cell ;

[0045] Constructing feature-related cross-correlation reduces cell-level cross-correlation in channels, thus reducing the loss function. Where N represents the latent embedding dimension and I represents the identity matrix;

[0046] Use the Readout function on the first feature representation matrix. Second feature representation matrix Projection is performed to obtain the first gene-level feature representation matrix. Second-level gene-level feature representation matrix

[0047] By calculating the first gene-level feature representation matrix The i-th dimension feature and the second gene-level feature representation matrix Cosine similarity between the j-th dimension features To construct matrix S gene ;

[0048] Constructing feature-related cross-correlation information reduces gene-level cross-correlation information in channels, thus reducing the loss function. Where d is the dimension of the Readout function. Represents the identity matrix.

[0049] Furthermore, based on the two-layer cross-correlation reduction mechanism, a cell-level cross-correlation reduction loss function is constructed in the graph cross-correlation reduction channel. Reducing loss function based on gene-level mutual information include:

[0050] By calculating the first symmetric graph representation matrix The i-th dimension feature and the second symmetric graph representation matrix Cosine similarity between the j-th dimension features To construct matrix T cell ;

[0051] Constructing graph-based cross-correlation information reduces cell-level cross-correlation information in channels, thus reducing the loss function. Where N represents the latent embedding dimension and I represents the identity matrix;

[0052] The Readout function is used to characterize the matrix of the first symmetric graph. Second symmetric graph representation matrix Projection is performed to obtain the first gene-level graph representation matrix. Second-level gene-level graph representation matrix

[0053] By calculating the first gene-level graphical representation matrix The i-th dimension feature and the second gene-level graph representation matrix Cosine similarity between the j-th dimension features To construct matrix T gene ;

[0054] Constructing feature-related cross-correlation information reduces gene-level cross-correlation information in channels, thus reducing the loss function. Where d is the dimension of the Readout function. Represents the identity matrix.

[0055] Furthermore, the overall negative likelihood loss of the FAE decoder based on the denoising ZINB model... in, This represents the pretreated gene expression matrix, π. ij Let μ be the probability of the loss event. ij and dispersion θ ij These represent the negative binomial components.

[0056] Furthermore, the self-supervised clustering loss function constructed based on the optimal transmission strategy... include:

[0057] Obtain the clustering result distribution matrix Q;

[0058] Based on the clustering result distribution matrix Q, a function is used to construct the transportation plan matrix P using optimal transmission;

[0059] By utilizing the Sinkhorn distance and solving for the function of the transportation plan matrix P through entropy constraints, the optimal transportation plan matrix can be obtained.

[0060] Based on the clustering result distribution matrix Q and the optimal transmission plan matrix Constructing a self-supervised clustering loss function in, For the optimal transmission plan matrix The element in q ij represents the elements in the clustering result distribution matrix Q.

[0061] A twin clustering system for interpreting single-cell RNA sequencing data, the system comprising:

[0062] The preprocessing module is used to preprocess single-cell RNA sequencing data to obtain a preprocessed gene expression matrix. A K-nearest neighbor graph between cells is constructed based on the preprocessed gene expression matrix;

[0063] The data augmentation module is used to perform feature space augmentation on the preprocessed gene expression matrix to obtain a perturbation attribute matrix. Structural spatial enhancement is performed on the K-nearest neighbor graph between the cells to obtain an enhanced graph; wherein, the enhanced graph includes: an enhanced graph after edge removal. And the enhanced image after diffusion operation A r This represents the enhanced image after edge removal. The adjacency matrix, A d This represents the enhanced graph after the diffusion operation. The adjacency matrix;

[0064] The twin information fusion module is used to perform information fusion and encoding / decoding based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data.

[0065] Compared with existing technologies, this invention, based on a Siamese graph autoencoder, provides a general clustering framework for interpreting scRNA-seq data, effectively addressing the challenges faced by existing clustering methods. It not only boasts higher clustering accuracy but also stronger noise resistance and the ability to handle large-scale data, providing a solid data foundation for various downstream biological tasks. This invention effectively solves the problems of insufficient exploration of intercellular information, low data noise resistance, and poor scalability for large datasets faced by existing scRNA-seq data clustering methods. Specifically, it has the following advantages:

[0066] (1) This invention has faster training speed and more accurate clustering results;

[0067] (2) The present invention has stronger noise resistance and the ability to process large-scale datasets;

[0068] (3) In the process of clustering scRNA-seq data, the present invention can effectively learn the data representation matrix and has a powerful data representation capability;

[0069] (4) It can be invented with good data visualization capabilities;

[0070] (5) The characterization matrix learned in this invention can support downstream tasks in a variety of biological fields. Attached Figure Description

[0071] Figure 1 A framework diagram for interpreting twin clustering methods for single-cell RNA sequencing data.

[0072] Figure 2 A flowchart of the overall process for interpreting twin clustering methods for single-cell RNA sequencing data.

[0073] Figure 3 Visualization results of this invention and four competing methods on human liver cells. Detailed Implementation

[0074] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0075] The present invention provides a twin clustering method for interpreting single-cell RNA sequencing data, such as... Figure 1 and Figure 2 As shown, it includes three steps:

[0076] Step 1: Data Standardization and Preprocessing. The raw scRAN-seq data is preprocessed using the Python SCANPY package: First, this invention filters out genes expressed in fewer than one cell, and also filters out cells with fewer than one gene expression. Next, this invention standardizes the data and transforms it using logTPM, filtering out any genes not counted in any cells. Then, the gene expression matrix of the top 1000 highly expressed genes is selected. A K-nearest neighbor graph between cells is constructed based on the preprocessed gene expression matrix. Finally, the preprocessed gene expression matrix and the K-nearest neighbor graph between cells are used as input to the core algorithm of scSCF.

[0077] Step 2: Core Algorithm. This part consists of three key components: (1) Data Augmentation Module: Enriching the data foundation for representation learning through data augmentation in the feature space and structure space respectively; (2) Siamese Information Fusion Module: First, reducing the correlation between cells and genes in the two channels with shared weights respectively. Second, integrating the complex information inherent in cells through a cross-modal information dynamic fusion mechanism; (3) Self-Supervised Learning Module: Constructing a self-supervised clustering loss function based on the optimal transmission strategy; (4) Joint Optimization: Jointly optimizing the cell representation space through a self-supervised module based on the clustering loss, ZINB loss, and loss from the Siamese Information Fusion Module based on the optimal transmission.

[0078] Step 3. Visualize clustering results and evaluate the model.

[0079] This invention achieves efficient clustering of scRNA-seq data through the above-mentioned technical means, solves the scalability and noise resistance problems faced by existing methods when processing large-scale data, and provides a reliable data foundation for downstream biological tasks.

[0080] Specifically, the contents of each implementation step of the present invention are as follows.

[0081] 1. Data standardization preprocessing:

[0082] The raw scRNA-seq counting matrix was obtained and preprocessed using the Python SCANPY package. First, this invention filters out genes expressed in fewer than one cell and cells with fewer than one gene expression count. Next, the data was standardized and transformed using logTPM to filter out genes not counted in any cell. Then, the expression matrices of the top 1000 highly expressed genes were selected to obtain the final preprocessed gene expression matrix. x ij(1≤i≤N, 1≤j≤D) represents the expression status of the j-th gene in the i-th cell; finally, a K-nearest neighbor (KNN) graph is constructed based on the preprocessed gene expression matrix. Where V = {v1, v2, ... v} n Let} and E represent the set of nodes and the set of edges between nodes, respectively. Each node is a single cell, and the edges represent the similarity between each cell and other cells.

[0083] Mapping method: A KNN graph is constructed using the Pearson correlation coefficient, characterized by the preprocessed gene expression matrix X and adjacency matrix A = (a ij ) N×N , here (v i ,v j When )∈E, a ij =1, otherwise a ij =0. The corresponding degree matrix is and Through D -1 (A+I) normalizes the adjacency matrix A to obtain in It is an identity matrix.

[0084] Assumption It is a preprocessed gene expression matrix The representation matrix, c = [c1, c2, ..., c N ] represents the cluster assignment. This method learns both the representation matrix Z and the cluster assignment c simultaneously.

[0085] 2. Data Augmentation Module:

[0086] Without affecting the information conveyed by the original scRNA-seq data, this invention creates novel and reasonable data that helps enhance the robustness and generalization ability of the model. The data augmentation module is divided into feature space augmentation and structure space augmentation.

[0087] 2.1 Feature Space Enhancement: Feature space enhancement, also known as feature perturbation, involves perturbing the feature representation of nodes in the preprocessed gene expression matrix. For example, adding Gaussian noise can yield a perturbed attribute matrix. As shown in Formula 1.

[0088]

[0089] Where ⊙ represents the Hadamard product, and N represents the distribution from a Gaussian distribution. The random noise matrix, It is a preprocessed gene expression matrix.

[0090] 2.2 Structure Space Enhancement: Structure space enhancement operates on graph structures by removing edges or using a diffusion matrix to generate a global view.

[0091] (1) Edge Removal: To achieve edge removal, a mask matrix is ​​generated based on the pairwise cosine similarity matrix of samples computed in the latent space. Edges with the lowest similarity values ​​(10%) are manually removed. Finally, the edges are normalized to obtain the edge-removed adjacency matrix. The calculation process is as follows remember To achieve an enhanced image after edge removal.

[0092] (2) Graph diffusion: This invention uses Personalized PageRank (PPR) to normalize the adjacency matrix. Convert to graph diffusion adjacency matrix A d α is the transmission probability, allowing messages to propagate through higher-order neighbor nodes. In this invention, α is set to 0.1. The diffusion adjacency matrix A... d Calculated as remember Enhanced graph after diffusion operation.

[0093] 3. Twin Information Fusion Module

[0094] This invention introduces a Siamese Information Fusion Mosule (SIFM), which includes an autoencoder, a feature decoder based on a denoising ZINB model, a dual cross-correlation reduction mechanism, and a cross-modal dynamic information fusion mechanism. It aims to fully explore intercellular information and learn discriminative sample representations across features.

[0095] 3.1 Automatic Encoder:

[0096] As shown in the framework diagram of the scSCF algorithm, the encoder part of this invention adopts a Siamese network structure. Through parameter sharing, each branch is equipped with a Feature Autoencoder (FAE) and a Symmetric Graph Autoencoder (SGAE) to extract attribute and structural information, respectively. In the decoder part, this invention uses a clustering-oriented latent embedding Z and a normalized adjacency matrix. As input, the inputs of the two sub-networks are reconstructed through the decoders of FAE and SGAE, respectively. Unlike traditional autoencoders, the FAE decoder has three output layers, which use a ZINB-based distribution to model highly sparse and overly dispersed scRNA-seq data. Furthermore, SGAE utilizes a symmetric architecture to reconstruct the weighted attribute matrix and adjacency matrix. This approach enhances SGAE's ability to capture the structural and attribute features of the graph, thereby achieving more accurate and comprehensive information reconstruction and significantly improving the network's generalization ability.

[0097] Feature autoencoder:

[0098] A feature autoencoder is proposed, which learns the discriminative representation matrix of scRNA-seq data by stacking three dense networks in the encoder and decoder. This is achieved using a perturbation attribute matrix. As input, reconstruct the weighted attribute matrix. This invention represents the feature encoder and decoder functions as follows:

[0099]

[0100] Where b and b′ represent the bias vectors of the encoder and decoder, respectively, W FAE and W′ FAE The weights of the layers are learned from network training, f enc (·) and f dec (·) represent the encoder and decoder, respectively.

[0101] Symmetrical graph autoencoder:

[0102] Symmetric graph autoencoders learn from perturbation attribute matrices and the adjacency matrix A with edges removed r (or the diffusion adjacency matrix A) d The network learns a representation matrix containing inter-cell topological information. During the decoding process, the weighted attribute matrix needs to be reconstructed simultaneously. and adjacency matrix Furthermore, the formula for one layer in the encoder and decoder is:

[0103]

[0104] in and Let represent the learnable parameters of the l-th encoder layer and the h-th decoder layer, respectively. σ is a non-linear activation function, such as ReLU or Tanh.

[0105] To minimize the reconstruction loss function on both FAE and SGAE, this invention aims to minimize the hybrid loss, as shown in Equation 4:

[0106]

[0107] Where η is a predefined hyperparameter, set to 0.1. It describes the joint mean squared error reconstruction loss of node attributes and graph structure, which is minimized. To optimize the objective function.

[0108] 3.2 FAE Decoder Based on Denoising ZINB Model:

[0109] The ZINB distribution can model highly sparse and overly dispersed counting data. To explain the characteristics of scRNA-seq data, this invention assumes that the data follows a ZINB distribution, i.e.:

[0110]

[0111] in, This represents the preprocessed gene expression matrix, which is also the feature matrix input to the encoder; π ij Let μ be the probability of a missing event (i.e., the probability that a non-zero gene expression value is observed to be 0). ij and dispersion θ ij They represent the negative binomial components, and δ0(·) is the Dirac function (when X... ij When δ = 0, δ0(X) ij =1, otherwise 0).

[0112] The decoder has three output layers, which are used to estimate the probability values ​​of parameters π, μ, and θ, respectively. Mean estimate The variance estimates are as follows:

[0113]

[0114] W π W μ W θ These are the corresponding weights, the size factor s i This is the ratio of the total number of cells to the median. Finally, this invention minimizes the overall negative likelihood loss of the feature decoder based on the ZINB model, i.e.

[0115]

[0116] 3.3 Dual-relational information reduction mechanism:

[0117] The learning process of the two-layer cross-correlation information reduction mechanism includes two channels: Feature Cross-correlation Information Reduction (FCR) and Graph Cross-correlation Information Reduction (GCR). Each channel operates identically except for the input data, including Cell-level Cross-correlation Information Reduction (cell-CR) and Gene-level Cross-correlation Information Reduction (gene-CR). The aim is to minimize the deviation between the cross-correlation matrix and the ideal identity matrix, reducing redundant information between nodes in the latent space, thereby leveraging structural information to improve the performance of the scRNA-seq data clustering algorithm.

[0118] Next, we will take the Feature Cross-correlation Information Reduction (FCR) channel as an example for explanation:

[0119] (1) Reduced cell-level cross-correlation information:

[0120] The cell-CR learning process consists of two steps. First, given the feature representation matrices learned by the feature encoder under the two views... and This invention firstly through Calculate the cosine similarity at the cell dimension. This represents the cosine similarity between the i-th dimension feature in the first view and the j-th dimension feature in the second view.

[0121] Then, this method makes the cell-CR matrix S cell Approximately equal to the identity matrix The target loss function for cell-CR, used as the FCR channel, is as follows:

[0122]

[0123] Where N is the latent embedding dimension (also the number of samples). By ensuring S cellThe diagonal elements are approximately equal to 1, and the off-diagonal elements are approximately equal to 0. This invention can guarantee the embedding alignment of each cell in two views and minimize the embedding consistency of different cells in different views. This helps the method reduce redundant information in the latent space, thereby learning a more discriminative representation.

[0124] (2) Reduced gene-level cross-correlation information:

[0125] In addition to constructing a meaningful representation matrix by reducing cellular correlations between two views, this invention also refines information correlations at the gene level. Specifically, gene-CR is performed in three steps.

[0126] The first step in this method is to use the Readout function. Cell characterization matrix in two views and Projected onto the gene-level characterization matrix and The calculation formula is as follows:

[0127]

[0128] In the second step, the invention again... Calculate the gene-CR matrix S gene The elements in the text, each element This represents the feature similarity between the i-th dimension feature in one view and the j-th dimension feature in another view.

[0129] The third step, according to the objective equation of the gene-CR target loss function of the FCR channel listed in Formula 10, is to make the gene-CR matrix S gene equal to the identity matrix

[0130]

[0131] Where d is the latent embedding dimension. In Equation 10, The emphasis is on making similar features even closer together. It emphasizes separating different characteristics.

[0132] Subsequently, similar to the operation of the FCR channel, this invention is based on a symmetric graph representation matrix. and symmetric graph representation matrix Two target loss functions for the GCR channel are obtained. and

[0133] This invention linearly combines the objective functions of FCR and GCR to obtain... and

[0134] 3.4 Cross-modal dynamic information fusion mechanism:

[0135] To filter redundant features and retain more discriminative features from scRNA-seq data, this invention introduces a cross-modal dynamic information fusion strategy, enabling the scSCF network to learn more discriminative representations to avoid collapse, thereby improving clustering performance.

[0136] First, this invention linearly combines the representation matrices learned by FAE and SGAE to obtain an initial fused representation matrix. As shown below:

[0137]

[0138] Where β is a learnable coefficient, which is initialized to 0.5 in this invention. It selectively determines the importance of the two information sources based on the characteristics of their respective datasets.

[0139] This invention employs graph convolution-like operations (i.e., message passing) to process combined information, enhancing the initial fusion representation matrix Z by incorporating local structural information from the data. I ,Right now

[0140] Next, this invention introduces an autocorrelation learning mechanism to utilize initial information to fuse nonlocal relationships between samples in the space. Specifically, this invention first uses... Calculate the normalized autocorrelation matrix S. Then, using S as coefficients, recombine Z by considering the global correlation between samples. L Z G =SZ L .

[0141] Finally, this invention employs skip connections to facilitate the smooth transfer of information within the fusion mechanism, thus obtaining a cluster-oriented latent representation. It can be used for clustering:

[0142] Z = γZ G +Z L (Formula 12)

[0143] Here, γ is a scale parameter, which is initialized to 0 in this invention and its weights are learned during network training.

[0144] In this way, redundant features can be filtered out while more discriminative features are retained in the latent space, enabling scSCF to learn meaningful representations, avoid collapse, and thus improve clustering performance.

[0145] 3.5 Propagation Regularization

[0146] This invention introduces propagation regularization to alleviate the oversmoothing problem during network training, as shown in the following formula:

[0147]

[0148] Where JSD(·) represents the Jensen-Shannon divergence. Using... Even in shallow network architectures, networks can capture long-distance dependencies, thereby reducing oversmoothing as information propagates deeper within the framework.

[0149] In summary, the target of the SIFM module is calculated using the following formula:

[0150]

[0151] Where δ is a trade-off parameter, set to 0.001.

[0152] 4. Self-Supervised Learning Module

[0153] To better adapt to the complex structure of scRNA-seq data, especially its high dimensionality and high sparsity, this invention employs an optimal transmission strategy to address this complexity. It uses a Student's t-distribution as the kernel to measure the embedding point h. i and cluster center c j The similarity between them.

[0154]

[0155] Where h i =f(x) i )∈Z corresponds to The representation of , ε=1 is the degree of freedom of the student's t-distribution, q ij This is the probability of assigning sample i to cluster j (i.e., soft assignment). In this invention, Q = [q...] ij ] is defined as the distribution of all sample assignments.

[0156] After obtaining the clustering result distribution Q, the learning objective is to optimize the data representation by learning high-confidence assignments. Most traditional methods employ a clustering-guided loss function to force the generated sample representations to have minimal distortion relative to pre-learned cluster centers. This approach aims to minimize the discrepancy between the data representation and the cluster centers, thereby improving cluster cohesion. However, to prevent degenerate solutions that assign all data points to a single (arbitrary) label, this invention introduces a constraint aligning the label distribution with the mixing ratio. This approach ensures improved clustering accuracy and a balanced contribution of each data point to the loss calculation. Therefore, this invention constructs the target probability matrix P by solving the following optimal transfer:

[0157]

[0158] Here, this invention treats the target distribution P as the transportation plan matrix in optimal transportation theory, and -logQ as the relevant cost matrix. Constraints are applied to P. T 1 N =N π , where π represents the proportion of points in each cluster, which can be estimated using intermediate clustering results. This step successfully adds the constraint that the cluster distribution must be consistent with the mixing proportion. Considering the computational burden of direct optimization, this invention utilizes the Sinkhorn distance and performs fast optimization through entropy constraints. The optimization process with entropy-constrained Lagrange multipliers is as follows:

[0159]

[0160] Where H is the entropy function and λ is the smoothing parameter controlling cluster balance, set to 5. At this point, the uniqueness of P is guaranteed, and it can be solved using fixed-point iteration using Equation 17:

[0161]

[0162] Where t represents the number of iterations, and in each iteration u (t) =1 N / (Q λ v (t-1) ) and v (t) =N π / (Q λ u (t) The starting point is v. (0) =1 N Through fixed-point iteration, this invention can obtain the optimal transmission plan matrix. During training, the present invention is fixed. and make Q with To maintain consistency, the clustering loss function is expressed as:

[0163]

[0164] The goal of the learning process is to minimize the clustering loss.

[0165] 5. Joint optimization

[0166] The overall optimization objective of the proposed method consists of three parts: SIFM loss, ZINB loss, and clustering loss.

[0167]

[0168] Where ρ and σ are adjustment parameters for balancing different losses, and in this invention, ρ is initialized to 10 and σ is initialized to 2.

[0169] The training process of scSCF includes three steps: (i) pre-training FAE and SGAE for at least 50 epochs respectively by minimizing the reconstruction loss function; (ii) integrating the two sub-networks into a unified framework and then training for 100 epochs; (iii) training the entire network for at least 50 epochs under the guidance of Equation 19 until convergence. Since the number of cell clusters k is unknown in practical applications, this invention uses the K-means method to determine the optimal value of k on the latent representation Z for clustering.

[0170] 6. Model Evaluation

[0171] To evaluate the effectiveness of the method, this invention employs three performance metrics widely used in unsupervised learning scenarios: Accuracy (ACC), Normalized Mutual Information (NMI), and Adjusted Rand Index (ARI). Assume U represents the cell type with true labels, and V represents the clustering result predicted by the model. Higher values ​​for these evaluation metrics indicate greater consistency between the two cluster assignments U and V, and thus better clustering performance.

[0172] ACC is defined as measuring the best match between two cluster assignments U and V. If there is a data point i, l i It's a real tag, u i If it's an assignment based on a clustering algorithm, then ACC is defined as:

[0173]

[0174] Here, m covers all possible one-to-one mappings between U and V, and the optimal mapping relative to n cells can be found efficiently using the Hungarian algorithm.

[0175] NMI is a function that measures the consistency between the predicted labels and the true labels of n cells. Specifically...

[0176]

[0177] ARI is a function used to evaluate the similarity between two assignments. Given a set X of size n and two clustering results U = {u1, u2, ..., u...} r}, V={v1,v2,…,v s}, and n ij =|U i ∩V j |,i=1,2,…r,j=1,2,…,s, The formula for calculating ARI is as follows:

[0178]

[0179] As can be seen from the above, this invention provides a general clustering framework for interpreting scRNA-seq data, which not only has higher clustering accuracy, but also stronger noise resistance and the ability to process large-scale data, providing a solid data foundation for various downstream biological tasks.

[0180] 7. Visualization of clustering results

[0181] To obtain an intuitive biological interpretation of the clustering results, this invention uses t-SNE to extract and visualize the distribution of the learned clustering embeddings Z in two-dimensional space. Each point on the visualization represents a cell, and each color represents a predicted cell type.

[0182] In summary, this invention provides a general clustering framework for interpreting scRNA-seq data, achieving efficient clustering of scRNA-seq data, solving the scalability and noise resistance problems faced by existing methods when processing large-scale data, and providing a reliable data foundation for downstream biological tasks.

[0183] The following is a specific experiment illustrating the twin clustering framework provided by this invention for interpreting single-cell RNA sequencing data.

[0184] 1. Dataset Overview

[0185] To evaluate the performance of scSCF, seven real-world scRNA-seq datasets of known cell types downloaded from the GEO database were clustered. Table 1 summarizes the details of the seven real-world datasets. These datasets originated from different sequencing methods. Of the seven datasets, four were from human samples and three were from mouse samples. These datasets cover a wide range of cell types, including retinal, lung, liver, kidney, pancreas, and peripheral blood mononuclear cells.

[0186] Table 1 Summary of scRNA-seq datasets

[0187] Shekhar mouse retina cells Drop-seq 27499 13166 19 Macosko mouse retina cells Drop-seq 14653 11422 39 QS mouse lung cells Smart-seq2 1676 23341 11 CITE_CMBC 10X Genomics 8617 2000 15 Human liver cells 10X Genomics 8444 4999 11 Human kidney cells 10X Genomics 5685 25125 11 Human pancreas cells CEL-seq2 2122 19046 9

[0188] 2. Experimental Environment and Parameter Settings

[0189] In the experiment, this invention implemented scSCF in Python 3.7 using PyTorch. See the detailed steps for the parameter settings of this invention.

[0190] 3. Evaluation Indicators

[0191] In the experiment, the present invention used three metrics for performance evaluation: accuracy (ACC), normalized mutual information (NMI), and adjusted Rand index (ARI).

[0192] 4. Analysis of Experimental Results

[0193] (1) Here, the model performance of scSCF is compared with that of five competing methods (scDeepCluster, scNAME, scDSC, scGNN, scCDCG), and experiments are conducted on seven real datasets. The experimental results are as follows:

[0194] Table 2 shows the clustering performance of all models on ten datasets (average values; experimental results for scSCF are bolded).

[0195]

[0196]

[0197] Based on the results, the present invention can observe that:

[0198] 1) The proposed method, scSCF, outperforms competing methods in all three metrics across all datasets. Unlike previous methods, scSCF learns more discriminative scRNA-seq data representations through reduced cross-correlation at both the cellular and gene levels, as well as cross-modal dynamic information fusion, thereby improving clustering performance.

[0199] 2) Referring to Table 1, taking Shekhar mouse retina cells and Macosko mouse retina cells as examples, scSCF performs exceptionally well on large-scale datasets, achieving a significant performance improvement. scSCF effectively utilizes the rich information inherent in scRNA-seq data and alleviates the over-smoothing problem, especially on large-scale datasets.

[0200] 3) This invention significantly improves upon existing technologies, demonstrating superior representation learning capabilities on scRNA-seq data while successfully resolving the representation collapse problem. Specifically, by thoroughly extracting higher-order structural information between cells and integrating gene expression attributes, scSCF effectively learns discriminative features from scRNA-seq data. The method of this invention enhances data interpretability and provides strong support for revealing complex biological relationships between cells, thereby promoting further development in biological research.

[0201] (2) To obtain an intuitive biological interpretation of the clustering results, this invention uses t-SNE to extract and visualize the distribution of the learned clustering embeddings Z in two-dimensional space. Each point on the visualization represents a cell, and each color represents a predicted cell type. From... Figure 3 As can be seen, regardless of cluster size, scSCF effectively distinguishes different cell types, forming clear boundaries between cell populations. In contrast, other methods fail to clearly delineate clusters. These visualizations confirm that scSCF can detect various cell types and cluster similar cells together while separating different cell types, thereby improving the accuracy and biological interpretability of cluster analysis in scRNA-seq data.

[0202] The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Those skilled in the art can modify or make equivalent substitutions to the technical solutions of the present invention without departing from the scope of the present invention. The scope of protection of the present invention should be determined by the claims.

Claims

1. A twin clustering method for interpreting single-cell RNA sequencing data, characterized in that, The method includes: Preprocessing of single-cell RNA sequencing data yields a preprocessed gene expression matrix. And based on the preprocessed gene expression matrix, construct intercellular [cell expression]... Nearest neighbor graph; Feature space enhancement is performed on the preprocessed gene expression matrix to obtain the perturbation attribute matrix. ; For the intercellular The nearest neighbor graph is subjected to structural space enhancement to obtain an enhanced graph; wherein, the enhanced graph includes: an enhanced graph after edge removal. And the enhanced image after diffusion operation , This represents the enhanced image after edge removal. The adjacency matrix, This represents the enhanced graph after the diffusion operation. The adjacency matrix; Information fusion and encoding / decoding are performed based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data; Specifically, information fusion and encoding / decoding are performed based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data, including: Construct a twin information fusion network autoencoder and a FAE decoder based on a denoising ZINB model; wherein each branch of the twin information fusion network autoencoder is equipped with a feature autoencoder and a symmetric graph autoencoder; Perturbation attribute matrix As input to the feature autoencoder in each branch, the first feature representation matrix is ​​obtained respectively. Second feature representation matrix ; Perturbation attribute matrix and adjacency matrix As input to the symmetric graph autodecoder in one branch, and the perturbation attribute matrix and adjacency matrix As input to the symmetric graph autodecoder in another branch, the first symmetric graph representation matrix is ​​obtained. First adjacency matrix and the second symmetric graph representation matrix Second adjacency matrix ; Fusion first feature representation matrix Second feature representation matrix The feature representation matrix is ​​obtained. ; Fusion First Symmetric Graph Representation Matrix Second symmetric graph representation matrix The symmetric graph representation matrix is ​​obtained. ; For the feature representation matrix and symmetric graph representation matrix The initial fusion representation matrix is ​​obtained by performing a weighted summation calculation. ; Based on normalized adjacency matrix Merging local structural information between data to enhance the initial fusion representation matrix To obtain the expression matrix ; where, the normalized adjacency matrix Indicating the intercellular relationship Normalized adjacency matrix of the nearest neighbor graph; An autocorrelation learning mechanism is used to utilize the nonlocal relationships between data in the fusion space based on initial information to obtain the representation matrix. The initial information fusion space is based on the perturbation attribute matrix. Build; For the expression matrix and expression matrix The fusion representation matrix is ​​obtained by performing a weighted summation calculation. ; For the fusion representation matrix Self-supervised clustering based on the optimal transmission strategy is performed to obtain the clustering results corresponding to the single-cell RNA sequencing data; Extract and visualize the fusion representation matrix The distribution of the clustering results is visualized.

2. The method according to claim 1, characterized in that, The preprocessing of single-cell RNA sequencing data to obtain a preprocessed gene expression matrix includes: Filtering single-cell RNA sequencing data, wherein the filtering includes: filtering out genes expressed in fewer than 1 cell and cells with fewer than 1 gene expression; The filtered data is standardized and transformed using logTPM to filter out any genes that are not counted in the cells. Highly expressed genes were screened to construct a pre-processed gene expression matrix.

3. The method according to claim 1, characterized in that, For the intercellular Structural space augmentation is performed on the nearest neighbor graph to obtain the augmented graph after edge removal. ,include: After calculating the pairwise cosine similarity matrix of the samples based on the preprocessed gene expression matrix, a mask matrix is ​​generated. ; Calculate the adjacency matrix after removing edges ;in, For the intercellular The degree matrix of the nearest neighbor graph. For the intercellular The normalized adjacency matrix of the nearest neighbor graph. Represents the identity matrix; Generate an enhanced map after edge removal. .

4. The method according to claim 1, characterized in that, For the intercellular The nearest neighbor graph is augmented in the structural space to obtain the augmented graph after the diffusion operation. ,include: The normalized adjacency matrix is ​​obtained using the Personalized PageRank algorithm. Convert to adjacency matrix ; Generate an enhanced graph after diffusion operation .

5. The method according to claim 1, characterized in that, The training processes for the twin information fusion network autoencoder, the FAE decoder based on the denoised ZINB model, and the self-supervised clustering based on the optimal transmission strategy include: Based on the perturbation attribute matrix Reconstruct the weighted attribute matrix Reconstruct the weighted attribute matrix Normalized adjacency matrix and reconstructing the adjacency matrix Construct a reconstruction loss function based on feature autoencoders and symmetric graph autoencoders. ; wherein, the reconstructed weighted attribute matrix It is a feature autoencoder based on a fusion representation matrix The obtained output is the reconstructed weighted attribute matrix. and the reconstructed adjacency matrix It is a symmetric graph autoencoder based on a fusion representation matrix The output obtained; Based on a two-layer mutual information reduction mechanism, a cell-level mutual information reduction loss function is constructed. Reducing loss function based on gene-level cross-correlation information The cell-level mutual information information reduction loss function is mentioned above. Includes: feature cross-correlation reduction, cell-level cross-correlation reduction, and loss function. Reducing graph cross-correlation reduces cell-level cross-correlation in channels and decreases the loss function. The gene-level cross-correlation information reduction loss function This includes: feature-level cross-correlation reduction loss function to reduce gene-level cross-correlation in the channel. Reducing gene-level cross-correlation in channels and graph cross-correlation reduces the loss function. ; Based on the fusion representation matrix and normalized adjacency matrix Construct a propagation regularization loss function ; Integrating cell-level cross-correlation information to reduce loss function Gene-level cross-correlation information reduction loss function Reconstruct the loss function and propagation regularization loss function The loss function of the twin information fusion network autoencoder is obtained. ; By assuming that single-cell RNA sequencing data follows a ZINB distribution, the overall negative likelihood loss of the FAE decoder based on a denoising ZINB model is constructed. ; Construct a self-supervised clustering loss function based on the optimal transmission strategy. ; Based on the loss function The overall negative likelihood loss and the clustering loss function Complete the training of the twin information fusion network autoencoder and the FAE decoder based on the denoising ZINB model.

6. The method according to claim 5, characterized in that, The method is based on a two-layer mutual information reduction mechanism, which constructs a cell-level mutual information reduction loss function in the feature mutual information reduction channel. Reducing loss function based on gene-level cross-correlation information ,include: By calculating the first feature representation matrix The Middle 1st-dimensional feature and second-dimensional feature representation matrix The Middle Cosine similarity between features To construct a matrix ; Constructing feature-related cross-correlation reduces cell-level cross-correlation in channels, thus reducing the loss function. ;in, Represents the potential embedding dimension. Represents the identity matrix; use Function on the first characteristic matrix Second feature representation matrix Projection is performed to obtain the first gene-level feature representation matrix. Second-level gene-level feature representation matrix ; By calculating the first gene-level feature representation matrix The Middle 1D features and second-level gene-level feature representation matrix The Middle Cosine similarity between dimensional features To construct a matrix ; Constructing feature-related cross-correlation information reduces gene-level cross-correlation information in channels, thus reducing the loss function. ;in, for Functional dimension Represents the identity matrix.

7. The method according to claim 5, characterized in that, The method is based on a two-layer cross-correlation reduction mechanism, which constructs a cell-level cross-correlation reduction loss function in the graph cross-correlation reduction channel. Reducing loss function based on gene-level cross-correlation information ,include: By calculating the first symmetric graph representation matrix The Middle 3D features and second symmetric graph representation matrix The Middle Cosine similarity between dimensional features To construct a matrix ; Constructing graph-based cross-correlation information reduces cell-level cross-correlation information in channels, thus reducing the loss function. ;in, Represents the potential embedding dimension. Represents the identity matrix; use Function representation matrix of the first symmetric graph Second symmetric graph representation matrix Projection is performed to obtain the first gene-level graph representation matrix. Second-level gene-level graph representation matrix ; By calculating the first gene-level graph representation matrix The Middle 3D features and second-level gene-level graph representation matrix The Middle Cosine similarity between dimensional features To construct a matrix ; Constructing feature-related cross-correlation information reduces gene-level cross-correlation information in channels, thus reducing the loss function. ;in, for Functional dimension Represents the identity matrix.

8. The method according to claim 1, characterized in that, The overall negative likelihood loss of the FAE decoder based on the denoising ZINB model. ;in, This represents the gene expression matrix after pretreatment. The probability of a loss event, mean and dispersion These represent the negative binomial components; The self-supervised clustering loss function is constructed based on the optimal transmission strategy. ,include: Obtain the distribution matrix of clustering results ; Combining the clustering result distribution matrix Using optimal transmission to construct the transportation planning matrix The function; Using Sinkhorn distance and applying entropy constraints to the transportation plan matrix Solving the function yields the optimal transmission plan matrix. ; Based on the distribution matrix of clustering results and optimal transmission plan matrix Constructing a self-supervised clustering loss function ;in, For the optimal transmission plan matrix The elements in Distribution matrix of clustering results The elements in.

9. A twin clustering system for interpreting single-cell RNA sequencing data, characterized in that, The system includes: The preprocessing module is used to preprocess single-cell RNA sequencing data to obtain a preprocessed gene expression matrix. And based on the preprocessed gene expression matrix, construct intercellular [cell expression]... Nearest neighbor graph; The data augmentation module is used to perform feature space augmentation on the preprocessed gene expression matrix to obtain a perturbation attribute matrix. ; for the intercellular The nearest neighbor graph is subjected to structural space enhancement to obtain an enhanced graph; wherein, the enhanced graph includes: an enhanced graph after edge removal. And the enhanced image after diffusion operation , This represents the enhanced image after edge removal. The adjacency matrix, This represents the enhanced graph after the diffusion operation. The adjacency matrix; The twin information fusion module is used to perform information fusion and encoding / decoding based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data; Specifically, information fusion and encoding / decoding are performed based on the perturbation attribute matrix and the enhancement graph to obtain the clustering results corresponding to the single-cell RNA sequencing data, including: Construct a twin information fusion network autoencoder and a FAE decoder based on a denoising ZINB model; wherein each branch of the twin information fusion network autoencoder is equipped with a feature autoencoder and a symmetric graph autoencoder; Perturbation attribute matrix As input to the feature autoencoder in each branch, the first feature representation matrix is ​​obtained respectively. Second feature representation matrix ; Perturbation attribute matrix and adjacency matrix As input to the symmetric graph autodecoder in one branch, and the perturbation attribute matrix and adjacency matrix As input to the symmetric graph autodecoder in another branch, the first symmetric graph representation matrix is ​​obtained. First adjacency matrix and the second symmetric graph representation matrix Second adjacency matrix ; Fusion first feature representation matrix Second feature representation matrix The feature representation matrix is ​​obtained. ; Fusion First Symmetric Graph Representation Matrix Second symmetric graph representation matrix The symmetric graph representation matrix is ​​obtained. ; For the feature representation matrix and symmetric graph representation matrix The initial fusion representation matrix is ​​obtained by performing a weighted summation calculation. ; Based on normalized adjacency matrix Merging local structural information between data to enhance the initial fusion representation matrix To obtain the expression matrix ; where, the normalized adjacency matrix Indicating the intercellular relationship Normalized adjacency matrix of the nearest neighbor graph; An autocorrelation learning mechanism is used to utilize the nonlocal relationships between data in the fusion space based on initial information to obtain the representation matrix. The initial information fusion space is based on the perturbation attribute matrix. Build; For the expression matrix and expression matrix The fusion representation matrix is ​​obtained by performing a weighted summation calculation. ; For the fusion representation matrix Self-supervised clustering based on the optimal transmission strategy is performed to obtain the clustering results corresponding to the single-cell RNA sequencing data; Extract and visualize the fusion representation matrix The distribution of the clustering results is visualized.

Images