The innate immune response is primarily mediated by the Toll-like receptors functioning through the MyD88-dependent and TRIF-dependent pathways. Despite being widely studied, it is not yet completely understood and systems-level analyses have been lacking. In this study, we identified a high-probability network of genes activated during the innate immune response using a novel approach to analyze time-course gene expression profiles of activated immune cells in combination with a large gene regulatory and protein-protein interaction network. We classified the immune response into three consecutive time-dependent stages and identified the most probable paths between genes showing a significant change in expression at each stage. The resultant network contained several novel and known regulators of the innate immune response, many of which did not show any observable change in expression at the sampled time points. The response network shows the dominance of genes from specific functional classes during different stages of the immune response. It also suggests a role for the protein phosphatase 2a catalytic subunit α in the regulation of the immunoproteasome during the late phase of the response. In order to clarify the differences between the MyD88-dependent and TRIF-dependent pathways in the innate immune response, time-course gene expression profiles from MyD88-knockout and TRIF-knockout dendritic cells were analyzed. Their response networks suggest the dominance of the MyD88-dependent pathway in the innate immune response, and an association of the circadian regulators and immunoproteasomal degradation with the TRIF-dependent pathway. The response network presented here provides the most probable associations between genes expressed in the early and the late phases of the innate immune response, while taking into account the intermediate regulators. We propose that the method described here can also be used in the identification of time-dependent gene sub-networks in other biological systems.
The innate immune system is the primary host response to invading pathogens. The innate immune response is characterized by germline-encoded pattern-recognition receptors (PRRs) that detect and bind to specific microbial components, also known as pathogen-associated molecular patterns (PAMPs). Toll-like receptors (TLRs) are a family of PRRs that are conserved from worm to mammals and expressed on different types of immune cells, such as macrophages, dendritic cells (DCs) and B cells, as well as non-immune cells, such as fibroblasts and epithelial cells. 10 and 13 TLRs have been identified in human and mouse, respectively, each with distinct microbial ligands. The binding of these ligands to their specific receptors triggers downstream signaling cascades causing the expression of pro-inflammatory cytokines, ultimately leading to systemic inflammation. TLRs primarily function through two pathways – the MyD88-dependent pathway which leads to the expression of proinflammatory cytokines, and the TIR-domain–containing adaptor protein-inducing IFN-β (TRIF) -dependent pathway which produces the type I interferons (IFNs) [1], [2]. Though much is known about the pathways activated during the innate immune response, recent perturbation studies have identified previously unknown regulators and transcription factors, highlighting the complexity of the innate immune system and the incompleteness of our current knowledge [3]–[5]. While these studies provide important information about the genes affected on perturbation of a causal gene, they do not explain the cause of the observed expression changes. Additionally, these studies are inherently limited to genes which show changes in expression at the time of observation thus providing an incomplete representation of the activated pathways. The complexity of the innate immune system, the ease of monitoring transcriptional changes, and the availability of large amounts of regulatory and interaction information, all facilitate its analysis using computational methods. An initial computational study mapped all the known interactions associated with the immune response from literature [6]. This study provided a high confidence signaling network and identified the “bow-tie” structure of the immune response. However, it was limited in size and coverage. Li et al. used this signaling map to identify 10 distinct input-output pathways [7]. The resultant modules were further used by Richard et al. to identify a minimum set of genes whose deletion affects the fidelity of the TLR signaling pathways [8]. Though these methods used novel approaches to analyze the TLR signaling pathways, they did not take the temporal changes of the immune response into account. Using a different approach, Seok et al. studied the regulatory networks of 10 transcription factors and their targets using the Network Component Analysis approach [9]. While this study considered the dynamic nature of the immune response through the use of time-course gene expression profiles, it was limited to only 10 transcription factors. Thus, the computational analyses so far performed to study the innate immune response have either been limited by the size of the molecular network used, or by the lack of time-course gene expression profiles. In this study, we perform a comprehensive computational analysis of the dynamic aspects of the innate immune response in the context of a large-scale molecular network. Several methods using condition-specific genetic, transcriptional and epigenomic data in the context of large protein-protein interaction (PPI) and protein-DNA interaction (PDI) networks have been developed, and have led to the identification of novel regulators and pathways in several cellular systems [10], [11]. These include Network Component Analysis (NCA) [12], DREM [13] and its recent update SDREM [14], ResponseNet [15] and SteinerNet [16], [17]. Data from time-course gene expression profiles is particularly informative in this context since it can capture chronological events in the cellular system. However, some of the methods listed above, like ResponseNet and SteinerNet, are insensitive to the temporal aspect of gene expression, while others like NCA and SDREM use the temporal gene expression information only to identify transcription factors activated at various time points but not to predict active networks. Others have used time-course gene expression profiles either to identify time-specific protein-modules in PPI networks [18]–[21], or to infer transcription regulatory networks activated over time [12], [13], [22]. Though all the methods described so far are relatively successful in identifying network components and modules activated at specific time points, no attempt has been made to identify paths connecting genes expressed at different time points. Such temporal paths can show potential connections between genes expressed at different stages of a response thus providing information about intermediate, transiently expressed regulators that would otherwise have been overlooked. In this work, we studied the innate immune response in dendritic cells (DCs) stimulated by lipopolysaccharide (LPS). LPS is a component of the outer membrane of Gram-negative bacteria and specifically binds to the TLR4 receptor, triggering both the downstream MyD88 and TRIF-dependent pathways. We used time-course gene expression profiles collected at 8 time points after LPS stimulation in the context of a high-confidence PPI, PDI and post-translational modifications (PTM) network. We grouped the gene expression profiles into three groups – the initial response genes (greatest fold-change in expression between 0. 5–1 hour after stimulation), the intermediate regulators (greatest fold-change in expression between 2–4 hours after stimulation) and the late effectors (greatest fold-change in expression between 6–8 hours after stimulation). We then attempted to identify the most probable paths connecting the initial response genes to the late effectors in the interaction network, while taking into account the intermediate regulators. In order to do this, we used a network flow optimization approach allowing the flow to follow a time-dependent path within the molecular network. Using this method, we were able to identify an optimal gene sub-network for activated DCs. Based on this sub-network, we identified several known core components of the innate immune response, novel down-stream participants and pathways connecting these core components. We were able to identify genes playing an important role in the innate immune response but showing no observable change in expression. We also analyzed time-course gene expression profiles of MyD88-knockout cells and TRIF-knockout cells, and compared their gene sub-networks to that obtained for wild-type DCs in order to identify the components that are independently activated in each pathway. Finally, we identified the distinct functional classes of genes expressed during different stages of the immune response and how their patterns of expression change in MyD88 and TRIF-knockout DCs compared to those in wild-type DCs.
We used a minimum cost flow optimization approach to identify important components of the innate immune response over time on LPS stimulation. A network of PPI and regulatory interactions, including transcription factor-target gene, phosphorylation, dephosphorylation and ubiquitination relationships, was prepared. Network edges were scored based on interaction reliability as obtained from the protein-protein interaction database, HitPredict [23]. Time-course gene expression levels were obtained using RNA-seq from DCs before LPS stimulation and up to 8 hours after LPS stimulation. The genes with significant changes in expression after LPS stimulation were divided into 3 groups based on the time of their greatest change in expression: In order to identify potential paths through the molecular network connecting the genes within the three groups, we formulated the problem as a minimum cost flow optimization problem incorporating the gene expression levels in three stages. Figure 1 shows a schematic representation of the proposed method. We set our source nodes as the initial response genes. The target nodes of the network were the late effector genes. Edges of the network were assigned costs that were inversely proportional to their interaction reliability. Edges were also given a flow capacity proportional to the observed change in expression of the adjacent genes. A constraint was added to the flow optimization problem to force the flow to go through at least one intermediate regulator. We solved the optimization problem to identify the path of minimum cost for the flow to pass through the network using linear programming techniques (see Materials and methods for the problem formulation). The method found the most probable paths in the network between genes expressed in the initial response and those expressed at a later time while taking into account the genes expressed during the intermediate stage. Each edge of the optimal sub-network was assigned a flow signifying its importance. This resulted in a weighted gene sub-network where the edges were scored according to their importance. Flows were calculated for nodes, or genes, as the sum of the flows of their incoming edges. Genes with high flows were considered important due to their connection to high-flow edges. The reliability of the optimal solution was confirmed and statistical significance was calculated for each gene in the optimal sub-network by randomizing the source and target nodes (see Materials and methods). The flow assigned to a gene within the sub-network shared an inverse relationship with its statistical significance, demonstrating that a high flow was a good indicator of reliability (Figure S1). The genes with the highest flows – Socs3, Nfκb1, Jak2, Jun, Fos, Cxcl10 and Stat1 are well-known components of the innate immune response. Table 1 shows 20 genes with the highest flows in the optimal sub-network for activated wild-type DCs (See Table S1 for the list of all predicted genes and their statistical significance). As shown by the results, the method not only predicted essential genes expressed within each of the 3 groups, but also genes for which no significant change in expression was detected but were connected to others with significant changes in expression over time. In order to evaluate the reliability of the gene network resulting from the paths identified by solving the flow optimization problem, we compared the genes in the optimal sub-network with the experimentally identified regulators of the innate immune response from previous perturbation experiments [3], [4]. Of the 125 regulators identified by Amit et al. [3], our sub-network contained 62 (49. 6%), all of which had a flow greater than 1 (Table S2). In a similar study by Chevrier et al. [4], our sub-network contained 30 of the 43 known or novel regulators identified (69. 8%), and 56 of the 102 (54. 9%) TLR target genes affected by the perturbation of these regulators (Table S3). The sub-network also contained the gene, Polo-like kinase 2 (Plk2), which activates a distinct signaling cascade. Thus, our sub-network contained a significant number of the regulators of the innate immune response that were recently experimentally identified. We further confirmed the quality of the predicted gene network through Gene Ontology (GO) and KEGG pathway enrichment analysis. The genes having flows greater than 1 in the sub-network, were enriched for the Toll-like receptor signaling pathway (p = 5. 10e-41), Jak-STAT signaling pathway (p = 4. 88e-45), pathways in cancer (p = 2. 50e-41) and chemokine signaling pathway (p = 5. 16e-40) among others (See Table S4 for full list). The association of the predicted genes with the innate immune response is further confirmed by the GO Biological Process terms enriched for these genes. Protein amino acid phosphorylation (p = 7. 80e-36), immune response (p = 1. 35e-32) and regulation of programmed cell death (p = 1. 72e-29) were some of the most enriched terms (See Table S5 for full list). 49. 7% of the genes identified in the optimal sub-network did not show significant change in their expression levels on LPS stimulation. In order to confirm that these genes contribute to the enrichment of functional terms associated with the innate immune response, we compared the enrichment of the KEGG pathways and the GO terms in all predicted genes with those that showed differential expression after LPS stimulation (Table S6, S7). Including predicted genes lacking differential expression significantly improved the enrichment of the KEGG pathways and the GO terms associated with the innate immune response over that observed for differentially expressed genes only. This further confirmed the association of the genes predicted in optimal sub-network with the innate immune response. Additional analysis of GO term enrichment of genes identified in the sub-network at each time point showed the distinct processes active during different stages of the immune response. Table 2 shows the most significant GO Molecular Function and Cellular Component terms enriched in genes identified at each time point. The most significant term enriched for genes expressed between 0. 5–1 hour is “transcription regulator activity” (p = 1. 18e-09) for 20% of the genes indicating an upregulation of transcription factors during the first hour of the immune response. On the other hand, genes predicted at 2–4 hours are enriched for “nucleotide binding” (p = 9. 33e-04,28. 5% genes) and “protein kinase activity” (p = 1. 27e-03,13% genes) suggesting a role for signal transducers. Finally, the terms enriched for genes predicted between 6–8 hours are “proteasome complex” (p = 2. 98e-11,7%) and “peptidase activity” (p = 5. 2e-08,13%) highlighting the activity of the immunoproteasome during this phase of the innate immune response. Finally, genes that were identified in the optimal sub-network but which did not show change in expression during the sampled time points were enriched for GO terms such “protein kinase activity” (p = 7. 52e-31,16%), “cytokine binding” (p = 5. 9e-26,6%) and “transcription factor activity” (p = 1. 18e-07,12%) (Table 3). To check the quality of the network paths predicted by the method, we identified all the possible paths predicted in the optimal sub-network that matched a directed path of the same length in a KEGG pathway. Our method was able to predict directed paths of 3 edges or more in 13 KEGG pathways, including the Jak-STAT signaling pathway, the Chemokine signaling pathway, the Toll-like receptor pathway and the MAPK signaling pathway (Table 4, Table S8). The longest predicted directed path contained 7 edges and was part of the Jak-STAT signaling pathway. Thus, the method was able to partially recover known pathways in the form of short paths connecting genes expressed at consecutive time points. We also identified all shortest paths up to 3 edges (i. e. containing 4 nodes at most) between genes expressed at different stages of the immune response and checked how well they were represented in the same KEGG pathway. We found that 84. 9% of the predicted paths have at least 2 genes in the same KEGG pathway, while 11. 6% of the paths have all genes in the same KEGG pathway (Figure S2). Taken together, these results confirm the reliability of the optimal gene sub-network identified for activated wild-type DCs. To demonstrate the utility of our algorithm, we compared the optimal sub-network identified by our method to that identified using a non-temporal minimum cost flow optimization method, ResponseNet [15]. Using minimum cost flow optimization through our initial network, ResponseNet identified paths from the initial response genes to the late effectors without taking the intermediate regulators into account (Table S9). Table 5 shows the results of the comparison between the optimal sub-networks predicted by our method and ResponseNet. ResponseNet identified fewer genes and interactions in the predicted sub-network. More significantly, since there was no constraint for the flow to pass through the intermediate regulators, it identified only 49 of these as compared to the 154 by the current method. Our method also identified significantly higher number of known regulators in the innate immune response in addition to longer paths in associated pathways. On the other hand, ResponseNet failed to identify a directed path of 3 or more edges within any KEGG pathway associated with the innate immune response. These results clearly demonstrate that including the intermediate regulators into the problem formulation, as we propose here, improves the ability of the method to predict candidate genes and associated networks using time-course gene expression profiles. The gene predicted with the highest flow in the optimal sub-network was Suppressor of cytokine signaling 3 (Socs3) followed by Nuclear factor κb1 (Nfκb1). Both genes were significantly upregulated between 2–4 hours and are well-known regulators of the innate immune response. Socs3, along with Socs1 and Socs2, is an inhibitor of cytokine signaling pathways. It is a key regulator of interleukins 6 and 10 (Il6 and Il10) [24]. In the identified sub-network, Socs3 is induced by the primary regulators of the immune response such as Nfκb1 and inhibits a large number of proteins, specifically interleukin receptors (Figure 2a). Nfκb1 is induced both in the early and late phase of the innate immune response and is primarily responsible for the expression of inflammatory cytokines. Other genes identified with high flows were the Janus kinase 2 (Jak2), Rous sarcoma oncogene (Src) and phosphoinositide-3-kinase, regulatory subunit 5 (Pik3r5), all of which have been implicated in the TLR response pathway. Src, a protein tyrosine kinase that modulates a large number of signaling pathways during the innate immune response, was upregulated between 2–4 hours. Along with Src, other tyrosine kinases from the Src family, such as Hck and Lyn, were also identified (Figure 2b). Syk, another protein tyrosine kinase of the Syk-ZAP70 family that is found in innate immune cell types, was also identified as part of the network though no significant change in gene expression levels was detected at the tested time points (Figure 2c). Several other components of the Src signaling pathways like Card9, Cblb, Fcerγ and various integrins were also identified within the gene sub-network. Among other known regulators, the induction of Ralgds by Ras proteins, and the further upregulation of the Rac genes, was also detected (Figure 2d). Gadd45b, an anti-apoptotic inhibitor induced by Nfκb [25] was also part of the sub-network. Gadd45b was significantly upregulated between 2–4 hours and was predicted to inhibit the cyclins B2, B3 and CDK (Figure 2e). Another anti-apoptotic inhibitor, the X-linked inhibitor of apoptosis (Xiap) was also identified. Xiap is regulated by Nfκb and in turn inhibits Casp3 and Casp7 thus controlling apoptosis (Figure 2f) [25]. Another class of proteins identified, were the Akt serine-threonine protein kinases Akt1, Akt2 and Akt3, which are downstream effectors of the PI3K pathway (Figure 2g). Expression level change was only observed for Akt1 which was down-regulated at 0. 5–1 hours followed by an up-regulation at 3 hours. Other predicted components include the Dual specificity phosphatases (DUSP proteins) which were significantly upregulated between 0. 5–1 hour, except Dusp6. The Dusp proteins regulate the immune response by dephosphorylating the Map kinases and repressing the LPS-induced inflammatory response (Figure 2h). Interestingly, the network indicated that the Dusp genes were expressed within the early stages of the innate immune response suggesting that control of inflammation begins soon after its induction. Many of the genes identified in the network do not show any significant change in expression after activation of the DCs, but are known to be essential for the response. An example is the protein phosphatase 2a catalytic subunit α (ppp2ca) which has a high flow in the sub-network. A serine threonine phosphatase required for the dephosphorylation of the 20S proteasome subunits, ppp2ca is known to affect the ability of the proteasome to degrade substrates, along with protein kinase A (PKA) [26]. Ppp2ca has also been recently shown to play an important role in the regulation of endotoxin tolerance through the regulation of MyD88 activity [27]. The identified gene sub-network indicated extensive interactions between ppp2ca and the subunits of the immunoproteasome, suggesting a role of ppp2ca in the regulation of the immunoproteasome (Figure 3). The immunoproteasome is induced by interferons and is central to the regulation of the immune response and in the prevention of auto-inflammatory diseases through its ability to degrade toxic protein aggregates during cytokine-induced oxidative stress [28]. We applied the method described above to time-course gene expression profiles obtained from DCs of MyD88 and TRIF-knockout mice in the context of the comprehensive molecular interaction network. MyD88 and TRIF are essential components of the innate immune response and trigger distinct pathways that result in the activation of early and late phase Nfκb, respectively. Previous studies have shown that Nfκb and Mapk8 (JNK) are activated in a delayed manner in MyD88-knockout cells. However, inflammatory cytokines like IL12 or TNFα are not produced [29]. In order to identify the MyD88-independent response network, we used gene expression levels from MyD88-knockout DCs to assign edge capacities, and removed the MyD88 gene and its links within the network prior to solving the minimum cost flow optimization problem (See Table S10 for identified genes and edges). We performed a similar analysis on the data from TRIF-knockout DCs by removing TRIF and its links from the network and predicting a MyD88-dependent response network on LPS stimulation (See Table S11 for identified genes and edges). A comparison of the genes and their flows in the identified sub-networks suggests that the response pathways active in the wild-type and TRIF-knockout sample are similar (Figure 4a). The active sub-networks identified for both these samples are enriched in the KEGG pathways “Cytokine-cytokine receptor interaction” (p = 1. 13e-29), “Jak-STAT signaling pathway” (p = 2. 34e-15) and “Toll-like receptor signaling pathway” (p = 6. 07e-11). These findings suggest the dominance of the MyD88 pathway in the wild-type response. Indeed, this dominance has been previously observed during pulmonary infection [30]. On the other hand, the most enriched pathways in the genes exclusively identified in the MyD88-knockout network are the “Circadian rhythm” (p = 6. 29e-5) and “Ubiquitin mediated proteolysis” (p = 3. 2e-4) suggesting an association between these pathways and the MyD88-independent, TRIF-dependent pathway (Table 6, Tables S12 and S13). In order to identify the dominant changes in the immune response over time, we classified the genes from the optimal sub-networks obtained for the wild-type, MyD88-knockout and TRIF-knockout DCs into functional classes. Global changes in the expression patterns of genes identified as part of the optimal sub-network at each of the 3 stages showed a dominance of functionally distinct groups at different times during the immune response (Figure 4b). In wild-type DCs, transcription factors and enzyme modulators were predominantly expressed during 0. 5–1 hour after LPS stimulation. On the other hand, kinases and signaling molecules were abundant between 2–4 hours after stimulation. Finally, proteases and defence/immunity proteins along with receptors showed the greatest changes in expression in the late phase of the immune response between 6–8 hours. TRIF-knockout DCs showed similar changes in the expression patterns of genes. However, these patterns were significantly different in the MyD88-knockout DCs. Transcription factors were not as significantly upregulated in the early phase, but more so in the late phase, when the expression of proteases and defence/immunity genes was significantly reduced. Thus, the identified sub-networks suggest a pattern in the global change in gene expression during the different stages of the immune response. The similarity of the patterns of gene expression in the TRIF-knockout DCs and wild-type DCs further support the dominant role of the MyD88-dependent pathway in the innate immune response. An analysis of the functional distribution of the genes predicted in the network, but not showing significant differential expression on activation, illustrates their similarity to the intermediate regulators in the wild-type as well as knockout DCs. Several important components of the innate immune response were identified in both knockout sub-networks, however, with significantly different flows. Nfκb1, Jak2 and Socs1 were genes with the highest flows (>40) in the TRIF-knockout network. These genes were also identified in the MyD88-knockout network, but with flows just above 1. This disparity in the flows possibly indicates their changing levels of expression and significance within the two sub-networks. The sub-network associated with MyD88-knockout DCs had different genes with high flows – Akt3, Casp8 and Stat2. Interestingly, the kinase Pik3r5 had similar levels of predicted flow in both knockout networks. It was upregulated in both instances but much more so in the MyD88-knockout DCs. Git1 and Cry1 were two of the important candidates identified only in the MyD88-knockout gene network. Git1 (G-protein coupled receptor kinase interacting protein 1) acts in the formation of a scaffold to bring together molecules to form signaling modules and increase the speed of cell migration. Its role in the innate immune response is currently not known. However, it was significantly upregulated in the MyD88-knockout sample and found to interact with Pxn, Arhgef6 and Arhgef7 (Figure 5a). The other important gene identified, Cry1, is a key component of the circadian core oscillator complex. The role of Cry1 in the negative regulation of the activation of Nfκb and further induction of proinflammatory cytokines has been recently elucidated [31]. Cry1 was significantly upregulated in the MyD88-knockout DCs between 6–8 hours after stimulation and could potentially be regulating the activation of Nfκb signaling. Though Cry1 was part of the gene network associated with the activation of wild-type DCs, it was not identified in the optimal gene sub-network associated with TRIF-knockout DCs, suggesting that the upregulation of Cry1 and its role might be controlled by the MyD88-independent, TRIF-dependent pathway (Figure 5b). The MyD88-knockout associated gene network also contained a number of genes from the E2 and E3 ubiquitin-conjugating enzyme families, including several members of the Trim family, which are known for their role in suppressing the immune response by increasing the ubiquitination and subsequent degradation of regulatory genes [32]. The selective prediction of these ligases in the MyD88-knockout response network suggests that proteolytic degradation might also be predominantly affected by the TRIF-dependent pathway. The response network identified for the TRIF-knockout sample highlights the wild-type MyD88 pathway wherein MyD88 triggers the activation of Nfκb which in turn induces the inflammatory cytokines, further inducing the Jaks and Stats and finally upregulating the Socs genes which repress the inflammatory response (Figure 5c).
The innate immune response is complex and occurs through multiple pathways. The interplay within the activated pathways makes the identification of novel components and their associations difficult. In this study, we addressed this issue by using time-course gene expression profiles of activated dendritic cells in combination with a comprehensive molecular interaction network. We developed a method based on minimum cost flow optimization in a large interaction network to identify paths between genes expressed at different time points of the immune response. Using this method, we identified an optimal gene sub-network activated during the innate immune response. We confirmed the role of several known and novel components in the identified network and suggest a role for the protein ppp2ca in the regulation of the immunoproteasome. A flow value was assigned to each identified gene and interaction within the network indicative of its importance. We also compared the response of the wild-type DCs with DCs from MyD88-knockout mice and TRIF-knockout mice and identified the global changes in expression patterns of genes in distinct functional classes. Our results are consistent with previous studies suggesting the dominant role of the MyD88-dependent pathway. We further showed that genes related to proteasomal degradation and circadian rhythms are primarily associated with the MyD88-independent, TRIF-dependent pathway. The method proposed here is independent of the biological system and can be used to identify time-dependent gene sub-networks with the help of time-course gene expression profiles related to any other cellular conditions. Future work in this area will be aimed at developing methods to accurately predict longer pathways while incorporating time-course gene expression profiles from multiple time points without the necessity of grouping them.
GM-CSF-induced bone marrow-derived dendritic cells (GM-DCs) were prepared from C57BL6/J mice (purchased from Japan Clea Inc.) as described previously [3]. The cells were stimulated with LPS from Salmonella minessota Re-595 (purchased from Sigma) at a concentration of 100 ng/ml. Stimulated cells were harvested at 0,0. 5,1, 2,3, 4,6, 8,16,24 hours after stimulation. Total RNA was extracted from the cells using TRIzol (Invitrogen) according to the manufacturer' s instruction. The RNA was subjected to RNA-seq as described in a previous study [35]. Mice deficient in MyD88 or TRIF were prepared as described in an earlier study [36], [37]. The RNA-seq data is available in the Sequence Read Archive under the accession number DRA001131. The RNA-Seq data for the wild type, MyD88-knockout and TRIF-knockout DCs at 10 time points, from 0 to 24 hours after LPS stimulation, were obtained in the form of 35 bp single-end reads. The reads were mapped to the RefSeq mm9 mouse reference genome [38] using Bowtie [39]. Exon-exon junctions were found using TopHat [40] with each read having at most 2 mismatches and 20 mappings to the reference genome, and a minimum intron length of 70 bp. For each read, the mapping with the highest alignment score was selected. The mapping statistics are shown in Table S14. Transcript abundances for all three samples at 10 time points were estimated using Cufflinks and Cuffdiff [41] using the –T option to treat the samples as a time-series. The data from the last two time points, 16 hours and 24 hours, was not used in this study because we were concerned about the effect of their large separation from the prior time points on the quality of the sub-network predicted. Maximum absolute log fold change in expression was calculated for each gene over all time points, as follows: (1) Where = maximum absolute log fold change for gene i over time j where j = {0. 5,1, 2,3, 4,6, 8}, Genes with at least 2 fpkm in 50% of the experiments, at least 10 fpkm for 2 or more time points and an absolute fold change greater than 2 for at least one time point in each sample were considered for further analysis. Each selected gene was assigned to one of the following groups depending on the time at which it showed the highest absolute fold change (Figure 1A): The genes and their expression levels are shown in Tables S15, S16, S17. A network of regulatory and physical interactions from mouse was prepared by combining the following datasets: Table S18 shows the counts of the different interaction types included in the network. PPIs were considered as bi-directional edges whereas all other associations (transcription factor–target gene, functional association, expression regulation, post-translational modification and inhibition) were considered uni-directional. Genes and their corresponding proteins were represented by a single node in the network. The edges of the network were weighted according to their reliability. Reliability scores provided by HitPredict and TRANSFAC were used. Innatedb core PPIs and interactions from KEGG pathways were uniformly assigned a high reliability score of 999 since these were manually curated. All scores were scaled to values between 0 and 0. 8 as shown in Table S19. (2) Where = scaling functionThe complete network of 103218 interactions among 12856 proteins, or protein complexes, along with the data source, reliability scores and edge weights is given in Table S20. The network was denoted by a graph G = (V, E) with E edges and V nodes (including the auxiliary source S and the auxiliary sink T). The auxiliary source, S, was connected to the set of initial response genes (GT1), while the auxiliary sink, T, was connected to the late effector genes (GT3). Direct edges between GT1 and GT3 were excluded. The intermediate regulators (GT2) were also a part of the network but not connected to the S or T nodes. All edges, E, were assigned a capacity and a cost (See Figure 1B). The genes identified as part of the optimal gene sub-network were assigned a statistical significance. This was done by solving the minimum cost flow optimization problem 5000 times using the original molecular interaction network but with randomly selected genes in the GT1, GT2 and GT3 sets i. e. initial response genes, intermediate regulators and late effectors in numbers equal to those from the real sample. The p-value was calculated as the fraction of solutions in which a gene was identified with an equal or higher flow than that in the optimal sub-network and with at least all the connecting edges in the optimal sub-network. We observed that the predicted flow in the final network increased with decreasing p-value (Figure S1) suggesting that a high flow was a good indicator of high statistical significance and hence greater reliability. The ResponseNet algorithm was implemented as a non-temporal minimum cost flow optimization method. The problem formulation was changed to remove the constraint in equation (12) thus allowing the flow to go from the initial response genes (source nodes) to the late effectors (target nodes) without being constrained to pass the intermediate regulators. Additionally, the term involving was also removed from the optimization problem. The algorithm was run on the same network as our method with identical edge costs. The capacities of edges GT1 – S and GT3 – T were set as described in equations (3) and (6). The capacity of all other edges was set to 1. The optimal solution was calculated for and the identified genes were compared to known regulators and KEGG pathways as described in the Results. All possible paths within the optimal sub-network from the initial response genes (GT1) to the late effectors (GT3) were identified and compared to all KEGG pathways to determine their overlap. Paths were predicted between genes in the groups GT1, GT2 and GT3 by finding the weighted shortest paths [47] of up to 3 edges in the optimal sub-network. The edges were weighted as per the formula suggested by Opshal et al. [48]: (13) where = flow assigned to edge (i, j), The shortest weighted paths identified were then compared to all KEGG pathways. An optimal gene sub-network was identified by solving the flow optimization problem using the time-course genes expression profiles from MyD88 and TRIF-knockout DCs in a manner similar to that described above for the wild-type DCs. The MyD88 gene and its interactions were removed from the starting network when the MyD88-knockout gene expression levels were considered. Similarly, during the analysis of the TRIF-knockout sample, TRIF (Ticam1) and its interactions were removed from the network. Enriched Gene Ontology terms and KEGG pathways were obtained using DAVID [49] with all mouse genes used as the background. Networks were prepared and formatted using Cytoscape2. 7 [50]. Protein functional classes were identified using PANTHER [51].
The innate immune response is the first level of protection in organisms against invading pathogens. It consists of a large number of proteins functioning in signaling cascades triggered by the binding of fragments from microbes to specific cellular receptors. Disruptions in these pathways can lead to numerous diseases. As such, the innate immune system has been the subject of a large number of studies. However, due to its complexity and the interplay of a large number of pathways, it is not yet completely understood. In this study, we measured transcriptional changes in activated immune cells and used this information in the context of a large network of protein-protein and protein-DNA interactions to identify a smaller network of response genes. We did this by identifying the most probable network paths connecting genes showing large changes in their expression patterns at successive stages of the response. Analysis of this activated gene network revealed the associations between various temporal regulators of the innate immune response. We also identified response networks for immune cells lacking important mediators, MyD88 and TRIF, to clarify the distinct functional modules affected by their associated pathways in the innate immune response.
Linking Transcriptional Changes over Time in Stimulated Dendritic Cells to Identify Gene Networks Activated during the Innate Immune Response