A multimodal deep learning model to infer cell-type-specific functional gene networks

<jats:title>Abstract</jats:title><jats:sec>                 <jats:title>Background</jats:title>                 <jats:p>Functional gene networks (FGNs) capture functional relationships among genes that vary across tissues and cell types. Construction of cell-type-specific FGNs enables the understanding of cell-type-specific functional gene relationships and insights into genetic mechanisms of human diseases in disease-relevant cell types. However, most existing FGNs were developed without consideration of specific cell types within tissues.</jats:p>               </jats:sec><jats:sec>                 <jats:title>Results</jats:title>                 <jats:p>In this study, we created a multimodal deep learning model (MDLCN) to predict cell-type-specific FGNs in the human brain by integrating single-nuclei gene expression data with global protein interaction networks. We systematically evaluated the prediction performance of the MDLCN and showed its superior performance compared to two baseline models (boosting tree and convolutional neural network). Based on the predicted cell-type-specific FGNs, we observed that cell-type marker genes had a higher level of hubness than non-marker genes in their corresponding cell type. Furthermore, we showed that risk genes underlying autism and Alzheimer’s disease were more strongly connected in disease-relevant cell types, supporting the cellular context of predicted cell-type-specific FGNs.</jats:p>               </jats:sec><jats:sec>                 <jats:title>Conclusions</jats:title>                 <jats:p>Our study proposes a powerful deep learning approach (MDLCN) to predict FGNs underlying a diverse set of cell types in human brain. The MDLCN model enhances prediction accuracy of cell-type-specific FGNs compared to single modality convolutional neural network (CNN) and boosting tree models, as shown by higher areas under both receiver operating characteristic (ROC) and precision-recall curves for different levels of independent test datasets. The predicted FGNs also show evidence for the cellular context and distinct topological features (i.e. higher hubness and topological score) of cell-type marker genes. Moreover, we observed stronger modularity among disease-associated risk genes in FGNs of disease-relevant cell types. For example, the strength of connectivity among autism risk genes was stronger in neurons, but risk genes underlying Alzheimer’s disease were more connected in microglia.</jats:p>               </jats:sec>

section_content: Functional gene networks (FGNs) capture functional relationships among genes and provide a systemlevel understanding of gene function, that could further shed light on genetic mechanisms underlying human diseases. Each node in a gene network represents a gene and each edge represents a functional connection between a gene pair. Evidence of functional connections include gene co-expression, physical interaction of proteins, and text mining. Gene networks vary across tissues and cell types because some genes are functionally related only in certain tissues and cell types due to the speci city needed to give rise to an array of functions. However, most curated gene networks that have been developed, such as STRING [1] and HumanNet [2], incorporate evidence from many tissues, cell types, and organisms, and may not closely re ect the gene relationships within speci c tissues or cell types. Our understanding of gene functions and genetic mechanisms underlying human diseases could then be limited in diseaserelevant cell types. 

The recent advances in deep learning techniques and emerging single cell datasets have opened up new opportunities for in silico prediction of cell-type-speci c FGNs. A convolutional neural network model (CNNC) was recently developed to predict gene relationships in speci c cell types through a novel encoding scheme that transforms pairwise gene expression data to an image-like object [3]. The CNNC model provides a exible computational framework and opportunity for further improvement by incorporating domain knowledge for a particular task. In this study, we are interested in predicting cell type-speci c gene functional relationships in the human brain. We reason that integrating CNNC framework with gene features derived from a global protein interaction network could learn a more comprehensive representation of genes and enable more accurate predictions of gene functional relationships. Accordingly, we developed a multimodal deep learning model (MDLCN) for predicting celltype-speci c FGNs by integrating single-nuclei expression data of the human brain with a global protein interaction network. We systematically evaluated the prediction performance of MDLCN and showed it outperformed the CNNC and a conventional machine learning model (i. e. boosting tree model) that transforms the pairwise gene expression data to a correlation coe cient for predicting cell-type-speci c FGNs. We further evaluated the cellular context of the predicted cell-type-speci c FGNs through topological analyses of cell-type marker genes and genes underlying autism and Alzheimer's disease.
section_name: Background

section_content: We developed MDLCN, a multimodal deep learning model, for predicting cell-type-speci c FGNs by leveraging single-cell gene expression data with a global protein interaction network (Fig. 1 ). Gene expression signatures of a gene pair were rst transformed to a co-expression matrix that captures the joint density of co-expression patterns of the gene pair across the cells in a particular cell type. We computed a set of proximity features for each gene pair based on a global protein interaction network we assembled from protein physical interaction evidence. The co-expression matrix and global proximity features of each gene pair were integrated to predict its functional relationship status through the MDLCN model. The MDLCN model was trained for each cell type using the cell-type-speci c gold standards.
section_name: Methodology

section_content: We assembled a gold standard of gene pairs for building cell-type-speci c FGNs following the approach for building tissue-speci c FGNs in a previous study [4]. We rst constructed a cell-type-naive functional relationship gold standard from 564 expert-selected gene ontology (GO) terms and experimentally derived gene annotations. Gene pairs co-annotated to expert-selected terms were treated as positive examples (i. e. functionally related), and pairs not co-annotated to any of these terms were considered as negative examples. Next, we identi ed cell-type-speci c genes that are de ned as the top ranked genes by a speci city score that was computed by the average expression level of the gene in each cell type divided by the total expression values of the gene across all cell types [5]. The higher the speci city score of a gene in a cell type, the more speci c the gene was to that cell type. Then, we combined the cell-type-naive gold standard with the cell-type-speci c genes to construct the cell-type-speci c gold standard. 

Ultimately, our cell-type-naive gold standard included 3,619,063 positive gene pairs and 49,095,410 negative gene pairs. Table. S1 presents the positive and negative class construction procedure (see Additional le 1). The positive gene pairs in our cell-type-speci c gold standard were a subset of positive examples in the cell-type-naive gold standard, and the two genes of each pair were cell type-speci c or one was cell type-speci c and the other was a house-keeping gene as de ned previously [6]. The negative gene pairs in our cell-type-speci c gold standard were either: 1) positive examples in cell-type-naive gold standard, but one gene was speci c to the corresponding cell type, and the other gene was speci c to a different cell type; 2) positive examples in cell-type-naive gold standard, but one gene was speci c to a different cell type, and the other gene was a house-keeping gene; 3) negative examples in cell-type-naive gold standard, but the two genes of each pair are cell type-speci c or one was cell type-speci c and the other was a house-keeping gene; 4) negative examples in cell-type-naive gold standard, but one gene was speci c to the corresponding cell type, and the other gene was speci c to a different cell type; 5) negative examples in cell-type-naive gold standard, but one gene was speci c to a different cell type, and the other gene was a house-keeping gene. 

The single-nuclei gene expression data used in this study includes 14,873 nuclei from the human brain that were clustered and annotated to major brain cell types in a previous study [7], including 3,400 excitatory neurons, 1,715 inhibitory neurons, 1,897 astrocytes, 245 endothelial cells, 386 microglia, 2,963 oligodendrocytes, and 682 oligodendrocyte precursor cells. We normalized the single nuclei count data by the "LogNormalize" function of Seurat that normalizes each feature count by the total counts in each cell, multiplied by a scale factor (10,000) and transformed to log scale [8]. Cell type-speci c genes were the ones which ranked in the top 5% in the speci city scores for all cell types except for excitatory neurons for which the top 10% ranked genes were used so that a su cient number of labeled gene pairs for model training could be collected. To have a balanced number of positive and negative gene pairs in each training set, the negative classes were randomly down sampled [9]. The number of positive and negative gene pairs in each cell type is presented in Table. S2 (see Additional le 1).
section_name: Training dataset

section_content: We encoded single-nuclei gene expression data of the gene pairs to an image-like 2D co-expression matrix. Speci cally, we rst imputed and smoothed single nuclei expression data from the human brain described above using MAGIC, a Markov a nity-based graph imputation method [10]. Then, the range of the expression values of each gene was divided into equal bins. Then, for each pair of genes, 2D coexpression matrix was constructed by counting the number of cells that each gene expressed in the corresponding bins [3]. As the number of bins ( ) plays an important role in the model performance in our experiments, it was tuned to 10, at which the model achieved the best prediction accuracy.
section_name: Transforming single-nuclei gene expression data to 2D coexpression matrices

section_content: Global protein interaction networks contain protein structural and functional information that are informative for cell-type-speci c gene functional relationships. We assembled a global protein interaction network based on experimentally validated protein physical interaction evidence from multiple resources including Biogrid [11], IntAct [12, 13], APID [14] and Inweb [15]. 

After overlapping with genes in the single-nuclei gene expression dataset, the global protein interaction network contained 16,873 genes and 142,340,628 pairs of physical interactions. We used ve metrics to measure the degree of similarity or proximity between a protein pair in the global protein interaction network, including Common Neighbors (CN), Jaccard's Coe cient (JC), Preferential Attachment (PA), Adamic-Adar Coe cient (AA), and Path Distance (PD). These metrics measure different topological relationships between two proteins in the network. In particular, CN counts the number of common neighborhoods between the two proteins, JC quanti es the similarity between their neighborhoods, PA calculates the likelihood of link existence by measuring the strength of the hubness of the two proteins, AA computes the proportion of their shared links to the total number of their neighbors, and PD measures the length of the shortest path between the two proteins [16].
section_name: Gene proximity features from global protein interaction network

section_content: We developed a multimodal deep learning model to predict cell-type-speci c gene functional relationships from co-expression matrices and proximity features between two proteins in the global protein interaction network. For each pair of genes, the co-expression matrix and the vector of proximity features were exploited as two modalities in our model, including a co-expression-processor modality to extract representations from the co-expression matrix and a proximity-processor modality to extract representations from proximity features as shown in Fig. 2. In the co-expression-processor modality, the input layer is a co-expression matrix for each gene-pair. The modality consists of three convolutional layers which map the local conjunctions of features from previous layers to a feature map. Immediately after each convolutional layer, there is a max pooling layer which down samples the output of convolutional layers by taking the maximum value over an input window. At the end of this modality, a attened layer is used to switch 2D features extracted from convolutional process to 1D features by retaining the weight orders and a densely connected layer is employed to compile the features extracted from previous layers to form the representations. The proximity-processor modality consists of an input layer for ve proximity features between each gene pair, four densely connected layers and a attened layer. The representations output from the two modalities are concatenated to a high-dimensional feature vector in a fusion layer and transformed through three densely connected layers. Finally, the feature vector is used in output layer to predict the probability of cell-type-speci c functional connection between a gene pair. We used recti ed liner activation function (ReLU) as the activation function across the whole network except the output layer where sigmoid function was used for binary classi cation. The dropout regularization is used in the multimodal deep learning model for preventing over tting. The model was implemented using the Keras library in Python. We chose the Binary-Cross Entropy as the loss function and the Adam optimizer to update weights. The hyperparameters in the model, including the number of lters in the convolutional layer, the kernel size of the convolutional layer, the kernel size of the max pooling, the size of the dense layer in co-expression modality, the size of the dense layers in proximity modality, the position and rate of dropout, and the optimizer's type, were tuned using the validation set and summarized in Table. S3 in Additional le 1.
section_name: Predicting cell-type-speci c FGNs

section_content: We evaluated the performance of the MDLCN model in predicting functional relationships of gene pairs in three different test sets, including a 1) dependent test set that included gene pairs with both genes appearing in the training set; 2) partially dependent test set that included gene pairs with only one gene appearing in the training set; and 3) independent test set that included gene pairs with both genes not appearing in the training set. We used the area under ROC curve (AUC-ROC) and the area under Precision-Recall curve (AUC-PRC) as the evaluation metrics with ve-fold cross validation. To evaluate the effects of the features from global protein interaction network, we compared the MDLCN with the CNNC that only uses 2D co-expression matrix to predict the interactions among the genes. To determine the strength of the 2D co-expression matrix in predicting the interaction among the genes, we also compared the MDLCN with a boosting tree model that uses the Pearson correlation to measure the co-expression between two genes.
section_name: Evaluation of model performance

section_content: To evaluate the constructed cell-type-speci c networks, we further examined whether cell-type marker genes have a distinctive topological structure across different cell type-speci c FGNs. Cell-type marker genes were based on a previous study [17] and were de ned as genes with at least one log-fold change in expression levels when cells of a given cell-type were compared against all other cells. The hubness of each marker gene was computed as the summation of edge weights that are directly connected to that gene. We also computed a topological speci city score for each marker gene to test whether cell-type marker genes have distinctive localization compared to random networks as done previously [18]. The topological speci city score represents the hubness of a gene in the predicted network normalized by its hubness distribution in random networks that were created by re-shu ing edge weights in the predicted network. The topological speci city score (topS) for each marker gene is calculated as
section_name: Topological analysis of cell-type marker genes

section_content: To evaluate whether disease genes show cell-type-speci c modularity in the constructed cell-type-speci c networks, we assessed the connectivity strength between disease genes in each network. We considered two brain disorders: autism spectrum disorder (ASD) and Alzheimer's disease (AD). We collected 408 high con dence ASD risk genes from the SFARI database [19] and 1,611 genes implicated in AD from the DisGeNET database [20]. We calculated the average connectivity over all pairs of disease genes for each network and compared the average connectivity with a background distribution from 1,000 random gene sets matched by gene numbers and gene length with the disease genes. We used the Z-score to test the deviation of connectivity of disease genes
section_name: Connectivity of disease genes in predicted cell-typespeci c FGNs

section_name: Results

section_content: We rst compared the 2D co-expression matrices between the positive and negative gene pairs for each cell type (Fig. 3 ). In all cell types, the positive gene-pairs show on average higher values in the bins corresponding to higher co-expression levels than the negative pairs, suggesting that the 2D coexpression matrices capture the complicated correlation among the genes to distinguish positive and negative classes. The proximity features extracted from global protein interaction network also show signi cant difference between the positive and negative gene pairs in all cell types, with the positive gene pairs having higher values for CN, JC, PA, and AA scores and having lower values in the PD score Prediction performance of the MDLCN model is illustrated in Fig. 4. The MDLCN model showed higher prediction accuracy in dependent, partially dependent, and independent test sets in all cell types compared to the CNNC and boosting tree models. For example, the AUC-ROC achieved by the MDLCN model is around 15% higher than the boosting tree model and 4% higher than the CNNC model in all three different testing scenarios. The improvements on AUC-PRC achieved by the MDLCN model are more substantial, suggesting the advantage of the model in identifying true cell-type-speci c gene functional relationships. 

For all the cell types, the prediction accuracy of the MDLCN model is satisfactory. The AUC-ROC (AUC-PRC) was higher than 0. 90 (0. 72) for the dependent test set, 0. 87 (0. 67) for the partially dependent test set, and 0. 87 (0. 63) for the independent test set (Fig. 4 ). independent and partially dependent test sets consist of functional relationships among new genes, which were not seen in the training data set and more challenging to predict. Nonetheless, the MDLCN model still achieved good prediction performance in these test sets, which indicates that our model is promising in its ability to predict functional relationships between new genes.
section_name: MDLCN model performance

section_content: We evaluated the topological features of cell-type makers genes in each constructed cell-type-speci c FGN. We observed that marker genes had higher hubness than non-marker genes in the cell type that corresponds to the marker genes (Fig. 5 ). Furthermore, marker genes had higher topological speci city score in their corresponding cell type than the rest of the cell types (Fig. 6 ). These observations demonstrated that constructed cell-type-speci c FGNs show distinct topological features for cell-type marker genes, re ecting the cellular context of predicted networks.
section_name: Downstream analysis Topological analysis of cell-type maker genes

section_content: We compared the strength of connectivity among risk genes that underlie ASD and AD across different cell types (Fig. 7 ). We observed that disease genes tended to be functionally related in each cell-typespeci c network, but the strength of connectivity among disease genes varied across different cell types. 

For example, the strength of connectivity among ASD risk genes is stronger in astrocytes, and neurons, and lower in microglia, endothelial cells, and oligodendrocytes. On the other hand, risk genes underlying AD were more connected in microglia, astrocytes, and endothelial cells, but less connected in neurons and oligodendrocytes. These results were consistent with the literature as ASD etiology is more related to the dysfunction of neurons and astrocytes [21], while microglia cells play a key role in AD pathogenesis [22], providing further evidence for the cellular context of predicted cell type-speci c FGNs. The hubness of marker genes versus non-marker genes in seven cell types: endothelial cells (Endo), excitatory neurons (Ex), oligodendrocyte precursor cells (Opc), astrocytes (Ast), inhibitory neurons (In), microglia (Mg), and oligodendrocytes (Oli). The P-values for even cell types are between 10 -41 and 10 -5 Figure 6 Topological speci city scores (topS) of marker genes in seven cell types: endothelial cells (Endo), excitatory neurons (Ex), oligodendrocyte precursor cells (Opc), astrocytes (Ast), inhibitory neurons (In), microglia (Mg), and oligodendrocytes (Oli). The P-values for seven cell types are between 10 -10 and 10 -4.
section_name: Connectivity of disease genes in predicted cell type-speci c network

section_content: Figure 7
section_name: Discussion And Conclusions

section_name: Acknowledgements

section_content: Not Applicable.
section_name: Acknowledgements

section_name: Funding

section_content: This research was funded by NIMH, R01MH121394.
section_name: Funding

section_name: Availability of data and materials

section_content: The protein physical interaction evidence was downloaded from multiple resources including Biogrid (https://wiki. thebiogrid. org/doku. php/biogrid_tab_version_2. 0), IntAct (https://ftp. ebi. ac. uk/pub/databases/intact/current/psimitab), and APID (http://cicblade. dep. usal. es:8080/APID/init. action#tabr1). The single-nuclei gene expression data and the list of house-keeping genes are available at http://www. gtexportal. org/home/ datasets and https://m. tau. ac. il/~elieis/HKG/HK_genes. txt, respectively. The curated training data and cell-typespeci c FGNs predicted from this study are available from the corresponding author on reasonable request.
section_name: Availability of data and materials

section_name: Declarations Supplementary materials

section_content: We proposed a multimodal deep learning model MDLCN to predict cell type-speci c FGNs by integrating single-cell gene expression data with global protein interaction networks. We showed the superior performance of the MDLCN model compared to both the CNNC and boosting tree models. We further demonstrated evidence for the cellular context of predicted cell type-speci c FGNs through the distinct topological features of cell-type marker genes and risk genes underlying two brain disorders: autism and Alzheimer's disease. 

While the MDLCN model holds the promise to predict functional gene relationships underlying a diverse set of cell types, the model should be viewed in light of two limitations. First, model performance may depend on the quality of the training dataset since the MDLCN employs a supervised approach. Second, the model can only predict the functional associations of two genes but not the direction of association. Further work includes extending the prediction to directional relationship of two genes by integrating more functional genomic datasets, such as those from ChIP-Seq or gene perturbation experiments.

section_content: Supplementary Tables and Figures.
section_name: Declarations Supplementary materials

section_content: The dataset was used in study belongs to public databases. Users can download relevant data for free for research and publish relevant articles. Our study is based on open access data, so there are no ethical issues and other con icts of interest. All methods were carried out in accordance with relevant guidelines and regulations.
section_name: Ethics approval and consent to participate

section_content: Not applicable.
section_name: Consent for publication

section_content: The authors declare that they have no competing interests. Z-scores for AD and ASD risk genes in endothelial cells (Endo), excitatory neurons (Ex), oligodendrocyte precursor cells (Opc), astrocytes (Ast), inhibitory neurons (In), microglia (Mg), and oligodendrocytes (Oli).
section_name: Competing interests

section_content: This is a list of supplementary les associated with this preprint. Click to download. Additional le1. pdf
section_name: Supplementary Files