Bayesian Networks Predict Neuronal Transdifferentiation

Abstract

We employ the language of Bayesian networks to systematically construct gene-regulation topologies from deep-sequencing single-nucleus RNA-Seq data for human neurons. From the perspective of the cell-state potential landscape, we identify attractors that correspond closely to different neuron subtypes. Attractors are also recovered for cell states from an independent data set confirming our models accurate description of global genetic regulations across differing cell types of the neocortex (not included in the training data). Our model recovers experimentally confirmed genetic regulations and community analysis reveals genetic associations in common pathways. Via a comprehensive scan of all theoretical three-gene perturbations of gene knockout and overexpression, we discover novel neuronal trans-differrentiation recipes (including perturbations of SATB2, GAD1, POU6F2 and ADARB2) for excitatory projection neuron and inhibitory interneuron subtypes.

Keywords: gene regulation network; neuroscience; systems biology.

Figure 1

Workflow pipeline. A.ii. Hierarchical clustering of 3,227 quality-filtered single neuron data sets from previous single-nucleus RNA sequencing study. Number of 10-fold up-($\uparrow$) and down-regulated (↓) DEGs given at each junction. These are defined as up-regulated in the subsequent left hand branch and down-regulated in the left hand branch each relative to the right hand branch. Terminal clusters I - IV at Level 2 used in this work.(Lake et al. 2016). A.i. 77 ten-fold DEGs (of which 74 are unique) used to train networks. Splits refer to junctions in the hierachical tree in Figure 1 A.ii. Genes in bold occur in multiple splits. B. Data discretization on the weighted arithmetic mean of the log₂(TPM) for each gene across 3,203 Level 2 single neurons. Data downsampled to 1,920 samples. C. Structure learning directed acyclic gene reulation networks using the discretized downsampled data with local and global optimization routines. 20 random seeds used to generate 20 different structures. 5 structures randomly chosen for further calculations. D.i. 450,000 random initialisations of the nodes in the continuous interval [0,1] for the 5 network structures. 2 TBN DBN inference performed for each network in each initial state. Converged attractor states subsequently clustered. D.ii. Nodes in 5 network structures initialised in the continuous interval [0,1] corresponding to the four neuronal subtype cell states. For each network structure in each initial state, 2 TBN DBN inference carried out for all 3-gene perturbation combinations (clamping nodes as overexpressed or knocked out for duration of inference). Subsequent node probabilities averaged over the 5 structures for each 3-gene perturbation in each initial state.

Figure 2

Experimental barcodes for clusters I - IV post data processing. For gene expression probabilities $\ge 0.5$ code displayed as green and for gene expression probabilities $<0.5$ code displayed as red. Genes that display the same pattern across all four clusters are displayed on the same row. Genes are colored red, purple, green and blue for those expressed in clusters I - IV respectively. DEG that are expressed in both PN subtypes I and II are colored magenta and DEG expressed in both IN subtypes III and IV are colored cyan.

Figure 3

Acceptance criterion for edges that are unfavorable to $BIC$ score during simulated annealing stage of structure learning. Representative temperatures in the range $0.1\le \text{T}\le 20$ plotted.

Figure 4

Convergence of $\Delta BIC$ (referenced from structure with 1.1M Nsteps) with NSteps (defined as number of steps during simulated annealing phase). Each value is the arithmetic mean over 5 independent structure runs.

Figure 5

Attractor heat map for all 1082 unique attractors of the unperturbed network model. Cell states defined using quaternary intervals $0.00-0.25,0.25-0.50,0.50-0.75,0.75-1.00$ and hierarchically clustered using the heatmap.2 function in R. Basin sizes for the 3 representative cell state clusters A, B and C given in legend.

Figure 6

Merged GRN for 20 independent BN structures. 228 unique edges with frequency $>0.40$ included (from a total of 870 unique edges summed over all 20 BN). Displayed using Yifan Hu algorithm as implemented in Gephi version 0.9.1. Node size proportional to out-degree and edge width proportional to frequency. Community detection algorithm (Blondel et al. 2008) run with resolution of 1.55. Modularity = 0.421. Number of communities = 4.

Figure 7

The 3 best transdifferentiaion recipes as ranked by RMSD between source cluster I perturbed states and experimental target states. Subplots are target state node probabilities vs. source state node probabilities in the unperturbed (experimental vs. experimental) and perturbed (target experimental vs. source theoretical perturbed) states, represented by cyan and dark blue data points respectively. Best recipes are highlighted in green (overexpression clamp) and red (knockout clamp) and arrows are drawn to show the direction and magnitude of node probability for the given clamp.

Figure 8

Cell state potential changes along a representative path for a single structure under perturbed conditions. Each line represents a different sampling of the conditional probabilities. The best 3-gene recipe for the interconversion between S-I→ T-IV, ERBB4↑ / ADARB2↑ / SATB2↓, is shown.