Caleb C ReagorpersName, Nicolas Velez-AngelpersName, A J HudspethpersName (2023date). “Depicting pseudotime-lagged causality across single-cell trajectories for accurate gene-regulatory inference.” PNAS Nexustitle·j.
Object of study — inferring direct gene-regulatory interactions from pseudotime-lagged single-cell trajectories with a convolutional network
Research question — how to identify causal regulatory interactions from time-lagged transcription factor-target relationships
Contribution — DELAYterm·ai, a supervised network that overcomes Granger-based limits such as inferring cyclic regulation
phenomenon · claim-type · cross-domain-transfer · mechanism · claim-relation · replication-status · limitation — grounded in the full text| phenomenon | pseudotime-lagged gene coexpression patterns between causally related genes — “In other words, causally related genes share important pseudotime-lagged patterns of gene coexpression with nearby cells in single-cell trajectories.” | |
| claim-type | method — “We show that combining supervised deep learning with joint probability matrices of pseudotime-lagged trajectories allows the network to overcome important limitations of ordinary Granger causality-based methods, for example, the inability to infer cyclic relationships such as feedback loops.” | |
| method — evidence | — “Here we describe DELAY (short for De picting La gged Causalit y ), a convolutional neural network for the inference of gene-regulatory relationships across pseudotime-ordered single-cell trajectories.” | |
| object of study — evidence | — “Here we describe DELAY (short for De picting La gged Causalit y ), a convolutional neural network for the inference of gene-regulatory relationships across pseudotime-ordered single-cell trajectories.” | |
| research question — evidence | — “Identifying the causal interactions in gene-regulatory networks requires an accurate understanding of the time-lagged relationships between transcription factors and their target genes.” | |
| contribution — evidence | — “We show that combining supervised deep learning with joint probability matrices of pseudotime-lagged trajectories allows the network to overcome important limitations of ordinary Granger causality-based methods, for example, the inability to infer cyclic relationships such as feedback loops.” | |
| mechanism | lagged coexpression features DELAY outperforms because it better learns regulatory features from pseudotime-lagged coexpression matrices. — “Together, these results suggest that DELAY outperforms other methods because it can better learn important gene-regulatory features from pseudotime-lagged gene coexpression matrices.” | mechanisms[0] |
| mechanism | Slingshot dependence Sensitivity to cell order is suspected to come from Slingshot's similarity-based minimum spanning trees. — “We suspect that DELAY is sensitive to the specific ordering of adjacent cells in trajectories because pseudotime inference methods such as Slingshot ( 6 ) infer lineages from minimum spanning trees that directly depend on cell-to-cell similarities in gene expression values” | mechanisms[1] |
| sense | causality - here a direct regulatory interaction of a transcription factor on a target, as opposed to mere correlation — “DELAY apparently recognizes the differences between input matrices of causally related genes and those of merely correlated genes” | senses[0] |
| relation | scRNA-seq —same-as→ single-cell RNA sequencing — (scRNA-seq) — “Our network outperforms several common methods for inferring gene regulation and, when given partial ground-truth labels, predicts novel regulatory networks from single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) data sets.” — “from single-cell RNA sequencing (scRNA-seq)” | concept_edges[0] |
| weight | scRNA-seq —same-as→ single-cell RNA sequencing — 0.95 — “Our network outperforms several common methods for inferring gene regulation and, when given partial ground-truth labels, predicts novel regulatory networks from single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) data sets.” | concept_edges[0] |
| relation | DELAY —contrasted-with→ Granger causality — Unlike — “Unlike ordinary Granger causality, DELAY performed well across cyclic interactions including 2-cycles (MIs) and 3-cycles (FBLs), which are ubiquitous among gene-regulatory networks” — “Unlike ordinary Granger causality, DELAY performed well across cyclic interactions” | concept_edges[1] |
| weight | DELAY —contrasted-with→ Granger causality — 0.85 — “Unlike ordinary Granger causality, DELAY performed well across cyclic interactions including 2-cycles (MIs) and 3-cycles (FBLs), which are ubiquitous among gene-regulatory networks” | concept_edges[1] |
| effect direction | neighbor genes —causes→ false-positive predictions — reduces — “Because the neighbor genes can present stronger alternatives to the primary hypothesis that A directly regulates B , including these matrices for A and B versus their highly correlated transcription factors reduced false-positive predictions.” — “reduced false-positive predictions” | concept_edges[2] |
| relation | neighbor genes —causes→ false-positive predictions — reduced — “Because the neighbor genes can present stronger alternatives to the primary hypothesis that A directly regulates B , including these matrices for A and B versus their highly correlated transcription factors reduced false-positive predictions.” | concept_edges[2] |
| weight | neighbor genes —causes→ false-positive predictions — 0.8 — “Because the neighbor genes can present stronger alternatives to the primary hypothesis that A directly regulates B , including these matrices for A and B versus their highly correlated transcription factors reduced false-positive predictions.” | concept_edges[2] |
| effect direction | pseudotime-lagged matrix —causes→ performance — increases — “Training DELAY on at least one pseudotime-lagged matrix or with at least one neighbor gene greatly increased the network's performance across all data sets.” — “greatly increased the network's performance” | concept_edges[3] |
| relation | pseudotime-lagged matrix —causes→ performance — greatly increased — “Training DELAY on at least one pseudotime-lagged matrix or with at least one neighbor gene greatly increased the network's performance across all data sets.” | concept_edges[3] |
| weight | pseudotime-lagged matrix —causes→ performance — 0.8 — “Training DELAY on at least one pseudotime-lagged matrix or with at least one neighbor gene greatly increased the network's performance across all data sets.” | concept_edges[3] |
| evidential status | mechanism — these results suggest — “Together, these results suggest that DELAY outperforms other methods because it can better learn important gene-regulatory features from pseudotime-lagged gene coexpression matrices.” | mechanisms[0] |
| evidential status | mechanism — We suspect — “We suspect that DELAY is sensitive to the specific ordering of adjacent cells in trajectories because pseudotime inference methods such as Slingshot ( 6 ) infer lineages from minimum spanning trees that directly depend on cell-to-cell similarities in gene expression values” | mechanisms[1] |
| replication | original | |
| limitation | Deep learning methods untrained on a new data set performed like random predictors. — “We inferred the network by the previous methods and found that tree- and mutual information-based methods performed slightly better than deep learning-based methods, which were not initially trained on the new data and consequently performed comparably to random predictors” | limitations[0] |
| limitation | DELAY remains sensitive to the ordering of cells in pseudotime. — “Although our convolutional neural network DELAY remains sensitive to the ordering of single cells in pseudotime, the network can nevertheless accurately infer direct gene-regulatory interactions from both time course and snapshot data sets” | limitations[1] |
| limitation | More than two neighbor-gene matrices sometimes decreased performance. — “training the network on more than two neighbor gene matrices per gene decreased performance in some instances” | limitations[2] |
| limitation | Thirteen reported interactions were not detected by DELAY. — “We additionally identified 13 reported interactions that were not detected by DELAY” | limitations[3] |
| limitation | DELAY is better at discriminating false positives than selecting true positives. — “DELAY performed slightly worse by the metric of area under the PR curve than by area under the ROC curve, indicating that the network is better at discriminating false positives than selecting true positives” | limitations[4] |
finding · quantity · variable · definition · key-term — grounded in the full text| finding | DELAY outperformed eight popular gene-regulatory inference methods. — “After training DELAY on the combined 70% splits, the network outperformed eight of the most popular approaches for inferring gene-regulatory relationships” | content_tags.findings[0] |
| finding | Fine-tuning on 1% of new examples matched or beat unsupervised methods. — “After training each network on only 1% of the new examples and validating it on the remaining 99%, the networks’ performance matched or exceeded that of the unsupervised methods” | content_tags.findings[1] |
| finding | DELAY relies heavily on the ordering of adjacent cells. — “This result suggests that the network relies heavily on the specific ordering of adjacent cells in each trajectory.” | content_tags.findings[2] |
| finding | DELAY is more sensitive to cell down-sampling than to dropouts in low-expressing cells. — “we also discovered that the network was more sensitive to down-sampling of single cells than to gene expression losses in low-expressing cells alone” | content_tags.findings[3] |
| finding | DELAY handled cyclic interactions that ordinary Granger causality cannot. — “Unlike ordinary Granger causality, DELAY performed well across cyclic interactions including 2-cycles (MIs) and 3-cycles (FBLs), which are ubiquitous among gene-regulatory networks” | content_tags.findings[4] |
| finding | DELAY performs best on zinc-finger proteins, bHLH factors and chromatin remodelers. — “We found that DELAY performs best on zinc-finger proteins, bHLH factors, and other chromatin remodelers.” | content_tags.findings[5] |
| finding | Motif analysis uncovered a likely Fiz1 binding sequence. — “Discriminative motif analysis of these genes’ enhancer sequences (Fig. 5 C) uncovered a likely Fiz1 consensus DNA-binding sequence (5′-CGCTGC-3′) similar to that of other Sox2 regulators from the Zic family of C 2 H 2 zinc-finger proteins” | content_tags.findings[6] |
| quantity | 347 — — “The resulting network consisted of 347 predicted transcription factor–target interactions with gene regulation probabilities >0.5.” | content_tags.quantities[0] |
| quantity | 22 of 28 — — “Twenty-two of 28 predicted interactions were confirmed by the accessibility of transcription factor-binding sites.” | content_tags.quantities[1] |
| quantity | 230 ± 1 min — — “With previously determined hyperparameters for scRNA-seq data, training DELAY on two graphics processing units (GPUs) required 230 ± 1 min (mean ± SD) per model.” | content_tags.quantities[2] |
| quantity | P > 0.5 — — “where a probability of P > 0.5 indicates a true gene-regulatory interaction for the given gene pair example” | content_tags.quantities[3] |
| variable | neighbor genes — “training the network on more than two neighbor gene matrices per gene decreased performance in some instances” | content_tags.variables[0] |
| definition | Granger causality tests — “Several methods for gene-regulatory inference rely on pseudotime in Granger causality tests, which try to determine whether new time series can add predictive power to inferred models of gene regulation” | content_tags.definitions[0] |
| definition | snapshots — “Single-cell sequencing technologies can provide detailed data for the investigation of heterogeneous populations of cells collected at specific times—so-called snapshots—during cellular differentiation or dynamic responses to stimulation” | content_tags.definitions[1] |
| key-terms | DELAY, Granger causality, pseudotime, gene-regulatory networks, convolutional neural network, transfer learning, scRNA-seq, scATAC-seq, hair cells, transcription factors |
| concept | DELAY — short for De picting La gged Causalit y — “Here we describe DELAY (short for De picting La gged Causalit y ), a convolutional neural network for the inference of gene-regulatory relationships across pseudotime-ordered single-cell trajectories.” | concepts[0] |
| concept | convolutional neural network — not defined in this paper's own words; DELAY's architecture, an inverted VGGnet — “Here we describe DELAY (short for De picting La gged Causalit y ), a convolutional neural network for the inference of gene-regulatory relationships across pseudotime-ordered single-cell trajectories.” | concepts[1] |
| concept | Granger causality — not defined in this paper's own words; its tests are glossed below — “Unlike ordinary Granger causality, DELAY performed well across cyclic interactions including 2-cycles (MIs) and 3-cycles (FBLs), which are ubiquitous among gene-regulatory networks” | concepts[2] |
| concept | Granger causality tests — try to determine whether new time series can add predictive power to inferred models of gene regulation — “Several methods for gene-regulatory inference rely on pseudotime in Granger causality tests, which try to determine whether new time series can add predictive power to inferred models of gene regulation” | concepts[3] |
| concept | pseudotime — not defined in this paper's own words; the axis along which “cells’ approximate temporal progressions along inferred lineages” are reconstructed — “Many algorithms exploit these cell-to-cell differences to infer dynamic trajectories and reconstruct cells’ approximate temporal progressions along inferred lineages in pseudotime” | concepts[4] |
| concept | snapshots — heterogeneous populations of cells collected at specific times — “Single-cell sequencing technologies can provide detailed data for the investigation of heterogeneous populations of cells collected at specific times—so-called snapshots—during cellular differentiation or dynamic responses to stimulation” | concepts[5] |
| concept | neighbor genes — the two transcription factors with the highest cross-correlation with A and B along the single-cell trajectory — “that is, the two transcription factors with the highest cross-correlation with A and B along the single-cell trajectory” | concepts[6] |
| concept | false-positive predictions — not defined in this paper's own words — “Because the neighbor genes can present stronger alternatives to the primary hypothesis that A directly regulates B , including these matrices for A and B versus their highly correlated transcription factors reduced false-positive predictions.” | concepts[7] |
| concept | pseudotime-lagged matrix — not defined in this paper's own words; a gene coexpression matrix built at a lag in pseudotime — “Training DELAY on at least one pseudotime-lagged matrix or with at least one neighbor gene greatly increased the network's performance across all data sets.” | concepts[8] |
| concept | performance — not defined in this paper's own words; validation metrics over PR and ROC curves — “Training DELAY on at least one pseudotime-lagged matrix or with at least one neighbor gene greatly increased the network's performance across all data sets.” | concepts[9] |
| concept | scRNA-seq — single-cell RNA sequencing — “Our network outperforms several common methods for inferring gene regulation and, when given partial ground-truth labels, predicts novel regulatory networks from single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) data sets.” | concepts[10] |
| concept | single-cell RNA sequencing — not defined in this paper's own words — “Our network outperforms several common methods for inferring gene regulation and, when given partial ground-truth labels, predicts novel regulatory networks from single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) data sets.” | concepts[11] |
| concept | guilt by association — compared sequences enriched in the enhancers of cofactors’ predicted targets to several databases of known DNA-binding motifs to identify the cofactors’ DNA-binding partners — “we compared sequences enriched in the enhancers of cofactors’ predicted targets to several databases of known DNA-binding motifs to identify the cofactors’ DNA-binding partners (so-called guilt by association)” | concepts[12] |
Full text — tagged
The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.
Depicting pseudotime-lagged causality across single-cell trajectories for accurate gene-regulatory inference Abstract Identifying the causal interactions in gene-regulatory networksterm·key requires an accurate understanding of the time-lagged relationships between transcription factorsterm·key and their target genes. Here we describe DELAYterm·ai (short for De picting La gged Causalit y ), a convolutional neural networkterm·ai for the inference of gene-regulatory relationships across pseudotime-ordered single-cell trajectories. We show that combining supervised deep learning with joint probability matrices of pseudotime-lagged trajectories allows the network to overcome important limitations of ordinary Granger causality-based methods, for example, the inability to infer cyclic relationships such as feedback loops. Our network outperforms several common methods for inferring gene regulation and, when given partial ground-truth labels, predicts novel regulatory networks from single-cell RNA sequencingterm·mth (scRNA-seqterm·mth) and single-cell ATAC sequencing (scATAC-seqterm·key) data sets. To validate this approach, we used DELAYterm·ai to identify important genes and modules in the regulatory network of auditory hair cellsterm·key, as well as likely DNA-binding partners for two hair cellterm·key cofactors (Hist1h1c and Ccnd1) and a novel binding sequence for the hair cell-specific transcription factorterm·key Fiz1. We provide an easy-to-use implementation of DELAYterm·ai under an open-source license at https://github.com/calebclayreagor/DELAYterm·ai .
The sequencing of genes expressed in single cells provides detailed information about the developmental programs that define the identities and states of various cell types, but few computational methods can use the dynamic information encoded in these representations to identify causal mechanisms. By exploiting advances in machine learning, we develop a deep neural network that learns from temporal features of gene regulation to identify direct regulatory interactions between transcription factorsterm·key and their target genes. Our method provides mechanistic insights into the development and function of cells by generating high-confidence predictions for interactions in complex gene-regulatory networksterm·key.
Introduction Single-cell sequencing technologies can provide detailed data for the investigation of heterogeneous populations of cells collected at specific times—so-called snapshotsterm·dom—during cellular differentiation or dynamic responses to stimulation ( 1 ). However, owing to inherent delaysterm·ai in molecular processes such as transcription and translation, static measurements from individual cells cannot reveal the causal interactions governing cells’ dynamic responses to developmental and environmental cues ( 2–4 ). Because population-level heterogeneity in tissues often reflects the asynchronous progression of single cells through time-dependent processes, observed patterns of gene expression can nonetheless indicate the stages of development to which individual cells belong ( 5 ). Many algorithms exploit these cell-to-cell differences to infer dynamic trajectories and reconstruct cells’ approximate temporal progressions along inferred lineages in pseudotimeterm·mth ( 6 , 7 ).
Several methods for gene-regulatory inference rely on pseudotimeterm·mth in Granger causality teststerm·mth, which try to determine whether new time series can add predictive power to inferred models of gene regulation ( 8 , 9 ). However, Granger causality-based methods can be error-prone when genes display nonlinear or cyclic interactions or when the sampling rate is uneven or too low ( 9–12 ). Because pseudotimeterm·mth trajectories exhibit these problems, Granger causality-based methods often underperform model-free approaches that exploit pure statistical dependencies in gene expression data ( 9 , 13 , 14 ).
By contrast, deep learning-based methods make no assumptions about the temporal relationships or connectivity between genes in complex regulatory networks; instead, these data-driven approaches learn general features of regulatory interactions ( 15 , 16 ). Here, we describe a deep learning-based method termed De picting La gged Causalit y (DELAYterm·ai) that learns gene-regulatory interactions from discrete joint probability matrices of paired, pseudotime-lagged gene expression trajectories. Our data suggest that DELAYterm·ai can address many shortcomings of current Granger causality-based methods and provide a useful, complementary approach to overcome common limitations in the inference of gene-regulatory networksterm·key from single-cell data.
Results A convolutional neural networkterm·ai predicts pseudotime-lagged gene-regulatory relationships To predict gene-regulatory relationships from single-cell data, we developed a convolutional neural networkterm·ai based on Granger causalityterm·mth ( 17 ). The input to DELAYterm·ai consisted of stacks of two-dimensional joint probability matrices for pairs of transcription factorsterm·key, putative target genes, and highly correlated “neighbor” genes ( 18 ). We constructed the input matrices by aligning the gene expression trajectories of a transcription factorterm·key A at several lagged positions in pseudotimeterm·mth relative to a target gene B to generate joint probability distributions from two-dimensional histograms of gene coexpression at each lag (Fig. 1 A). The value of a given lag indicated the ordinal difference between cells’ rank-ordered positions in pseudotimeterm·mth. Each input matrix consisted of the L2-normalized cell number counts from a two-dimensional gene coexpression histogram with 32 fixed-width bins in each dimension, spanning each gene's minimum and maximum expression values. Although the marginal probability distributions for both A and B remained essentially unchanged at each lag except for cells lost at the leading and lagging edges of the shifted trajectories, realigning the gene expression in pseudotimeterm·mth altered key features of the resulting joint probability matrices. In other words, causally related genes share important pseudotime-lagged patterns of gene coexpression with nearby cells in single-cell trajectories.
Pseudotime-lagged causality allows accurate inference of gene-regulatory networksterm·key. A) Shifting the gene expression trajectory of transcription factorterm·key A forward in pseudotimeterm·mth with respect to that of target gene B generates a series of unique joint probability distributions. The resulting joint probability matrices contain gene coexpression signatures—in the upper-left, upper-right, and lower-right regions—that indicate A directly regulates B . B) The inverted architecture of DELAYterm·ai is wide at the beginning of the network and progressively narrows to a one-dimensional vector followed by a linear classifier and sigmoid activation function to generate gene regulation probabilities. The network uses leaky ReLU activations and padded 3 × 3 convolutions throughout. C and D) DELAYterm·ai outperforms eight of the most popular methods for inferring gene-regulatory interactions across several benchmark scRNA-seqterm·mth data sets, including six unsupervised methods (C) and two supervised methods (D). E) All three supervised methods perform slightly worse upon cross-validation, but DELAYterm·ai still outperforms all other methods on average, with the exception of PPCOR for a single metric. F and G) DELAYterm·ai does not immediately generalize to a new testing data set (F) but outperforms all other methods if fine-tuned on a small percentage of the new data (G). Values for the area under the precision–recall curve (AUPRC) in C), D), F), and G) are normalized by the proportion of positive examples per data set (horizontal lines, E) and averaged across five model replicates for supervised methods. Boxes in E) show the first quartiles, means, and third quartiles for the per- transcription factorterm·key AUPRC values across the k = 5 validation folds. The statistical significance between DELAYterm·ai and the next-best neural network was assessed using a one-sided Wilcoxon signed-rank test (* P ≤ 0.05; ** P ≤ 0.01; *** P ≤ 0.001).
Using ground-truth labels from cell type–specific chromatin immunoprecipitation sequencing (ChIP-seq) data ( 14 ), we conducted supervised learning to train our neural network to predict whether A directly regulates B . This procedure resembled a regression in a Granger causality testterm·mth ( 17 ), in which values of a time series Y at timepoints y t are regressed against values from another time series X at timepoints x t , x t − 1 , x t − 2 , … , x t − T up to some maximum lag T to determine whether any time-lagged values of X add explanatory power to Y ’s autoregressive model. DELAYterm·ai likewise learned higher weights for gene coexpression matrices at specific pseudotimeterm·mth lags that indicated the true regulatory relationship between genes. After comparing several neural network architectures, we selected a six-layered convolutional network trained on pseudotime-aligned ( T = 0 ) and five pseudotime-lagged ( T ∈ { 1 , 2 , 3 , 4 , 5 } ) gene coexpression matrices to predict direct gene-regulatory relationships (Fig. 1 B). We also trained the network on lagged coexpression matrices of two highly correlated neighbor genesterm·mth per gene, that is, the two transcription factorsterm·key with the highest cross-correlation with A and B along the single-cell trajectory. Because the neighbor genesterm·mth can present stronger alternatives to the primary hypothesis that A directly regulates B , including these matrices for A and B versus their highly correlated transcription factorsterm·key reduced false-positive predictionsterm·mth.
DELAYterm·ai outperforms several common methods of gene-regulatory inference We trained DELAYterm·ai on gene expression data sets from human embryonic stem cells (hESCs) ( 19 ), mouse embryonic stem cells (mESCs) ( 20 ), and three lineages of mouse hematopoietic stem cells (mHSCs) ( 21 ). We generated separate training data sets for each of the hematopoietic lineages, and all data sets contained at least 400 cells per lineage ( Table 1 ). Each trajectory was oriented according to known experimental timepoints or precursor cell types and lineages, and pseudotimeterm·mth values were inferred separately with Slingshot ( 6 ) for each lineage (Fig. S1 ). We chose these data sets because the cell types and trajectories are well characterized and offer cell type–specific ChIP-seq data to generate ground-truth networks ( 14 ). Although the three hematopoietic data sets contained similar numbers of examples of true regulation and no interaction, both of the embryonic data sets were class-imbalanced and contained fewer examples of true regulation. To maximize the network's generalizability, we trained DELAYterm·ai simultaneously on all five gene expression data sets but validated it on each data set individually. We first generated randomly segregated 70–30% splits of all possible gene pair examples for each data set and then merged the 70% splits to create a combined training data set. In a separate analysis, we also cross-validated DELAYterm·ai using inductive splits wherein specific transcription factorsterm·key appeared in only the training or validation splits.
Summary of the data sets used to train, test, and evaluate the neural network.
After training DELAYterm·ai on the combined 70% splits, the network outperformed eight of the most popular approaches for inferring gene-regulatory relationships, including six unsupervised methods ( 8 , 22–26 ) (Fig. 1 C) and two deep convolutional neural networksterm·ai ( 18 , 27 ) (Fig. 1 D). We measured the performanceterm·mth of the six unsupervised methods across the combined training and validation splits for each data set and then compared the results with the performanceterm·mth of the supervised methods on the held-out examples alone. With one exception, the deep learning-based methods outperformed all others according to the areas under both the precision–recall (PR) and receiver operating characteristic (ROC) curves (Fig. S2 ). Moreover, DELAYterm·ai outperformed all the other methods according to both metrics, even though one of the deep learning-based methods, DeepDRIM, was trained on 5-fold as many neighbor geneterm·mth matrices. Upon separate cross-validation, DELAYterm·ai performed slightly worse on inductive splits but with a single exception outperformed all other methods on average (Fig. 1 E). Across three of the training data sets, DELAYterm·ai also outperformed a modified version of DeepDRIM that was trained on the same pseudotime-lagged input matrices as DELAYterm·ai (Fig. S2 ). Together, these results suggest that DELAYterm·ai outperforms other methods because it can better learn important gene-regulatory features from pseudotime-lagged gene coexpression matrices.
Transfer learningterm·key allows DELAYterm·ai to predict novel gene-regulatory networksterm·key from new single-cell data sets To test whether DELAYterm·ai generalizes to new data sets, we examined the human hepatocyte (hHep) gene-regulatory networkterm·key using an additional data set with over 400 single cells and known ground-truth interactions from ChIP-seq data ( 14 , 28 ). We inferred the network by the previous methods and found that tree- and mutual information-based methods performed slightly better than deep learning-based methods, which were not initially trained on the new data and consequently performed comparably to random predictors (Fig. 1 F). To determine whether this lack of generalizability arises from batch effects in the single-cell data, we used transfer learningterm·key to fine-tune the three deep learning-based methods on a small number of randomly segregated examples from the new data set while validating on the rest of the examples (Fig. 1 G). After training each network on only 1% of the new examples and validating it on the remaining 99%, the networks’ performanceterm·mth matched or exceeded that of the unsupervised methods; training on up to 20% of the new examples yielded further performanceterm·mth increases. DELAYterm·ai again outperformed every other deep learning-based method. During fine-tuning, we also observed a monotonic increase in the value of the mean validation metric across PR and ROC curves, suggesting that the networks can learn useful features of gene regulation after longer training periods. Moreover, these results demonstrate that DELAYterm·ai can accurately predict gene-regulatory relationships from new data sets with partially known ground-truth labels.
DELAYterm·ai performs well with new input configurations, augmented matrices, and modified data sets Using different numbers of pseudotime-lagged gene coexpression matrices or neighbor geneterm·mth matrices, as well as examples from modified data sets with fewer cells or additional gene dropouts, we next examined the performanceterm·mth of DELAYterm·ai across various input configurations. We employed the original data sets to train new models on gene coexpression matrices of up to 10 pseudotimeterm·mth lags (Fig. 2 A) or up to 10 neighbor genesterm·mth, as well as on gene coexpression matrices with varying dimensions and resolutions (Fig. S3 ). Training DELAYterm·ai on at least one pseudotime-lagged matrixterm·mth or with at least one neighbor geneterm·mth greatly increased the network's performanceterm·mth across all data sets. Although training DELAYterm·ai on up to 10 pseudotime-lagged matrices resulted in comparable or slightly better performanceterm·mth across all data sets, training the network on more than two neighbor geneterm·mth matrices per gene decreased performanceterm·mth in some instances. Adding channel masks for specific lagged matrices suggested that DELAYterm·ai relies on redundancies across all available lagged inputs (Fig. S3 ).
The neural network relies on cell order in pseudotimeterm·mth and gene expression strength. A) Training new models of DELAYterm·ai on increasing numbers of pseudotime-lagged matrices gives the largest performanceterm·mth increase when using up to a single lag. B) Reordering single cells in pseudotimeterm·mth sharply decreases performanceterm·mth, suggesting that DELAYterm·ai relies on the specific ordering of adjacent cells in each trajectory. C and D) The network is sensitive to random down-sampling of cells across data sets (C) but relatively more robust to induced, additional gene dropouts in weakly expressing cells (D), suggesting that DELAYterm·ai relies heavily on highly expressing cells. E and F) The network learns larger input weights for lagged matrices of the first pseudotimeterm·mth lag (E), which also contain more cells in the combined “ON” region ( A ON ∪ B ON ; dotted outline, F) on average across training data sets (dotted line, E). The combined “ON” region is comprised of the upper-left “ON–OFF” quadrant ( A ON ∩ B OFF ), upper-right “ON–ON” quadrant ( A ON ∩ B ON ), and lower-right “OFF–ON” quadrant ( A OFF ∩ B ON ). G) Masking different regions of the input matrices shows that the network relies heavily on the combined “ON” region. The lines and shaded regions in A–E) show the average and full range of values across five model replicates, and the markers in G) show the average values across model replicates. The statistical significance in G) was assessed with a Kruskal–Wallis test (*** P ≤ 0.001).
We also characterized DELAYterm·ai's performanceterm·mth on pseudotime-shuffled trajectories and observed a sharp decrease in performanceterm·mth after reordering nearby cells in each trajectory (Fig. 2 B). This result suggests that the network relies heavily on the specific ordering of adjacent cells in each trajectory. Upon examining DELAYterm·ai's performanceterm·mth on modified data sets containing fewer single cells (Fig. 2 C) or additional gene-dropout noise (Fig. 2 D), we also discovered that the network was more sensitive to down-sampling of single cells than to gene expression losses in low-expressing cells alone. These results indicate that DELAYterm·ai relies more heavily on highly expressing cells and therefore assigns larger input weights for the first lagged input because on average that input contains stronger features of gene activation than other lags (Fig. 2 E and F). To investigate this hypothesis, we used augmented input matrices with masked regions to show that DELAYterm·ai relies heavily on the combined upper-left, upper-right, and lower-right regions ( A ON ∪ B ON ) across all lagged matrices (Fig. 2 G). By performing a post hoc analysis of the statistical dispersion across correctly inferred gene pair examples, we additionally found that DELAYterm·ai generally performs better on transcription factorsterm·key with stable gene expression along single-cell trajectories (Fig. S3 ).
DELAYterm·ai recognizes causal relationships in hierarchical and cyclic gene-regulatory interactions To further investigate how the selection of neighbor genesterm·mth can alter DELAYterm·ai's internal representations and subsequent gene-regulatory inferences, we modified gene pair examples in the original 30% validation splits to exclude all neighbor genesterm·mth C with known interactions across several classes of three-gene motifs. First, we used DELAYterm·ai to infer gene regulation for potential hierarchical and cyclic interactions across validation set examples (Fig. 3 A), including mutual interactions (MIs), feedback loops (FBLs), and three classes of feedforward loops (FFLs). Unlike ordinary Granger causalityterm·mth, DELAYterm·ai performed well across cyclic interactions including 2-cycles (MIs) and 3-cycles (FBLs), which are ubiquitous among gene-regulatory networksterm·key ( 29 ). We next excluded from the same validation set examples all input matrices for any neighbor geneterm·mth of either A or B known from the ChIP-seq data to be involved in potential FBLs or FFLs, replacing the inputs for all omitted neighbors C with the matrices of other highly correlated neighbor genesterm·mth. Because we saw a significant decrease in performanceterm·mth upon excluding these neighbors C in most cases of shared and sequential gene regulation, but not for downstream targets, DELAYterm·ai apparently recognizes the differences between input matrices of causally related genes and those of merely correlated genes. Moreover, these results suggest that although our network uses principles of Granger causalityterm·mth, it can achieve true causal inference for genes involved in several classes of hierarchical and cyclic regulatory interactions by avoiding several limiting assumptions of ordinary Granger causalityterm·mth.
DELAYterm·ai uncovers causal, cyclic, and context-specific gene-regulatory relationships. A) The network's performanceterm·mth when inferring putative transcription factor-target gene pairs across two- and three-gene motifs in the validation set suggests that DELAYterm·ai can distinguish between different types of hierarchical and cyclic gene regulation. Upon exclusion of one or more neighbor genesterm·mth C from the input features of the validation examples, the network's performanceterm·mth declines significantly if C is a shared or sequential regulator of A and B , but not if C is a shared target. B) GO-term enrichment indicates that DELAYterm·ai performs best on zinc-finger proteins (PHD-type, GATA-type, C2H2-type, NHR-type, and C5HC2-type), bHLH factors (Myc-type), and other chromatin remodelers (SNF2-related, chromodomain, bromodomain, HDA domain, and helicase). Markers and error bars in A) show the average values and full range of performanceterm·mth across five model replicates. The statistical significance in A) was assessed with a one-sided Wilcoxon signed-rank test (* P ≤ 0.05) and in B) with a Fisher's exact test ( P adjusted ≤ 0.05).
DELAYterm·ai performs best on zinc-finger proteins, bHLH factors, and other chromatin remodelers To explore how DELAYterm·ai performs on classes of transcription factorsterm·key containing different types of DNA-binding domains, we analyzed the enrichment of gene ontology (GO) terms across correctly inferred transcription factorsterm·key (Fig. 3 B). We found that DELAYterm·ai performs best on zinc-finger proteins, bHLH factors, and other chromatin remodelers. Although DELAYterm·ai performs well on C2H2-type zinc-finger proteins in terms of the total number of correct predictions, we found that it performs better on other chromatin remodelers and plant homeodomain (PHD) zinc-finger proteins by the overall term enrichment (Table S1 ). Interestingly, these results indicate that training the network on cell type–specific ChIP-seq data allows the network to identify regulatory relationships involving some non sequence-specific transcription factorsterm·key and cofactors that nevertheless associate with specific targets at preferred chromatin conformations.
Predicting the gene-regulatory networkterm·key of auditory hair cellsterm·key through multi-omic transfer learningterm·key We next sought to generate predictions for a cell type with complex but incompletely characterized gene-regulatory dynamics. We devised a pipeline for multi-omic transfer learningterm·key to infer the gene-regulatory networkterm·key for developing mouse outer hair cellsterm·key, the mechanical amplifiers of the inner ear. We used both gene expression ( 30 ) and chromatin accessibility ( 31 ) data sets to fine-tune our network. By first calculating the cell-by-gene accessibility scores across annotated genes, we were able to generate lagged input matrices for the single-cell ATAC sequencing (scATAC-seqterm·key) data set (Fig. 4 A). Because the underlying gene expression and gene accessibility distributions are both zero-inflated ( 8 ), the resulting coaccessibility matrices are qualitatively similar to gene coexpression matrices.
DELAYterm·ai accurately predicts gene-regulatory interactions in the auditory hair cellterm·key network through multi-omic transfer learningterm·key. A) An example of an empirical joint probability matrix from scATAC-seqterm·key trajectories. B) After fine-tuning the network on lagged scATAC-seqterm·key input matrices, DELAYterm·ai performs comparably with other neural networks and to previous testing on hHep scRNA-Sseq data. C) Training DELAYterm·ai separately on scRNA-seqterm·mth and scATAC-seqterm·key data sets of hair cellterm·key development reveals a stronger correlation between predictions from model replicates than between data sets across the transcription factor-only network. D) Average gene regulation probabilities P across the hair cellterm·key gene-regulatory networkterm·key accurately predict interactions between transcription factorsterm·key (rows) and targets (columns), when comparing known binding sites to open-chromatin peaks in the target genes’ enhancer sequences. Hierarchical clustering with WPGMA reveals distinct gene modules for prosensory genes ( Sox2 , Id2 , Hes1 , and Prox1 ) and hair cellterm·key genes ( Atoh1 , Pou4f3 , Gfi1 , Lhx3 , and Barhl1 ). Up- or down-regulation was deduced from the correlation in the gene expression data. Markers and error bars in B) show the average values and full range of performanceterm·mth across five model replicates for each neural network. Markers in C) show the median target gene rank correlations across comparisons. The statistical significance in C) was assessed with a Kruskal–Wallis test (*** P ≤ 0.001).
To determine whether pseudotime-lagged gene coaccessibility matrices contain features that indicate direct gene-regulatory relationships, we fine-tuned DELAYterm·ai on the scATAC-seqterm·key matrices using ground-truth examples collected from two cell type–specific data sets of Sox2 and Atoh1 target genes ( 32 , 33 ). Seventy percent of the randomly segregated ground-truth examples were used for fine-tuning, and the remaining 30% were held out for validation (Fig. 4 B). DELAYterm·ai performed slightly worse by the metric of area under the PR curve than by area under the ROC curve, indicating that the network is better at discriminating false positives than selecting true positives. Although DELAYterm·ai did not outperform all other deep learning-based methods, the network's performanceterm·mth was comparable with previous training and validation on small fractions of the hHep single-cell RNA sequencingterm·mth (scRNA-seqterm·mth) data set, which suggests that gene coaccessibility matrices are also useful for inferring direct gene-regulatory interactions.
We separately fine-tuned the network on all available Sox2 and Atoh1 ground-truth examples from an additional gene expression data set for mouse outer hair cellsterm·key ( 30 ). With previously determined hyperparameters for scRNA-seqterm·mth data, training DELAYterm·ai on two graphics processing units (GPUs) required 230 ± 1 min (mean ± SD) per model. We then compared the target gene rank correlations between the inferred transcription factor-only gene-regulatory networksterm·key to determine whether networks inferred from scRNA-seqterm·mth data and scATAC-seqterm·key data were similar (Fig. 4 C). Although we observed stronger correlations between predictions from data set-specific models than between average predictions across data sets, the target gene rank correlations between the two data sets were highly variable, and the predictions of some transcription factorsterm·key agreed better than others. Reasoning that the inferences with better agreement across data sets constituted the best predictions and highest confidence interactions, we derived the consensus hair cellterm·key gene-regulatory networkterm·key from the average predictions across both data sets. The resulting network consisted of 347 predicted transcription factorterm·key–target interactions with gene regulation probabilities >0.5.
DELAYterm·ai identifies important genes, interactions, and modules in the hair cellterm·key gene-regulatory networkterm·key We used hierarchical clustering to group transcription factorsterm·key and target genes in the transcription factor-only gene-regulatory networkterm·key for hair cellsterm·key by similarities in their predicted targets and regulators, respectively (Fig. 4 D). This procedure revealed two distinct developmental modules corresponding to prosensory genes such as Sox2 , Id2 , Hes1 , and Prox1 and hair cell-specific genes such as Atoh1 , Pou4f3 , Gfi1 , Lhx3 , and Barhl1 . We sought to validate the predicted interactions by comparing the locations of known transcription factor-binding sites ( 34 , 35 ) to the locations of open-chromatin peaks ( 31 ) within 50 kb of target genes’ transcription start sites. Twenty-two of 28 predicted interactions were confirmed by the accessibility of transcription factor-binding sites. Of the six remaining interactions, three were instances of predicted target-gene regulation by transcriptional cofactors lacking true DNA-binding domains. We additionally identified 13 reported interactions that were not detected by DELAYterm·ai ( 31 , 36–54 ).
The most notable feature of the inferred gene-regulatory networkterm·key for hair cellsterm·key is that Sox2 and Id2 —in addition to their proteins’ well-known roles in regulating target genes and maintaining a prosensory cell fate—are themselves the targets of a wide variety of transcription factorsterm·key including Sox proteins, homeobox factors, zinc-finger proteins, and Notch effectors. Eleven of the 22 interactions confirmed by our binding-site accessibility analysis represented direct regulation of Sox2 and Id2 , including mutually activating interactions between Sox2 and Sox9 , Sox11 , Prox1 , and Hes1 , and mutual inhibition with Hist1h1c . One other study has suggested that Hes1 directly regulates Sox2 ( 55 ).
Another notable feature of the inferred hair cellterm·key network is that the LIM-homeobox transcription factorterm·key Lhx3 regulates its own expression as well as that of two other LIM-only transcription factorsterm·key (Lmo4 and Lmo1). This result implies a role for Lhx3 in maintaining the later expression of these important early Sox2 inhibitors ( 31 , 47 , 49 ). Other key features of the network include Sox11 regulation by Hes1, Nfkbia regulation by Gfi1, and regulation of the splicing factor gene Srsf2 by the products of three different prosensory genes. Srsf2 was recently predicted to play a role in the splicing of an important deafness gene in humans ( 56 ).
Discriminative motif analysis of target-gene enhancer sequences enables de novo discovery of DNA-binding motifs DELAYterm·ai permits a complementary approach to typical gene-regulatory inference workflows such as SCENIC ( 57 ) that use cis -regulatory sequences to identify and discard false-positive interactions resulting from indirect gene regulation ( 11 , 30 , 31 ). Because these methods rely on databases of known DNA-binding motifs, they necessarily overlook predictions for cofactors and putative transcription factorsterm·key with unknown binding sequences. Instead of using cis -regulatory sequences to refine our initial predictions, we introduced enhancer sequence information post hoc to perform discriminative motif analysis within the enhancers of predicted targets in the hair cellterm·key network (Figs. 5 A and S4 ). We identified enriched motifs that closely resembled known DNA-binding motifs for nine of 11 transcription factorsterm·key with at least one predicted target gene in the transcription factor-only hair cellterm·key network (Table S2 ).
Accurate target-gene predictions enable de novo discovery of DNA-binding motifs. A) Three examples of motifs enriched in the enhancer sequences of predicted target genes (bottom row) closely resemble known DNA-binding motifs (top row) for transcription factorsterm·key with at least one predicted target gene in the transcription factor-only hair cellterm·key gene-regulatory networkterm·key. B and C) Sequences for putative transcription factorsterm·key or cofactors may represent novel DNA-binding motifs or indicate sequence-specific interactions through multi-protein complexes. Motifs for two hair cellterm·key genes closely match known DNA-binding sequences, suggesting that the linker histone Hist1h1c (B, top row) forms complexes with SoxB2/A factors (B, bottom row) and that the putative C2H2 zinc-finger protein Fiz1 (C, top row) recognizes motifs similar to Zic family C2H2 zinc-finger proteins (C, bottom row). The statistical significance of each motif alignment was estimated using the cumulative density function of all possible comparisons of known motifs across enriched sequences for A) or inferred motifs across all database motifs for B) and C) (** P ≤ 0.01; *** P ≤ 0.001).
We additionally sought to predict the mechanisms by which several identified hair cellterm·key cofactors accomplish sequence-specific gene regulation in the absence of true DNA-binding domains. Specifically, we compared sequences enriched in the enhancers of cofactors’ predicted targets to several databases of known DNA-binding motifs to identify the cofactors’ DNA-binding partners (so-called guilt by associationterm·mth). Through this analysis, we were able to determine that the most significantly enriched sequence in the enhancers of the histone Hist1h1c's predicted target genes closely matched known binding motifs for several SoxB2/A transcription factorsterm·key (Fig. 5 B). In addition, a sequence enriched in the enhancers of cyclin Ccnd1's targets accorded with known motifs for the Pbx family of homeobox transcription factorsterm·key (Fig. S5 ). These cofactors might form complexes with the identified transcription factorsterm·key to regulate their target genes in hair cellsterm·key.
As a final demonstration of DELAYterm·ai's high-confidence predictions, we used the fine-tuned model for the scRNA-seqterm·mth data to generate regulatory predictions for transcription factorsterm·key expressed in at least 20% of developing hair cellsterm·key (Table S3 ). After considering these additional transcription factorsterm·key, we identified the putative C2H2 zinc-finger protein Fiz1 as a likely regulator of Sox2 and Hes6 . Discriminative motif analysis of these genes’ enhancer sequences (Fig. 5 C) uncovered a likely Fiz1 consensus DNA-binding sequence (5′-CGCTGC-3′) similar to that of other Sox2 regulators from the Zic family of C 2 H 2 zinc-finger proteins ( 58 ).
Discussion Building upon several deep learning-based methods ( 18 , 27 ), we have demonstrated that combining fully supervised deep learning with joint probability matrices of pseudotime-lagged single-cell trajectories can overcome certain limitations of current Granger causality-based methods of gene-regulatory inference ( 9 , 12 ), such as their inability to infer cyclic regulatory motifs. Although our convolutional neural networkterm·ai DELAYterm·ai remains sensitive to the ordering of single cells in pseudotimeterm·mth, the network can nevertheless accurately infer direct gene-regulatory interactions from both time course and snapshotterm·dom data sets, unlike many supervised methods that rely strongly on the number of available time course samples ( 15 , 16 ). We suspect that DELAYterm·ai is sensitive to the specific ordering of adjacent cells in trajectories because pseudotimeterm·mth inference methods such as Slingshot ( 6 ) infer lineages from minimum spanning trees that directly depend on cell-to-cell similarities in gene expression values ( 7 ).
We have presented a multi-omic paradigm for fine-tuning DELAYterm·ai on both gene expression and chromatin accessibility data sets for the development of auditory hair cellsterm·key. The network's accurate predictions allow de novo inferences of known and unknown DNA-binding motifs, establishing DELAYterm·ai as a complementary approach to common methods of gene-regulatory inference. Our method's high-confidence predictions across the hair cellterm·key gene-regulatory networkterm·key also support it as an attractive option for experimentalists with limited resources to predict true gene-regulatory relationships from complex, large-scale gene-regulatory networksterm·key while avoiding spurious, indirect interactions often introduced through unsupervised methods ( 14 ).
We have additionally identified with high confidence interactions between cofactors and target genes in the regulatory network of hair cellsterm·key that methods such as SCENIC ( 57 ) would have overlooked. Using “guilt by associationterm·mth,” we predict that the cyclin Ccnd1—a known transcriptional cofactor found at the enhancers of several Id genes during retinal development ( 59 )—forms complexes with the Pbx family of homeobox transcription factorsterm·key, of which at least one is a known target of Prox1 in the inner ear ( 60 ). Moreover, we predict that the retina-specific linker histone Hist1h1c ( 61 ) forms complexes with the SoxB2 family of transcription factorsterm·key, which have known antagonistic effects on Sox2 expression in hair cellsterm·key ( 62 ). Finally, we predict that the putative hair cell-specific C2H2 zinc-finger protein Fiz1—which is also expressed during retinal development ( 63 )—is a regulator of Sox2 in hair cellsterm·key and identified its likely DNA-binding sequence in the enhancers of predicted target genes.
Believing that DELAYterm·ai will be a valuable resource to the community, we have provided an easy-to-use and open-source implementation of the algorithm as well as pre-trained model weights for subsequent fine-tuning on new single-cell data sets and partial ground-truth labels. Unlike other deep learning-based methods, our implementation of DELAYterm·ai maximizes usability and prioritizes model flexibility. Although we chose to fine-tune DELAYterm·ai on cell type–specific ChIP-seq data, future studies may choose to fine-tune DELAYterm·ai on curated interactions from databases or from gain- or loss-of-function experiments, especially in the absence of ChIP-seq data ( 14 ). These additional interactions can also supplement smaller ground-truth data sets, such as those of Sox2 and Atoh1 target genes, to mitigate false-negative predictions. We believe that the modest computational cost of training new models of DELAYterm·ai will prove a worthwhile investment for future investigations into complex gene-regulatory networksterm·key.
Materials and methods Preparing single-cell RNA-seq data sets to train and test the convolutional neural networkterm·ai To train our convolutional neural networkterm·ai, we used scRNA-seqterm·mth data sets from two time course studies of endodermal specification of hESCs ( 19 ) and mESCs ( 20 ) and an in vivo snapshotterm·dom study of erythroid (E), granulocyte–monocyte (GM), and lymphoid (L) specification in mHSCs ( 21 ). We additionally employed a data set from a time course study of the differentiation of hHeps from induced pluripotent stem cells to test our network ( 28 ). For each of these studies, we collected the normalized gene expression values, pseudotimeterm·mth values, and ground-truth networks from BEELINE’s supplementary data files ( 14 ). In BEELINE, Pratapa et al. collected the normalized gene expression values from the original studies or obtained the values themselves by log-transforming the transcripts per kilobase million. Moreover, those authors used Slingshot ( 6 ) to infer the pseudotimeterm·mth values separately for each data set, orienting the inferred trajectories by known experimental timepoints or precursor cell types and lineages. Finally, Pratapa et al. selected genes with differential expression across pseudotimeterm·mth by applying generalized additive models (GAMs) to compute gene variances and their associated P values. To train DELAYterm·ai, we utilized their gene variances and corrected P values to select differentially expressed transcription factorsterm·key ( P adjusted < 0 .01 ) and 500 additional genes with the highest variance for each data set. We also collected the ground-truth labels from BEELINE’s supplementary data files, which contained tables of transcription factorterm·key–target gene interactions curated from cell type–specific ChIP-seq experiments in the ENCODE ( 64 ), ChIP-Atlas ( 65 ), and ESCAPE ( 66 ) databases. Table 1 provides details for each data set including descriptions of the ground-truth networks.
We also characterized our network on modified data sets with shuffled pseudotimeterm·mth, fewer cells, and additional gene-dropout noise. We shuffled pseudotimeterm·mth values across single-cell trajectories using np . random . normal in NumPy (v1.20.2) to select the indices of single cells either leading or lagging successive cells in each trajectory at some length scale σ before swapping the cells’ positions in pseudotimeterm·mth. We also used np . random . choice to select the indices of single cells to retain when down-sampling the number of cells in a data set. Lastly, we induced additional gene-dropout noise in data sets by setting the gene expression values to 0 with a probability of p for the bottom p percent of cells expressing each gene in a given data set.
Generating examples of transcription factorterm·key–target gene pairs from single-cell data sets For each data set, we generated mini-batches of 512 transcription factorterm·key–target gene pair examples by first enumerating all possible combinations of differentially expressed transcription factorsterm·key and potential target genes—including both transcription factorsterm·key and highly varying genes. We then generated aligned and pseudotime-lagged gene coexpression matrices for up to five lags with the following configurations as separate input channels: transcription factorterm·key–target, transcription factorterm·key–transcription factorterm·key, target–target, transcription factorterm·key–neighbor (for two neighbors), and target–neighbor (for two neighbors). We concatenated the input channels for each example to form three-dimensional stacks of input matrices with dimensions 42 × 32 × 32 (channels × height × width). For each gene, we used np . correlate in NumPy to select the transcription factorsterm·key with the highest absolute cross-correlation at a maximum offset of five pseudotimeterm·mth lags to use as neighbor genesterm·mth. We trained five model replicates on unique, randomly segregated 70–30% splits of all possible gene pair examples generated with the random _ split function in PyTorch (v1.8.1). Each random split achieved a strict segregation of unique gene pair examples into either the training or validation sets, though gene pairs for the same transcription factorterm·key but different target genes appeared in both sets. For separate model cross-validation, we generated five inductive splits across each data set wherein specific transcription factorsterm·key appeared in only the training (80%) or validation (20%) splits. For both analyses, training splits from individual data sets were subsequently merged to create combined training data sets.
We later used the original randomly segregated validation splits to characterize our network on augmented coexpression matrices and matrices generated from modified single-cell data sets. To augment gene coexpression matrices, we masked with zeros either specific region across all input channels or entire input channels corresponding to specific pseudotimeterm·mth lags.
Constructing a convolutional neural networkterm·ai to classify lagged gene coexpression matrices We designed a convolutional neural networkterm·ai based on an inverted VGGnet ( 67 ) that uses five convolutional layers to first expand the input to 1,024 channels and then successively halve the number of channels to 64 before classifying examples as either true regulation or no interaction with a fully connected linear layer. Each convolutional layer sums over the two-dimensional cross-correlations ( ⋆ ) of 3 × 3 kernels and input channels i to find the features for a given output channel j , as shown in Equation 1 : in which the weights and bias are trainable parameters, N is the mini-batch size, C is the number of channels, C out = C in / 2 for layers 2 through 5, and all convolutions are zero-padded at the edges of matrices. To preserve the gradient flow during training, we used leaky rectified linear units (ReLUs) with a negative slope of 0.2 (Equation 2 ) as our nonlinear activation function after each convolutional layer: As shown in Equation 3 , we additionally used 2 × 2 , unpadded max-pooling layers between convolutional layers to identify important features and down-sample activation maps: in which x represents the 2 × 2 input windows and h , w are the height and width of the input channel j . After the final convolutional layer, we used global-average pooling (Equation 4 ) to reduce the remaining 64 feature maps j to a single, 64-dimensional vector x : We lastly used a fully connected linear layer (Equation 5 ) with a sigmoid activation function (Equation 6 ) to generate gene regulation probabilities: where a probability of P > 0.5 indicates a true gene-regulatory interaction for the given gene pair example.
Training and fine-tuning the network on pseudotime-lagged gene coexpression matrices We used PyTorch's implementation of stochastic gradient descent (SGD) to optimize our network with respect to the class-weighted binary cross-entropy loss, summed across each mini-batch and scaled by the overall batch size 512, as shown in Equation 7 : in which y and y ^ are the target and predicted values (respectively), w N is the fraction of true examples in the mini-batch, and N is the size of the current mini-batch (which might be <512). Prior to optimization, we used He initialization ( 68 ) with uniform priors to set the weights for all convolutional and linear layers. We used a learning rate (LR) of 0.5 to train each model for up to 100 epochs, validating performanceterm·mth after each epoch and stopping training after 10 or more epochs without an improvement in the mean validation metric across PR and ROC curves. For fine-tuning, we trained each model for up to a maximum of 5,000 additional epochs due to the observed monotonic increase in the value of the mean validation metric. We occasionally reduced the LR to 0.25 or 0.1 if the training became unstable, and we tried several mini-batch sizes ≥24 for fine-tuning and separate model cross-validation. All training and testing was performed on two Nvidia RTX 8000 GPUs.
Comparing DELAYterm·ai to other top-performing gene-regulatory inference methods We compared the performanceterm·mth of our neural network to two other convolutional neural networksterm·ai, as well as the top six best-performing methods identified in a previous benchmarking study. Using BEELINE, we inferred the gene-regulatory networksterm·key for all training and testing data sets with the tree-based methods GENIE3 ( 22 ) and GRNBoost2 ( 23 ), the mutual information-based method PIDC ( 24 ), and the partial correlation and regression-based methods PPCOR ( 25 ), SCODE ( 26 ), and SINCERITIES ( 8 ). We utilized the best parameter values identified in BEELINE for the partial correlation and regression-based methods. For the two deep learning-based methods [CNNC ( 27 ) and DeepDRIM ( 18 )], we used the original studies to reconstruct the neural networks in PyTorch before training models on the same randomly segregated training examples as DELAYterm·ai. In addition to training both neural networks on pseudotime-aligned gene coexpression matrices for primary transcription factorterm·key–target gene pairs, we also trained DeepDRIM on 10 neighbor geneterm·mth matrices per gene as in the original study. Moreover, we separately trained DeepDRIM on five pseudotime-lagged matrices per gene pair with 10 neighbors across the same inductive splits as DELAYterm·ai to compare the cross-validated performanceterm·mth between the networks.
Analyzing enrichment of GO terms across correctly inferred gene pair examples We used Enrichr ( 69 ) ( https://maayanlab.cloud/Enrichr/ ) to analyze GO-term enrichment across correctly inferred transcription factorsterm·key for terms related to InterPro DNA-binding domains ( 70 ). We used a Fisher's exact test with Benjamini–Hochberg correction for multiple-hypothesis testing ( P adjusted < 0 .05 ) to assess the statistical significance for each GO term.
Preparing single-cell multi-omics data sets to infer the hair cellterm·key gene-regulatory networkterm·key We used two single-cell data sets from developing mouse outer hair cellsterm·key ( 30 , 31 ) to predict the consensus hair cellterm·key gene-regulatory networkterm·key. First, we collected the normalized gene expression values from the original scRNA-seqterm·mth study ( 30 ) and then used reciprocal PCA to integrate single cells across four timepoints of sensory epithelium development into a single data set in Seurat ( 71 ) (v3.1.4). Then, we used Slingshot ( 6 ) (v1.0.0) to infer pseudotimeterm·mth values across the outer hair cellterm·key trajectory and applied LOESS regressions and GAMs to select the genes that were differentially expressed across both the sensory epithelium and pseudotimeterm·mth, respectively. We again selected all differentially expressed transcription factorsterm·key ( P < 0 .01 after Bonferroni correction for multiple-hypothesis testing) and 500 additional genes with the highest variance and used the integrated gene expression assay in Seurat to generate the corresponding gene coexpression matrices. In a separate analysis, we broadened our selection criteria to encompass all genes that were differentially expressed in at least 20% of the outer hair cellsterm·key, regardless of their expression across the full sensory epithelium.
For the scATAC-seqterm·key data set, we first collected the cell-by-peak chromatin accessibility data from the original study ( 31 ). Then, we used the function createGmatFromMat in SnapATAC ( 72 ) (v1.0.0) to calculate the cell-by-gene accessibility scores as the counts of each 5 kb bin per gene in the UCSC mouse genome mm10 ( 73 ) ( TxDb . Mmusculus . UCSC . mm 10 . knownGene ; v3.4.4) that contained at least one open-chromatin peak for a given cell. These values were then log-normalized. We again used Slingshot to infer pseudotimeterm·mth values across the outer hair cellterm·key trajectory and applied GAMs to the raw accessibility counts to select the differentially accessible genes across pseudotimeterm·mth. For the inferred scATAC-seqterm·key network, we selected all differentially accessible transcription factorsterm·key and variable genes ( P adjusted < 0 .01 ) that were also differentially expressed along the scRNA-seqterm·mth trajectory. Because the scATAC-seqterm·key data set had fewer cells than the scRNA-seqterm·mth data set, we used 24 fixed-width bins in each dimension to generate the corresponding gene coaccessibility matrices prior to fine-tuning.
Discriminative motif analysis to discover de novo transcription factor-binding motifs We used the UCSC mouse genome mm10 and twoBitToFa to download the 100 kb enhancer sequences spanning 50 kb upstream and downstream of all target genes’ transcription start sites in the transcription factor-only hair cellterm·key gene-regulatory networkterm·key and then divided them into 100 bp fragments and sorted the resulting sequences into primary (predicted targets) and control (predicted no interaction) groups for each transcription factorterm·key with at least one target gene in the network. We then used STREME ( 74 ) (v5.4.1) to perform discriminative motif analysis with a P value threshold of 0.05 to identify enriched motifs in predicted target genes’ enhancers, which we then compared with either the transcription factorsterm·key’ known binding motifs from CIS-BP ( 34 ), or to all motifs from the JASPAR ( 75 ), UniPROBE ( 76 ) (mouse), and Jolma et al. ( 77 ) databases using TOMTOM ( 78 ) (v5.4.1) with an E value threshold of 10 for significant alignments.
Supplementary Material Click here for additional data file.