Using single cell atlas data to reconstruct regulatory networks

Abstract

Inference of global gene regulatory networks from omics data is a long-term goal of systems biology. Most methods developed for inferring transcription factor (TF)-gene interactions either relied on a small dataset or used snapshot data which is not suitable for inferring a process that is inherently temporal. Here, we developed a new computational method that combines neural networks and multi-task learning to predict RNA velocity rather than gene expression values. This allows our method to overcome many of the problems faced by prior methods leading to more accurate and more comprehensive set of identified regulatory interactions. Application of our method to atlas scale single cell data from 6 HuBMAP tissues led to several validated and novel predictions and greatly improved on prior methods proposed for this task.

Figure 1.

The flowchart of the MTLRank framework. (A) Hubmap tissue-specific scRNA-seq and scATAC data sets were used. TF RPKM matrix was generated from tissue-specific scRNA-seq data and TF activity matrix was generated from the integration of Cistrome DB ChIP-seq TF prior information and scATAC-seq data. (B) Computation of TF activity scores from ChIP-seq data and scATAC-seq data. Blue blocks represent TF binding sites from ChIP-seq data and green block represents open chromatin region from scATAC-seq data. TF binding sites were weighted by the scATAC-seq open chromatin regions (see Materials and Methods for more details). (C) Multi-task based model training and deep SHAP based TF ranking. The final tissue-specific regulatory networks were constructed from the TF ranking results.

Figure 2.

Comprehensive evaluation of MTLRank models. (A) R² scores for the 500 randomly sampled genes in each tissue. TFA: TF activity score matrix, velo: RNA velocities; NN: baseline neural network model. (B) R² scores for the 500 randomly sampled genes in each tissue from SNARE-seq data. Left part represents R² for all 500 genes and right part represents R² for genes that have >200 available TF activity scores from the SNARE-seq data. The bars marked with ‘single cell’ stand for the scATAC-seq signals that were paired with the scRNA-seq data at single cell level and the bars marked with ‘average’ stand for the scATAC-seq signals that were averaged out across the cells. (C) Recall percentage of the TF-Marker database marker genes. The recall percentage was computed for the TFs in each tissue-specific network. The two numbers on top of each bar indicates the total number of recallable TFs versus recalled TFs by the method. P-values from hypergeometric test were also marked on top the bars. P-values <0.05 were marked by red color. Total number of recallable TFs are the intersection between the tissue markers from TF-Marker database and all available TFs in RPKM matrix.

Figure 3.

Tissue-specific subnetwork with distinctive organ functions. (A) Predicted drug metabolism network in liver. The top 10 regulators and top 10 target genes from the original tissue-specific network were marked by red circles. (B) Predicted MAPK-Th1/Th2 network in spleen. The top 10 regulators and top 10 target genes from the original tissue-specific network were marked by red circles. (C) Predicted calcium signaling network in kidney. The top 10 regulators and top 10 target genes from the original tissue-specific network were marked by red circles.