CN116386729A - scRNA-seq data dimension reduction method based on graph neural network - Google Patents

Abstract

The invention relates to data mining in bioinformatics, in particular to the mining of single-cell RNA sequencing data. Specifically, it involves dimensionality compression and clustering of single-cell RNA sequencing data through deep learning methods to achieve the purpose of effectively identifying cell populations. The method of the present invention includes collecting and preprocessing scRNA-seq data; constructing a graph neural network model; using the constructed model to reduce the dimension of the preprocessed data; and performing cluster analysis on the result after dimension reduction. Our model constrains the data structure and performs dimensionality reduction through a graph neural network module, and preserves both cell-cell relationships and gene-gene relationships in the dimensionality reduction results. Using standardized mutual information and adjusted Rand index as evaluation indicators, experiments on five real scRNA-seq datasets show that the method has good performance.

Description

A Dimensionality Reduction Method for scRNA-seq Data Based on Graph Neural Network

technical field

The invention relates to data mining in bioinformatics, in particular to the mining of single-cell RNA sequencing data. Specifically, it involves dimensionality compression and clustering of single-cell RNA sequencing data to achieve the purpose of effectively identifying cell populations.

Background technique

With the explosion of single-cell RNA sequencing (scRNAseq) technologies in recent years, unprecedented opportunities for single-cell transcriptional profiling have emerged. Traditional bulk RNA sequencing methods sequence mixtures of millions of cells. This results in the gene expression of a gene reflecting the average of gene expression in all cells, ignoring heterogeneity between cells. Unlike bulk RNAseq, scRNAseq isolates cells in the first step and sequences thousands of genes in each cell in the second step. Depending on the sequencing protocol, millions of expression values are collected for each gene, allowing the identification of new cell types, the determination of gene regulatory mechanisms, and the resolution of cellular dynamics during development.

Single-cell RNA sequencing (scRNA-seq) is an ideal method for studying cell-to-cell variation. Conventional dimensionality reduction techniques such as principal component analysis (PCA) and t-distributed stochastic neighborhood embedding (t-SNE) implemented on scRNA-seq data for visualization and downstream analysis significantly increase our knowledge of cellular heterogeneity and understanding of developmental progress. The recent emergence of massively parallel scRNA-seq (e.g., droplet platforms) enables the sequencing of millions of cells in complex biological systems, providing new insights into tissue and cellular microenvironment dissection, identification of rare/new cell types, developmental The inference of lineages and the elucidation of the mechanisms by which cells respond to stimuli offer excellent potential. However, the data generated by massively parallel scRNA-seq has the characteristics of high dropout, high noise, and complex structure, which brings a series of challenges to dimensionality reduction. In particular, preserving the complex topology between cells is a great challenge.

In the past few years, many dimensionality reduction methods for scRNA-seq data analysis have been developed or introduced. Recently developed competing methods include DCA, scVI, scDeepCluster, PHATE, SAUCIE, scGNN, ZINB-WaVE, and Ivis. Among them, deep learning shows the greatest potential. For example, DCA, scDeepCluster, Ivis, and SAUCIE adapted autoencoders to denoise, visualize, and cluster scRNA-seq data. However, these deep learning-based models only embed distinct cellular features while ignoring cell-to-cell relationships, which limits their ability to reveal complex intercellular topologies and also makes them difficult to elucidate developmental trajectories. The recently proposed graph autoencoder is very promising because it preserves long-distance relationships between data in the latent space.

However, studies have shown that gene interactions involved in gene regulatory networks or protein-protein interaction (PPI) networks are informative in different biological contexts. Furthermore, previous studies have shown that joint analysis of scRNA-seq data with prior gene interaction information can lead to meaningful understanding of the data. NetNMF-sc is a network regularized non-negative matrix factorization specially designed for scRNA-seq analysis, which utilizes prior gene networks to obtain more meaningful low-dimensional representations of genes. Correspondingly, scRNA-seq data also contain rich information to infer gene-gene interactions.

Inspired by the above understanding, we propose scTPGAE, a computational method based on graph neural networks, which utilizes two graph neural networks to simultaneously preserve cell-cell relationships, gene-gene relationships into dimensionality reduction results to achieve Better downstream analysis results.

Contents of the invention

Aiming at the problems existing in the above methods and the complexity of scRNA-seq data, the present invention proposes a dimensionality reduction method for scRNA-seq data based on a graph neural network. The method of the present invention can effectively solve the problems of important information loss and insufficient feature extraction existing in the existing dimensionality reduction methods, and simultaneously retain the cell-cell relationship and gene-gene relationship in the dimensionality reduction results, and obtain better clustering accuracy. The steps of the described method include:

1. Data preprocessing

First, assume we have a raw scRNA-seq count matrix C that filters out genes not counted in any cells. C can be expressed as a P by N dimensional matrix, where P is defined as the total number of genes, N is defined as the total number of cells, and _Cij represents the expression value of gene i in cell j.

In this work, we first preprocess the raw scRNA-seq count data, including log transformation and z-score normalization. We have a normalized output X, publicized as follows

X=zscore(X')

where _Sj is the size factor of each cell j. The advantage of data preprocessing is to preserve the effect of data size differences and convert discrete values to continuous values, thus providing greater flexibility for subsequent modeling.

In addition to the above-mentioned gene-cell relationship matrix, the input required by the graph neural network also requires a cell-cell relationship graph and a gene-gene interaction network.

Among them, the cell-cell relationship graph is constructed by the K nearest neighbor (KNN) algorithm in the Scikit-learn Python package. The default K is predefined as 35 in this study and adjusted according to the dataset in our experiments. The generated adjacency matrix is a 0-1 matrix, 1 means connected and 0 means disconnected.

The gene-gene interaction network can make use of existing data, and we collected seven different human gene interaction networks and one mouse gene interaction network to evaluate the performance of scTPGAE. One of the most famous gene interaction networks is the STRING database, a PPI network that collects and integrates protein-protein association information from various sources such as literature and experiments. HumanNet is a human functional gene network that integrates multiple types of omics data through a Bayesian statistical framework. HumanNet includes a hierarchical structure of human gene networks, i.e., human-derived PPIs, co-functional links, co-citations, and mutual exclusions from other species. Specifically, we use two versions of HumanNet, HumanNet-CF and HumanNet-PI, which consist of a co-functional network and a PPI network, respectively. FunCoup is a genome-wide functional association network that uses a unique redundancy-weighted Bayesian integral to combine 10 different types of functional association data. GeneMANIA creates combinatorial gene networks by weighting multiple functional genomic datasets. In addition, we collected two functional similarity matrices from pgWalk, which are from KEGG pathways and Gene Ontology biological processes, respectively. Next, we converted these two similarity matrices into a gene network by filtering out those gene pairs whose similarity value is less than a certain threshold (i.e., 0.9). These two networks are called pgWalk-kegg and pgWalk-gobp respectively.

2. Building a graph neural network for dimensionality reduction

(1) Graph neural network G1 that preserves cell-cell relationships

Graph autoencoders are a type of artificial neural network for unsupervised representation learning on graph-structured data. Graph autoencoders have low-dimensional bottleneck layers and thus can be used as dimensionality reduction models. Suppose the input is a cell-cell graph of node matrix X and adjacency matrix A. In our joint graph autoencoder, there is one encoder E for the whole graph and two decoders D _X and D _A for nodes and edges respectively. In practice, we first encode the input graph as a latent variable h=E(X,a), and then decode h into a reconstructed node matrix x _r =D _X (h) and a reconstructed adjacency matrix A _r =D _A (h). The goal of the learning process is to minimize the reconstruction loss

where weights are hyperparameters. In our experiments, it is set to 0.6.

We implement our model using the Python package Spektral32. There are many types of graph neural networks that can be used as encoders or decoders. Therefore, in order to extract features of a node with the help of its neighbors, we apply a graph attention layer as default in the encoder. Other graph neural networks such as GCN, GraphSAGE, and TAGCN can also be implemented as encoders in scTPGAE. The feature decoder _DX is a four-layer fully connected neural network with 64, 256, 512 nodes in the hidden layers.

The edge decoder consists of a fully connected layer followed by quadrantization and activations:

A _r ＝D _A (h)＝σ(ZZ ^T )

where Z = σ(Wh) as the output of a fully connected layer with a weight matrix W, and σ(x) = max(0,x) is the straight linear unit.

(2) Graph neural network G2 that retains gene-gene relationships

We noticed that when applying the gene interaction network to a certain dataset, only those interaction pairs where two interacting genes appear in the dataset are kept, and the rest are discarded. In other words, the number of interaction pairs of gene interaction networks of different datasets may be different from each other. To capture the two regulatory directions and their corresponding strengths in a pair of genes, the gene interaction network is considered as a directed graph, so for the edges of gene A and gene B from an undirected gene network, such as the STRING PPI network, we Think of it as a pair of edges (i.e. the edge from A to B and the edge from B to A).

The specific construction method of the graph neural network is the same as that of the graph neural network that preserves the cell-cell relationship, except that the input of the graph neural network is converted from a cell-cell relationship graph to a gene-gene relationship PPI interaction network. The interaction relationship between genes can be presented spontaneously in a graph format, where a graph neural network is applied to model this relationship. In the graph convolution layer, each node represents a gene, and the edge between two nodes represents the relationship between these two corresponding genes. The graph representation module is designed as a graph convolutional layer, which updates each node by aggregating the information of its neighbor nodes.

3. Dimensionality reduction of scRNA-seq data

Dimensionality reduction of the preprocessed scRNA-seq data was performed using the constructed graph neural network.

The gene-cell count matrix and cell-cell relationship are input into the graph neural network G1 to obtain the dimensionality-reduced cell feature θ1.

The gene-cell count matrix and the gene-gene interaction network are input into the graph neural network G2 to obtain the dimensionality-reduced cell feature θ2.

The learned cell features are concatenated as dimensionality reduction results for subsequent downstream analyses.

4. K-means algorithm clustering

This method uses the ZINB conditional likelihood to reconstruct the decoder output of scRNA-seq data. The ZINB distribution is proved to be a good model for scRNA-seq data and is a generally accepted gene expression distribution structure. .

In order to evaluate the effectiveness of the method, we applied the k-means clustering algorithm to cluster the data after dimensionality reduction, and evaluated it with the index of normalized mutual information. Assuming that X is the predicted clustering result and Y is the cell type with the true label, the NMI score is calculated as follows:

MI is the mutual entropy between X and Y, and H is the Shannon entropy.

It can be seen from the above that the dimensionality reduction method for scRNA-seq data based on the graph neural network provided by one or more embodiments of this specification preserves both cell-cell relationships and gene-gene relationships in the dimensionality reduction results. Our model constrains the data structure and performs dimensionality reduction through two graph neural network modules. Experiments on five real scRNA-seq datasets demonstrate that our method can provide more accurate low-dimensional representations of scRNA-seq data.

Detailed ways

In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be further described in detail in combination with experiments below. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention.

1. Dataset overview

To evaluate the performance of scTPGAE, we focus on relatively large datasets; five real scRNA-seq datasets with known cell types are selected. The table below summarizes the basic information of the five real datasets, which we describe below.

(i) 10X PBMC dataset, provided by the 10X scRNA-seq platform, collected from a single healthy human; (ii) mouse embryonic stem cell dataset describing the heterogeneous differentiation of mouse embryonic stem cells after withdrawal of leukemia inhibitory factor (LIF) (iii) The mouse bladder cell dataset is from the Mouse Cell Atlas project GSE108097. From the raw count matrix, we selected approximately 2700 cells from bladder tissue; (iv) the worm neuronal cell dataset was analyzed by single-cell combinatorial indexed RNA-sequencing from L2 larval stage C. elegans; (v ) Zeisel dataset contains 3005 cells from mouse cortex and hippocampus GSE60361.

2. Experimental environment and parameter settings

The hardware environment is mainly a PC host. Among them, the CPU of the PC host is 11th Gen Intel(R) Core(TM) i5-1135G7, 2.42GHz, the memory is 16GB RAM, and the 64-bit operating system. The software uses Windows 10 as the platform and is implemented in the Python language in the Pycharm environment. The python version is 3.5.0 and the Tensorflow version is 1.4.0.

We implement our model using the Python package Spektral32. There are many types of graph neural networks that can be used as encoders or decoders. Therefore, in order to extract features of a node with the help of its neighbors, we apply a graph attention layer in the encoder as default. Other graph neural networks such as GCN, GraphSAGE, and TAGCN can also be implemented as encoders in scTPGAE. The feature decoder _DX is a four-layer fully connected neural network with 64, 256, 512 nodes in the hidden layers.

The edge decoder consists of a fully connected layer followed by quadrantization and activations:

A _r ＝D _A (h)＝σ(ZZ ^T )

where Z = σ(Wh) as the output of a fully connected layer with a weight matrix W, and σ(x) = max(0,x) is the straight linear unit.

In addition to the above-mentioned gene-cell relationship matrix, the input required by the graph neural network also requires a cell-cell relationship graph and a gene-gene interaction network.

Among them, the cell-cell relationship graph is constructed by the K nearest neighbor (KNN) algorithm in the Scikit-learn Python package. The default K is predefined as 35 in this study and adjusted according to the dataset in our experiments. The generated adjacency matrix is a 0-1 matrix, 1 means connected and 0 means disconnected.

The gene-gene interaction network can make use of existing data, and we collected seven different human gene interaction networks and one mouse gene interaction network to evaluate the performance of scTPGAE.

3. Evaluation indicators

To make the results of different methods easy to compare, we employ K-means for cluster analysis and set the parameter K to the true number of clusters in each dataset. In our experiments, scTPGAE models are evaluated using two metrics, Normalized Mutual Information (NMI) and Adjusted Rand Index (ARI), which are widely used in model performance evaluation in unsupervised learning scenarios.

4. Analysis of experimental results

Here, this method is mainly tested on five real data sets, and the normalized mutual information and adjusted Rand index obtained are shown in the table below.

Normalized Mutual Information (NMI)

Adjusted Rand Index (ARI)

The above experimental results show that the scTPGAE method based on graph neural network is a promising new method. This method achieved good performance on five real datasets, which indicates that this method can provide a more accurate low-dimensional representation of scRNA-seq data.

It can be seen that the scTPGAE method we proposed is a method for dimensionality reduction and cluster analysis of single-cell RNA-seq data. This method has the following advantages. First, scTPGAE combines the potential space distribution with the selected prior Matching; secondly, scTPGAE retains the relationship between cells and cells in the dimensionality reduction results; thirdly, the scTPGAE method preserves the relationship between genes and genes while preserving the cell-cell relationship; finally, this method considers Parallel and scalable features in deep neural network frameworks. Our model constrains the data structure and performs dimensionality reduction through a graph neural network module. Using standardized mutual information and adjusted Rand index as evaluation indicators, experiments on five real scRNA-seq datasets show that the method has good performance.

Description of drawings

Figure 1: Schematic flow chart of scRNA-seq data dimensionality reduction method based on graph neural network;

Figure 2: Experimental results using normalized mutual information (NMI) as a measure;

Figure 3: Experimental results using the Adjusted Rand Index (ARI) as a measure.

Claims (5)

1. A scRNA-seq data dimension reduction method based on a graph neural network is characterized by comprising the following implementation steps:

(1) Preprocessing data; collecting scRNA-seq datasets from different species, different types, different cell numbers; preprocessing the collected original scRNA-seq data by adopting a logarithmic conversion and z fraction normalization method, and reconstructing the input data by utilizing zero expansion negative binomial distribution to obtain noiseless data;

(2) Constructing a graphic neural network for dimension reduction, which is an automatic encoder framework consisting of a depth encoder, an intermediate hidden layer and a depth decoder; the topological structure between cells and the topological structure between genes can be simultaneously reserved in the dimension reduction result;

(3) Reducing the dimension of the preprocessed scRNA-seq data by using the constructed graph neural network, learning a hidden layer feature vector by using an intermediate hidden layer of an automatic encoder, restraining prior distribution of the hidden layer feature vector, and matching the hidden layer feature vector with the selected prior distribution; connecting the hidden layer feature vectors learned in the two graph neural networks so as to facilitate subsequent downstream analysis;

(4) And clustering the dimensionality reduced data by using a k-means clustering algorithm to obtain a standardized mutual information score and adjust the Rand index.

2. The method for reducing dimension of scRNA-seq data based on graphic neural network according to claim 1, wherein the data is collected and the collected single cell RNA sequencing data is preprocessed:

we collected five scRNA-seq datasets from different species, different types, different cell numbers, and were then preprocessed using the method of logarithmic transformation and z-score normalization.

Specifically, we performed data preprocessing operations on the following five data sets.

(1) 10X PBMC dataset, provided by the 10X scRNA-seq platform, data collected from a healthy human;

(2) A mouse embryonic stem cell dataset describing a transcriptome of mouse embryonic stem cell heterodifferentiation following withdrawal of Leukemia Inhibitory Factor (LIF);

(3) The mouse bladder cell dataset was from the mouse cytogram project GSE108097. From the original count matrix, we selected about 2700 cells from bladder tissue;

(4) The worm neuron cell dataset was analyzed by single cell combinatorial indexing RNA sequencing from L2 larval stage caenorhabditis elegans;

(5) The Zeisel dataset contained 3005 cells from mouse cortex and hippocampal GSE60361.

3. The method for reducing dimension of scRNA-seq data based on graphic neural network according to claim 1, wherein the construction of a graphic neural network is an automatic encoder framework composed of a depth encoder, an intermediate hidden layer and a depth decoder, and specifically comprises:

(1) Graphic neural network G1 retaining cell-cell relationship

The graph automatic encoder is an artificial neural network for unsupervised representation learning of graph structure data. The graphic auto-encoder has a low-dimensional bottleneck layer and thus can be used as a dimension-reduction model. Assume that the inputs are a cell-cell relationship graph of node matrix X and adjacency matrix a. In our joint picture automatic encoder, there is one encoder E for the whole picture, two decoders D _X And D _A For nodes and edges, respectively. In practice, we first encode the input graph as the latent variable h=e (X, a), and then decode h into the reconstructed node matrix X _r ＝D _X (h) And a reconstructed adjacency matrix A _r ＝D _A (h) A. The invention relates to a method for producing a fibre-reinforced plastic composite The goal of the learning process is to minimize reconstruction losses

Wherein the weights are superparameters. In our experiments, set to 0.6.

We use Python package Spektral32 to implement our model. There are many types of graphic neural networks that can be used as encoders or decoders. Therefore, to extract the features of the nodes by means of their neighbors, we apply the graph attention layer as default in the encoder. Other graph neural networks such as GCN, graphSAGE and TAGCN may also be implemented as encoders in the scTPGAE. Feature decoder D _X Is a four-layer fully connected neural network with 64, 256 and 512 nodes in the hidden layer.

The edge decoder consists of one fully connected layer, then the components of quadrant and activation:

A _r ＝D _A (h)＝σ(ZZ ^T )

where z=σ (Wh) as the output of the fully connected layer with the weight matrix W, σ (x) =max (0, x) is a straight linear unit.

(2) Graph neural network G2 retaining gene-gene relationship

We note that when a gene interaction network is applied to a data set, only those interaction pairs in which two interacting genes occur in the data set are retained, and the remaining pairs are discarded. In other words, the number of interaction pairs of the gene interaction networks of different data sets may differ from each other. To capture both regulatory directions and their corresponding intensities in a pair of genes, the gene interaction network is considered a directed graph, so for the edges of the a and B genes from the undirected gene network, e.g., STRING PPI network, we consider it as a pair of edges (i.e., the edge from a to B and the edge from B to a).

The specific graph neural network construction method is the same as that of the graph neural network which retains the cell-cell relationship, except that the input of the graph neural network is converted from the cell-cell relationship graph into the PPI interaction network of the gene-gene relationship. The interaction relationships between genes can spontaneously be presented in a graphical format, where a graphical neural network is applied to model such relationships. In the graph roll stack, each node represents one gene, and the edge between two nodes represents the relationship of the two corresponding genes. The graph representation module is designed as a graph volume layer, updating each node by aggregating the information of its neighboring nodes.

4. The method for reducing the dimension of scRNA-seq data based on the graphic neural network according to claim 1, wherein the method for reducing the dimension of the preprocessed scRNA-seq data by using the constructed graphic neural network is characterized by comprising the following steps:

and (3) performing dimension reduction on the preprocessed scRNA-seq data by using the constructed graph neural network.

Inputting the gene-cell count matrix and the cell-cell relation into the graph neural network G1 to obtain the cell characteristics theta 1 after dimension reduction.

Inputting the gene-cell count matrix and the gene-gene interaction network into the graph neural network G2 to obtain the cell characteristics theta 2 after dimension reduction.

The learned cell characteristics are linked as a dimension reduction result of subsequent downstream analysis.

5. The method for reducing the dimension of scRNA-seq data based on the graphic neural network according to claim 1, wherein the k-means clustering algorithm is applied to cluster the dimension-reduced data. The method specifically comprises the following steps:

the present method uses the ZINB conditional likelihood to reconstruct the decoder output of the scRNA-seq data, and the ZINB distribution has proven to be a better model for describing the scRNA-seq data and is a widely accepted gene expression distribution structure.

In order to evaluate the effectiveness of the method, a k-means clustering algorithm is applied to cluster the data after dimension reduction, and standardized mutual information and an adjusted Rand index are used as evaluation indexes. Experiments performed on five real scRNA-seq datasets indicate that the present method can provide a more accurate low-dimensional representation of the scRNA-seq data.

Images