Differential expression analysis reveals the common differentially expressed genes among five NSC differentiation paths
Differential expression analysis was performed on the five preprocessed NSC differentiation paths (NSC_NPC, NPC_AST, NPC_RGL, OPC_OLI, OPC_MO). It yielded a set of DEGs with different numbers (ranging from 566 to 2582) for each path. Interestingly, 12 common DEGs co-occurred in the five NSC differentiation paths (Fig. 3A), among which some genes presented over-expressed (represented by the red bar) and down-expressed (represented by the blue bar) across different NSC differentiation paths, while some genes consistently down-expressed across different NSC differentiation paths. Commonly, the alterations of gene expression reflect on the functional changes in biological processes. Therefore, as shown in Fig. 3A, it implies that the specific genes play different roles in different cell differentiation paths.
To evaluate the biological processes related to these common DEGs, we performed the GO enrichment analysis. The GO enrichment analysis (Fig. 3B) showed that these 12 common DEGs were significantly enriched in differentiation and development of neural system. For example, Cdk5, App, Gsk3b, Ednrb and Ptprz1 were enriched in negative regulation of neuron differentiation(adjusted p-value = 1.14e-7); Cdk5, Dclk1, Gsk3b, Ptprz1were enriched in neuron migration(adjusted p-value = 1.64e-4); Cdk5, App, Dclk1, Gsk3b, Ptprz1 were enriched in axon development(adjusted p-value = 2.04e-4). These GO terms serve as further evidence to support that the identified common DEGs play important roles in NSC differentiation and development.
Although the result of differential expression analysis helped us narrow the candidate gene sets, it failed to assess the importance of these candidate genes and the association among them. In contrast to conventional differential expression analysis, network analysis, especially differential network analysis, allows us to study these complex life processes quantitatively and qualitatively from both bioinformatics and system biology perspectives [19]. To explore which genes, pathways or patterns are responsible for making the cell fate decision of NSC differentiation, comprehensive analyses of differential network were performed in the following sections.
Overlapping differential networks of five paths reveals the candidate regulatory pattern
In order to assess these common DEGs compared with all DEGs in a systematic perspective, three networks (two state networks and one differential network) were separately built for each path, which integrated single-cell transcriptome and proteome data (see Methods for detail). A state network characterizes the association of each pair of genes in a specific cell type, while a differential network characterizes the changed association of each pair of genes from one cell type to another. After overlapping these differential networks of each path, two topologically overlapped subnetworks were identified: subnetwork 1 (Slc1a3, Slc38a2) and subnetwork 2 (Gsk3b, App, Cdk5). Here we focused on the Gsk3b_App_Cdk5 as a candidate RP because there is significant evidence that GSK3 signaling has a critical role in the regulation of neurogenesis, neurodevelopment, and in neuroplasticity [34]. Therefore, the associated genes of Gsk3b in differential networks may play the similar role.
The weight of edges in a differential network measures the changed correlation of each pair of genes and indicates the alterations of regulation between genes. It can be observed that the differences of edges correlation in the RP present an inconsistency in each path (Fig. 4). In path NSC_NPC, NPC_RGL and OPC_OLI, the correlations among Cdk5, Gsk3b and App were weakened, while in path NPC_AST and OPC_MO, the correlations among Cdk5, Gsk3b and App were strengthened, which reflects the specificity of the RP.
According to the recently published literatures relating to the RP Gsk3b_App_Cdk5, all genes in the RP were separately found to be related to NSC differentiation in some degree. For instance, neuronal and glial differentiations of human neural stem cells are regulated by amyloid precursor protein (APP) levels [35]. GSK-3b inhibits degradation of b-catenin and thereby potentiates its downstream signaling, significantly enhancing neuronal differentiation [36]. Fine-tuning of Synapsin III expression and phosphorylation by CDK5 activation through Sema3A activity are essential for proper neuronal migration and orientation [37]. In total, all these literatures have confirmed that, each gene in the RP separately participates in specific stage of NSC differentiation. However, it is not clear whether the RP is involved in the regulation of all differentiation stages of the NSC lineage.
Network diffusion analysis reveals an indirect role of the regulatory pattern
In order to explore those important genes that indirectly influence topological differences (gene signaling propagation), the MC method based on perturbation test was performed on the differential network of each path [28], which is capable of quantifying the indirect influence of perturbating one gene on network propagation and network stability.
The median score of the network diffusion in each path is around 0.75 (the white dots of Fig. 5A), while the score of the most of the 12 common DEGs is greater than the median score (Fig. 5B). The result of network diffusion suggests that the scores of these three genes (Gsk3b, App, Cdk5) in the RP are consistently high, and the result of gene hierarchical clustering based on their scores on network diffusion (Fig. 5B) shows that these three genes play an interactive role in complex life processes.
Differential network analysis of common DEGs shows inconsistent topology in each path
To comprehensively assess the differences of network topology, five kinds of differential network measurement analysis methods (DDC, DEC, DPC, DBC and DiffRank) were performed on each path, and the average score of the five methods is regarded as the overall network topological differences (See Methods for details).
The topological differences of 12 common DEGs in five paths were shown in Fig. 6. There were inconsistencies of network topology of these 12 DEGs during different differentiation paths. When compared to the analysis of network diffusion, most of the 12 DEGs show no significant network topology differences (the score smaller than 0.2), which implies that they play indirect regulatory roles in the five differentiation paths. However, it should be noted that more than one of the top three genes with most significantly topology difference in each path were involved in the RP Gsk3b_App_Cdk5. Therefore, the inconsistent scores of topological differences and consistent high scores of network diffusion reveal that these 3 genes in the RP may alternate their direct or indirect regulatory roles in different differentiation paths in order to promote cell differentiation towards different trajectory development.
KEGG enrichment analysis reveals the common and specific pathways implicated in the regulatory pattern
In order to explore the pathways of the RP Gsk3b_App_Cdk5 and their interactive genes, we performed the KEGG enrichment analysis on the pattern and their first-order neighbors in each differential network. The top ten enriched pathways in each path are shown in (Fig. 7A). Different paths are enriched in pathways related to distinct chemical signaling, which suggests the heterogeneity of NSC lineage differentiation. After intersecting all enriched pathways in each path, four common pathways are identified (Fig. 7B), including axon guidance, Alzheimer’s disease, prostate cancer and endocrine resistance. Importantly, one of the common pathways, axon guidance, refers to the process by which growing neural axons follow specific, predictable paths to reach their target locations, which may play an important role in cell differentiation and development [38]. In this manner, specific regulators are able to establish cell type identity through shaping molecular programs for axon guidance [39], while aberrant axon guidance regularly results in most neurodegenerative diseases, such as Alzheimer’s disease [40]. The above analysis confirms that the identified pattern plays important roles in NSC lineage differentiation and development.