Intracellular and Intercellular Gene Regulatory Network Inference From Time-Course Individual RNA-Seq

Abstract

Gene regulatory network (GRN) inference is an effective approach to understand the molecular mechanisms underlying biological events. Generally, GRN inference mainly targets intracellular regulatory relationships such as transcription factors and their associated targets. In multicellular organisms, there are both intracellular and intercellular regulatory mechanisms. Thus, we hypothesize that GRNs inferred from time-course individual (whole embryo) RNA-Seq during development can reveal intercellular regulatory relationships (signaling pathways) underlying the development. Here, we conducted time-course bulk RNA-Seq of individual mouse embryos during early development, followed by pseudo-time analysis and GRN inference. The results demonstrated that GRN inference from RNA-Seq with pseudo-time can be applied for individual bulk RNA-Seq similar to scRNA-Seq. Validation using an experimental-source-based database showed that our approach could significantly infer GRN for all transcription factors in the database. Furthermore, the inferred ligand-related and receptor-related downstream genes were significantly overlapped. Thus, the inferred GRN based on whole organism could include intercellular regulatory relationships, which cannot be inferred from scRNA-Seq based only on gene expression data. Overall, inferring GRN from time-course bulk RNA-Seq is an effective approach to understand the regulatory relationships underlying biological events in multicellular organisms.

Keywords: bulk RNA-seq; gene regulatory network; intracellular; mouse; time course.

FIGURE 1

Pseudo-time analysis for time-course individual bulk RNA-Seq of mouse embryos in early development. (A) Transcriptomic trajectory of mouse embryos from E7.5 to E13.5. Each point on the PC1-PC2 plane indicates each individual RNA-Seq result. The solid line indicates an inferred trajectory by slingshot. (B) A scatterplot of pseudo-time and corresponding stage for each sample. Each point indicates each individual RNA-Seq result and is jittered along the Y-axis. (C) An example (Fabp7) of difference in gene expression dynamics along stage and pseudo-time. The black lines indicate smoothing curves for each data. Sum of squared residuals (SSRs) between the observed and fitted values were calculated. (D) A scatterplot of SSR of normalized gene expression of all genes along stage and pseudo-time. The red line indicates the same value between stage and pseudo-time.

FIGURE 2

Inference of downstream genes based on the values in (A). (A) Selection of a parameter D for SCODE. The sum of squared residuals (SSRs) for each D (from 1 to 10 by 0.5) was calculated. We used D = 4. (B) Pearson’s correlation coefficients between the values of each A from 20 optimizations and the meanA, which is the average of each value from all optimizations. (C) Calculation of Pearson’s correlations between the average expression level of each gene and the absolute values in each row and column of A. (D) The violin plots of Pearson’s correlations. (E,F) Scatter plot of the average expression level of each gene and the absolute values in the column (E) and row (F) of A for Sox8. (G) An example of the definition of thresholds for significant regulatory relationship between the regulator and downstream genes. A scatter plot of the values of A for a regulator, Sox8, in decreasing order. The solid line indicates a regression for the scatter plot. The dashed lines indicate the defined thresholds.

FIGURE 3

Validation of the inferred regulatory network of transcription factors with the TF2DNA database. (A) Scatter plot of area under the curve (AUC) target selection for each transcription factor (TF) based on the absolute values of A inferred using SCODE and weight.matrix inferred using dyGENIE3. (B) Scatter plots of the validated inferred target gene rate of E2f3 and Arid5b in the TF2DNA database. Downstream genes of TFs were selected based on the absolute values of A in decreasing order. Solid lines indicate genes at the thresholds. (C) Bar graph of validated and non-validated downstream genes of each TF in the TF2DNA database, in decreasing order of the total number of inferred downstream genes. Only the top 30 TFs are shown. (D) Scatter plot of the validated target gene rate of the inferred downstream genes of TFj and the background rate of target genes in the TF2DNA database. Histograms show the distribution of validated target gene rates of inferred downstream genes and the background rates of downstream genes in the TF2DNA database.

FIGURE 4

Overlap of inferred genes downstream of ligand- and receptor-related genes. (A) Bar plot of the number of inferred downstream genes that are common and unique for each ligand-receptor pair. (B–D) Upset plots of inferred downstream genes that are positively and negatively regulated by the representative ligand- and receptor-related genes. (B) Retinoic acid signaling pathway. (C) Delta/Notch signaling pathway. (D) Hedgehog signaling pathway.

FIGURE 5

Validation of time-course individual RNA-Seq-based GRN inference with a public organ-level RNA-Seq data. (A) Scatter plots of area under the curve (AUC) target selection for each transcription factor (TF) based on the absolute values of A inferred from public time-course organ-level RNA-Seq data using SCODE. (B) Hierarchical clustering of the inferred GRNs from each organ-level RNA-Seq data.