Abstract Toll-like receptors (TLRs) are critical components to stimulate immune responses against various infections. Recently, TLR agonists have emerged as a promising way to activate anti-tumor immunity. L-pampo, a TLR1/2 and TLR3 agonist, induces humoral and cellular immune responses and also causes cancer cell death. In this study, we investigated the L-pampo-induced signals and delineated their interactions with molecular signaling pathways using RNA-seq in immune cells and colon and prostate cancer cells. We first constructed a template network with differentially expressed genes and influential genes from network propagation using the weighted gene co-expression network analysis. Next, we obtained perturbed modules using the above method and extracted core submodules from them by conducting Walktrap. Finally, we reconstructed the subnetworks of major molecular signals utilizing a shortest path-finding algorithm, TOPAS. Our analysis suggests that TLR signaling activated by L-pampo is transmitted to oxidative phosphorylation (OXPHOS) with reactive oxygen species (ROS) through PI3K-AKT and JAK-STAT only in immune and prostate cancer cells that highly express TLRs. This signal flow may further sensitize prostate cancer to L-pampo due to its high basal expression level of OXPHOS and ROS. Our computational approaches can be applied for inferring underlying molecular mechanisms from complex gene expression profiles. Subject terms: Computational biology and bioinformatics, Gene regulatory networks, Bioinformatics, Pattern recognition receptors Introduction Toll-like receptors (TLRs) are pattern recognition receptors originally discovered to stimulate innate immune reactions against microbial and viral infection. TLRs are essential in bridging the innate and adaptive immune system and have multifaceted roles in immunology, inflammation, autoimmune diseases, and cancers^[36]1–[37]3. Almost all innate immune cells, such as macrophages, dendritic cells, neutrophils, and non-immune cells, including epithelial cells and cancer cells, express a variety of TLRs. The engagement of TLRs and their ligands triggers the activation of downstream signaling pathways through myeloid differentiation factor 88 (MyD88) and toll-IL-1 receptor structural domain (TRIF) that leads to the activation of NF-κB and IRF7 transcription factors^[38]4. These signal cascades initiate the production of pro-inflammatory cytokines, including tumor necrosis factor-α (TNF-α), IL-1 and IL-6, and type I IFNs^[39]4. Due to the robust immune stimulation of TLRs, various TLR-targeted strategies that boost the immune system to kill cancer cells have emerged in cancer treatment^[40]5. Currently, many preclinical and clinical studies provide the potential of TLR agonists as anti-tumor agents that regulate tumor microenvironment and induce tumor clearance^[41]6,[42]7. We have developed a proprietary adjuvant called L-pampo, composed of TLR1/2 and TLR3 agonists. L-pampo activates innate immune cells through TLR1/2 and TLR3 signal pathways and induces robust humoral and cellular immune responses^[43]8. Recently, we have discovered that L-pampo has immunotherapeutic properties that promote anti-tumor immunity. L-pampo activates dendritic cells and M1 macrophages critical for CD8 T cell expansion and infiltration. Interestingly, L-pampo can induce cancer cell death through alarmin signaling, such as HMGB-1, in cancer cells, promoting an antigen-spreading process that bolsters the innate immune cell activation and anti-tumor immune responses in the tumor microenvironment^[44]9. However, how L-pampo leads to cancer cell death by regulating complex TLR downstream of signal pathways in the cancer cells, how the L-pampo signal network is differentially involved in both cancer cells and immune cells and which signaling molecules are critical in intervening between TLR signaling and cancer-related signaling are not fully understood yet. In this study, we investigated the effect of L-pampo on the TLR signaling pathways in innate immune cells and cancer cells, then examined the interaction of TLR-related signaling with cancer-related signal pathways. To this end, we designed a network-based analysis framework to demonstrate the differences in transcriptional expression profiles among cell types (Fig. [45]1). We performed RNA sequencing (RNA-seq) on three different human cell lines based on TLR2/3 expression and the responsiveness to L-pampo treatment: a differentiated human monocyte cell line (THP-1), a prostate cancer cell line (PC-3), and a colon cancer cell line (SW620) at different time points (0h, 3h, and 6h). In Step 1, we identified differentially expressed genes (DEGs) and the influential genes using network propagation^[46]10. In Step 2, we identified the perturbed gene modules potentially affected by L-pampo using the weighted gene co-expression network analysis (WGCNA)^[47]11. Then, we constructed a template network to compare the cellular differences between immune cells and cancer cells in a more unbiased way. In Step 3, we detected the core subnetworks of TLR signaling genes and the major perturbed gene modules for signal reconstruction from L-pampo treatment to cell death. We also quantified how the cell lines differ based on the subnetworks in terms of differential responses to L-pampo. Consequently, we could summarize the signaling pathways found by the above steps into several key molecular functions. Figure 1. [48]Figure 1 [49]Open in a new tab Overview of comparative network-based analysis. In step 1, RNA-seq data of three cell lines (THP-1, SW620, and PC-3) with L-pampo treatment were used as input to identify differentially expressed genes (DEGs) and influential genes in each cell line. In step 2, the global template network was constructed from the gene co-expression profiles across the three cell lines. Additionally, perturbed modules with L-pampo were found using WGCNA. In step 3, subnetworks connected to the perturbed modules were detected and the relative differences between these subnetworks were quantitatively measured. Finally, key molecular functions inferred by this process were further researched. Results Identification of DEGs and influential genes perturbed by L-pampo (Step 1) Quantifying gene expression Our main objective in this study is to determine whether L-pampo drives the differential responses of cellular networks in different cell types. Particularly, to investigate gene signatures associated with L-pampo in immune cells and cancer cells, we chose a differentiated human monocyte cell line, THP-1, that highly expresses TLR2/3 and induces various immune responses to L-pampo^[50]12,[51]13. Additionally, we utilized two cancer cell lines, PC-3 and SW620, that can express TLR2/3^[52]14,[53]15 but manifest differential responses by L-pampo treatment. The cell death was induced to PC-3 in a dose-dependent manner, while L-pampo did not affect the cell death on SW620 (Figure S1a-b). Given that cancer cell death is an essential cellular process in the tumor microenvironment that leads to enhancing anti-tumor immune responses, we sought to examine transcriptomes related to cell death induced by L-pampo. We treated cells with L-pampo and isolated RNA samples at different time points (0, 3h, and 6h), and performed the RNA-seq to obtained gene expression profiles. As expected, the TLR2/3 are more highly expressed in immune cell (THP-1) compared to the cancer cells (Figure S1c). Also, we observed that transcriptomes of TLR2/3 were more highly expressed in PC-3 than SW620, suggesting that the signal transmission from the TLR receptors may depend on TLR2/3 expression in both immune cells and cancer cells (Figure S1c). Differentially expressed gene (DEG) and functional analyses First, we analyzed differentially expressed genes (DEGs) in three cell lines (THP-1, SW620, and PC-3) and three time points (0h, 3h, and 6h after treating L-pampo). The analyses revealed that 3,732 genes, 501 genes, and 1,415 genes upregulated or downregulated in THP-1, SW620 and PC-3 cell lines across different time points, respectively (Figure S2a). To understand the biological roles of the putative genes induced by L-pampo, we performed over-representation analysis for each DEG at THP-1, SW620, and PC-3 using KEGG pathway database^[54]16. Only 20 DEGs were identified among three cell lines and related to toll-like receptor (TLR) signaling pathway and cytokine-cytokine receptor interaction pathway. The 558 genes in the intersection between THP-1 and PC-3 were related with ribosome and oxidative phosphorylation (OXPHOS) (Figure S2b). The over-representation analysis showed that the common DEGs from 0 to 6h in THP-1 and PC-3 were related to JAK-STAT signaling and cytokine-cytokine receptor interaction besides ribosome and OXHOS, while no significant changes were observed in SW620 (Fig. [55]2). Figure 2. [56]Figure 2 [57]Open in a new tab KEGG pathway analysis for DEGs from THP-1, SW620, and PC-3. The columns represent cell lines and time changes, and the rows show the selected significant KEGG pathway terms from the cell lines. In each column, the length of each bar indicates the number of DEGs contained by the respective KEGG pathway. The bar color shows the significance of pathway enrichment on the -log10 scale for the adjusted p-value. The bold characters indicate the KEGG pathways of our interest perturbed by L-pampo in three cell lines. The gene signatures related to reactive oxygen species (ROS), NOD-like receptor signaling, IL-17 signaling, TGF-beta signaling, and lipid and atherosclerosis were also upregulated in PC-3 between 0 and 6h. Interestingly, RIG-1, a pathogen pattern recognition receptor, was also found to be significantly upregulated in SW620 and PC-3 compared to THP-1 (Figure S2c). On the other hand, the DEGs related with necroptosis which is a programmed inflammatory cell death^[58]17 were increased between 3 and 6h in PC-3 (Figure S2d). As an example, NLRP3, a well-known key regulator in necroptosis varied in the same way as necroptosis genes^[59]18. Influential gene analysis by network propagation DEG analysis is very powerful since it led us to a probably list of perturbed genes by L-pampo treatment across cell types and time points. However, those genes/proteins that play critical roles given a molecular stimulus (e.g., diseases) are often not caught by only DEG analysis. In many cases, DEGs may reflect significant transcriptomic changes as an effect of perturbation not a cause. On the other hand, conventional pathway enrichment analysis based on DEGs often yields a large number of pathways, making it challenging to organize the results and draw clear conclusions from fragmented evidence. Moreover, our DEG analysis only showed 20 DEGs common to the three cell lines in the 0h-6h comparison (Figure S2b), which account for just 0.5% of the total DEGs. This implies that hidden regulators may be responsible for the differences among the cell lines. To tackle this challenge, adopting gene–gene interaction networks is a potent strategy to augment DEG analysis and this approach has been used to identify functional genetic modules and even drug targets^[60]19. Fundamentally, using gene–gene interaction networks can allow us to search the vicinity of DEGs for hidden players that may cause or mediate important molecular signals. Therefore, we devised a network-based analysis and investigated the complex gene–gene interactions in immune cells and cancer cells to understand the L-pampo induced mechanism. To identify hidden players, we utilized network propagation that prioritizes genes by iteratively considering the influence of the seed genes (i.e., DEGs) on their neighbor genes in a protein–protein interaction (PPI) network^[61]19. We ranked the genes that may highly affect DEGs, naming the top-ranked genes as “influential genes” and this method as “influential gene analysis”. A PPI network is a collection of all known and curated interactions between genes/proteins. Analyzing the connections of DEGs and transcription factors on a PPI network can provide insights into key regulators at the transcriptome level^[62]10. However, all the interactions/correlations reported within a PPI network are not functionally meaningful to our specific problem solving, as a PPI network contains any known biological interactions accumulated to date. Thus, we also computed gene–gene correlations from our data and prune those edges (i.e., interactions reported in the PPI network) that did not have high correlation coefficients (Pearson correlation coefficient < 0.7)^[63]20 . As a result, we constructed our own gene–gene interaction network reflecting cellular specificity of our data (THP-1, SW620, and PC-3) and we then identified the top 200 influential genes using DEGs as seeds on the constructed network for the following steps. Instantiation of a template network for comparison (Step 2) Construction of a template network Although consideration of temporal dynamics and cell line characteristics is essential for our study design, the DEG and influential gene analyses are constrained in multifactorial comparisons due to distinct baselines for each condition. To address this limitation, we built and used a universal gene–gene interaction network to compare transcriptomic profiles of all different cells and time points on the same baseline and referred to this as “template network”. We utilized DEGs, influential genes, and transcription factors to construct the template network so that it can reflect molecular signaling flows induced by L-pampo on the three cell lines. In this network, we particularly added transcription factors from CollecTRI database^[64]20 to further consider transcriptional regulations over DEGs and influential genes. As a result, 7,252 genes were selected and used as input for WGCNA’s hard thresholding. In contrast to WGCNA’s soft thresholding, hard thresholding applies a cutoff to remove all weak correlations, simplifying the network structure and making it less sensitive to outliers. When we chose an optimal cutoff of 0.9 based on the scale-free topological fit, we obtained the final template network with 6,585 nodes and 1,214,811 edges. Identification of perturbed modules with L-pampo Next, we constructed hierarchical clustering using WGCNA on 4,582 DEGs identified across all conditions, showing the results using a dendrogram (Figure S3a). WGCNA requires a power parameter for a soft threshold because it assumes that increasing the power can reduce noise in the gene–gene correlation matrix. We chose a power of 18, which achieves a scale-free topological fit by R^2 larger than 0.9 (Figure S3b). To understand the relationship between gene modules and traits (i.e., cell type and time point), we performed module-trait relationship analysis, which informed the modules change in different directions depending on the cell type or time (Figure S3c). Different co-expression patterns over cell lines and time points guided us to select gene modules differentially perturbed by L-pampo. Among the ten identified modules, we finally chose five ‘perturbed modules’, named as A, B, C, D, and E, respectively. The correlation of Module A went up at 3h and stayed high until 6h only in THP-1. Module B had high correlation values in PC-3 compared to the other two cell lines. The correlation of Module C also increased at later time points in all cell lines, but it moves up faster in PC-3 than THP-1 (i.e., 3h vs. 6h). On the other hand, those of Modules D and E clearly decreased by time in THP-1 but not in the other two cancer cell lines. Next, we conducted a gene set enrichment analysis on the perturbed modules and found that there was an increase in the levels of PI3K-AKT signaling (Modules A and B) and oxidative phosphorylation (OXPHOS) (Modules B and C) (Figure S3d), respectively. To identify finer expression patterns within Modules A and B, we conducted a community detection analysis using the Walktrap algorithm to divide each module into submodules. This algorithm clusters densely connected gene communities given a network by comparing the neighborhood genes through simulated random walks^[65]21. We identified three submodules within Module A and named them as Modules A-1, A-2, and A-3. We repeated the same process for Module B, identifying two submodules and marking them as Modules B-1 and B-2. Upon a deeper analysis of the submodules, it was found that Module A-1 contained genes related to JAK-STAT, OXPHOS, and TLR signaling pathways. In contrast, Modules A-2 and A-3 were more associated with RAP-1 and PI3K-AKT pathways, respectively. JAK-STAT is closely linked to TLRs and involved in various cellular processes such as immune response, inflammation, or carcinogenesis^[66]22. Module B has a significantly large number of genes (676 genes), and most of these genes are highly expressed in PC-3. Modules B-1 and B-2 were related to PI3K-AKT and OXPHOS, whereas Module C was uniquely characterized by OXPHOS. These findings demonstrate that PI3K-AKT, JAK-STAT, and OXPHOS pathways can play essential roles in L-pampo-mediated response in THP-1 and PC-3. Oxidative phosphorylation (OXPHOS) and reactive oxygen species (ROS) as a key process in the L-pampo response by comparing three cellular networks (Step3) To identify the key molecular signals triggered by L-pampo, we used the TOPAS algorithm^[67]23 to analyze the shortest paths from the TLR pathway genes to the perturbed modules within the template network. We obtained six subnetworks, namely TLR-Module A, TLR-Module B-1, TLR-Module B-2, TLR-Module C, TLR-Module D, and TLR-Module E, that summarize the molecular interactions among the TLR pathway genes and the perturbed modules (Figure S4a). In addition to the qualitative assessment, we also attempted to quantify the differences of all subnetworks between cell lines. We used network-based Jensen-Shannon divergence (nJSD) to estimate the distance between probability distributions of the networks of interest^[68]24. Interestingly, our analysis revealed that TLR-Module A, TLR-Module B, and TLR-Module C subnetworks have longer distances between THP-1 and SW620 than THP-1 and PC-3. We also observed that these subnetworks showed larger differences between THP-1 and PC-3 compared to the other subnetworks (Fig. [69]3a). Figure 3. [70]Figure 3 [71]Open in a new tab Comparison of detected subnetworks. (a) Quantifying the relative differences of cell line pairs in terms of subnetwork between TLR genes and selected gene modules. The x-axis represents the cell line pairs and the y-axis the average of nJSD for the subnetworks. The average of nJSD was obtained by taking the mean of nJSDs at different time points. (b–d) Visualizing the expression changes in the subnetwork connected from TLR signaling genes to Module C (TLR-Module C) for the three cell lines: (b) THP-1, (c) SW620, (d) PC-3. Each circle represents a gene, and it has three time points of 0h, 3h, and 6h (counterclockwise). The node color represents the normalized TPM (Transcripts Per Million) using z-score across different cell lines and time points per gene. The color of the outer ring means the module ID (Module A: red, Module B: brown, Module C: magenta, Module D: green, and Module E: black) and the nodes marked by asterisks are TLR signaling genes. Especially, TLR-Module C subnetwork exhibited different expression patterns across the cell lines, which reflect the cell line specificity against L-pampo treatment (Fig. [72]3b-d). THP-1 showed the highest expression levels for TLR signaling genes that are related to pro-inflammatory cytokines, chemokines and activation markers compared to SW620 and PC-3 (Fig. [73]3b, 5 o’clock directed clusters with asterisk mark). SW620 displayed reduced expression levels in TLR-related genes within the subnetwork. However, some TLR signaling genes such as MAPK8, MAPK12, MAPK13, MAP2K6, TAB1, and AKT2 were relatively highly expressed in SW620 compared to THP-1 and PC-3, but their expression levels did not appear to be affected by L-pampo treatment (Fig. [74]3c, 1 o’clock directed cluster with asterisk mark). PC-3 showed increased gene expression levels in all genes within Module C and some Module B genes including MDH1, LSM3, and PSMA2 (magenta and brown color outer rings, respectively) compared to THP-1 and SW620 (Fig. [75]3d). Interestingly, MDH1 plays an essential role in the conversion of malate to oxaloacetate during the tricarboxylic acid (TCA) cycle in the OXPHOS process^[76]25. In particular, we observed that PIK3CD and IL-12A in Module B-1 played a crucial role in connecting the above major subnetworks such as TLR-Module A, TLR-Module B-2, and TLR-Module C (Figure S4a). TLR-Module A is enriched mostly by PI3K-AKT pathway genes, and PIK3CD and IL-12A acted as a bridge connecting Modules A-1 and A-3 (Figure S4a). These two genes are involved in TLR, PI3K-AKT, and JAK-STAT signaling pathways^[77]16,[78]26. Genes related to JAK-STAT signaling, such as IL-6, CSF2, CSF3, and IL-7R in TLR-Module A-1 (Figure S4a), were connected to TLR-Module C through hub genes IL-6 and S100A6 (Fig. [79]3b-d). Interestingly, PIK3CD and IL-12A simultaneously worked as a bridge connecting not only TLR-Module B-1 and TLR-Module B-2 but also TLR-Module C and TLR-Module B-2 (Figure S4a). Therefore, TLR-Module C can be bridged with TLR-Module B through these two genes as a hub. Given that OXPHOS and ROS pathways were enriched in Modules A-1, B-2, and C (Figure S3d), we examined the expression changes of OXPHOS and ROS genes in each module over time. This showed that the average expression levels of OXPHOS and ROS genes in Modules A-1, B-2, and C notably increased in PC-3 (Fig. [80]4a). When collecting all the OXPHOS and ROS genes from the modules detected by WGCNA, genes belonging to Modules A-1, B-2, and C in PC-3 exhibited distinct expression patterns compared to THP-1 and SW620 (Figure S5, the genes marked by stars). More specifically, those genes in Module C showed substantially increased expression levels at the earlier time point (3h) in PC-3. However, their expression levels were low but continued to increase until the latest time point (6h) in THP-1 (Figure S5). Figure 4. [81]Figure 4 [82]Open in a new tab Reconstructing functional networks by combining different modules. (a) Average of RNA expression of OXPHOS and ROS genes of each module. Each point and error bar indicates mean ± standard deviation of the normalized-TPM using z-score across different cell lines and time points per gene. Line colors denote three different cell lines: THP-1 (red), SW620 (green), and PC-3 (blue). The number under each module title indicates the number of OXPHOS and ROS genes in the that module. (b) The reconstructed network of Modules A-1, B-2, and C with OXPHOS and ROS genes (PC-3 DEGs). Each node represents a gene, and it has three time points of 0h, 3h, and 6h (counterclockwise). The node color represents the normalized TPM using z-score across different cell lines and time points per gene. The color of the outer ring means the module ID (Module A: red, Module B: brown, and Module C: magenta). The nodes marked by asterisks are OXPHOS and ROS genes. To elucidate the OXPHOS and ROS process in differential L-pampo response, we conducted network analyses using TOPAS on the OXPHOS and ROS pathway genes with the combination of Modules A-1, B-2, and C (Fig. [83]4b). Our analysis revealed that genes related to OXPHOS and ROS, such as NADH dehydrogenase (NDUF) and ATP synthase, which form complexes of the mitochondrial respiratory chain^[84]27, were predominantly upregulated in Modules A-1, B-2 and C (Fig. [85]4b, asterisk marked). As shown in Fig. [86]4a, these genes had higher basal expression levels in PC-3 than other cell lines and further activation of these genes due to L-pampo may have led PC-3 to cell death (Figure S5). Collectively, our study revealed that L-pampo treatment triggered PI3K-AKT and JAK-STAT pathways via TLRs and likely led to cell death through OXPHOS and ROS pathways. We found that PIK3CD and IL-12 in Module B-1 connected PI3K-AKT to OXPHOS and ROS. Also, these two genes worked as a hub by linking Modules B-2 and C as well. Furthermore, OXPHOS and ROS also communicated with JAK-STAT pathway through IL-6 and S100A6 connecting Modules A and C (Fig. [87]5). These findings suggest that the OXPHOS and ROS process may play a critical role in distinguishing the differential responses in cancer cell death by L-pampo. Figure 5. [88]Figure 5 [89]Open in a new tab An illustrative diagram represents the inferred signaling transduction for L-pampo-induced PC-3 cell death L-pampo engages Toll-like receptors (TLRs) of cancer cells and triggers two primary signaling pathways: PI3K-AKT and JAK-STAT. Activation of these pathways appears to increase oxidative phosphorylation (OXPHOS) and reactive oxygen species (ROS) production through hub genes, which leads to L-pampo-induced cancer cell death, especially in PC-3. The color of the box means the module ID (Module A: red, Module B: brown, and Module C: magenta). Discussion This study showed that L-pampo induces genes related to the ribosome, OXPHOS, ROS, JAK-STAT signaling and cytokine-cytokine receptor interaction in immune cells and cancer cells. These signaling pathways are essential for cell growth and proliferation, biochemical processes, and induction of immune responses. Among these pathways, it is important to narrow down potential molecular mechanisms leading cancer cells (i.e., PC-3) to cell death by L-pampo. Necroptosis in cancer cells can be a candidate pathway that can cause cellular toxicity and this pathway is related to TNF, RIPK1/3, and MLKL signaling cascades^[90]17. However, the member genes related to necroptosis showed low expression levels in PC-3 (data not shown) although they appear differentially expressed across the time points. This result suggests that necroptosis may not be the major cell death mechanism induced by L-pampo. On the other hand, OXPHOS and ROS signaling was located at the top in our pathway enrichment analysis (Fig. [91]2). The OXPHOS process is essential for cell death, proliferation, cellular redox balance, and metabolism^[92]28, and this pathway has been reported to be prominently activated in prostate cancer and neurodegenerative diseases^[93]29. It has been known that the OXPHOS produces more oxidative stress and induces apoptosis in cancer cells^[94]30. Given that, several anti-cancer drugs targeting OXPHOS have been developed to induce cancer cell death using this mechansim^[95]31,[96]32. In addition, as a downstream cell stress signaling of ROS can lead to cell death via NOD-like receptor pathway genes, which were identified as DEGs from our analysis^[97]33. NOD-like receptor pathway is also known to regulate proinflammation and cell death mechanism via apoptosis^[98]34. In addition to this mechanism, we also found that cell progression and survival seems to be suppressed by L-pampo in PC-3, as BCL2, an important progression marker and regulator in cancer, was significantly down-regulated at 3h and 6h in this cancer cell line. This suggests that the cell stress owing to OXPHOS with ROS may be augmented by the suppression in cancer cell progression. Then, we attempted to infer how OXPHOS and ROS is further activated by L-pampo. We found clues on the potential association between PI3K-AKT and OXPHOS and ROS pathways as it appears that PIK3CD and IL-12 bridge PI3K-AKT and OXPHOS/ROS. A recent study has reported that increased AKT activity leads to increased mitochondrial respiration and causes excessive accumulation of ROS, resulting in cell death in chronic lymphocytic leukemia^[99]35. These studies support the idea that increased PI3K-AKT signaling via TLRs can up-regulate OXPHOS and ROS and lead to cancer cell death. Additionally, dying cancer cells produce damage-associated molecular patterns (DAMPs) such as HMGB and calreticulin that stimulate immune cells, leading to additional cancer cell death in the tumor microenvironment^[100]36. As a whole, we suggested that OXPHOS and ROS pathways may play a critical role in cell death triggered by L-pampo. Although the signal conduction from TLRs to OXPHOS and ROS through PIK3-AKT is common to THP-1 (i.e., immune cells) and PC-3 (i.e., prostate cancer cells), the distinct expression patterns of OXPHOS and ROS genes may take special place in prostate cancer. It has been reported that OXPHOS and ROS pathways are considered to be important to prognosis and treatment in prostate cancer^[101]29,[102]37,[103]38. Thus, we promptly investigated whether the basal expression levels of OXPHOS and ROS genes in prostate cancer patients (TCGA-PRAD) could be associated with cancer prognosis. The overall survival rate of prostate cancer is very high with a five-year survival rate of 97.5%, yet approximately 30% of patients relapse, making biochemical recurrence (BCR) as a critical prognostic indicator^[104]39. We conducted a BCR-free survival analysis comparing the lowest quartile (Q1) and the highest quartile (Q4) of single sample gene set enrichment analysis (ssGSEA) scores of OXPHOS and ROS genes. As a result, the patients with the highest ssGSEA scores (Q4) have significantly lower recurrence rates (HR = 0.37, p-value = 0.04) (Figure S6). This data implies that highly expressed OXPHOS and ROS genes may be related with lower biochemical recurrence rates in prostate cancer. It would be worthwhile in future research to study the potential roles of OXPHOS and ROS in prostate cancer. Finally, further research is needed to reveal the detailed molecular underpinnings of how L-pampo enhances the efficacy of cancer immunotherapy, such as cancer vaccines and chemo/radio immunotherapy combination^[105]40. Conclusion Understanding key disease mechanisms or drug’s mode of action by analyzing molecular data such as RNA-seq expression profiles, may not be trivial because critical gene–gene interactions are often buried by other complex signaling flows in the cell. Moreover, interpreting time-course gene expression profiles is even more challenging with the increased number of comparisons between different time points. In this study, we devised a sophisticated network-based approach to compare three cell lines at different timepoints on the same reference (i.e., the template network). Using this approach, we discovered potentially important signaling pathways that may be lethal to TLR2/3 expressed cancer cells with L-pampo treatment (i.e., PC-3). From our analysis, the molecular stimulation due to OXPHOS and ROS pathways can lead cancer to cell death. The genuinely high expression levels of OXPHOS and ROS genes in PC-3 also supports the observation that PC-3 is more sensitive to L-pampo than SW620. In conclusion, we believe that our study provides valuable insights into the development of immune-based cancer therapies, which can produce potentially better outcomes for cancer patients. Moreover, our computational methods can generate more concise and testable hypotheses regarding disease or drug mechanisms. Therefore, we expect that these methods will also inspire computational biologists or AI experts who need to delve into notoriously complex molecular data. Materials and methods Cell culture and L-pampo treatment THP-1, PC-3, and SW620 cells were obtained from the American Type Culture Collection (ATCC, Rockville, MD, USA). THP-1 cells were grown in RPMI1640 medium (Gibco, Carlsbad, CA, USA) containing 10% fetal bovine serum (Gibco, Carlsbad, CA, USA), 1% antibiotics (100 mg/L streptomycin, 100 U/mL penicillin), and 0.4% 2-Mercaptoethanol (Sigma Aldrich, St. Louis, MO, USA). THP-1 cells were differentiated into macrophage-like cells in a medium with 1% phorbol 12-myristate 13-acetate (PMA) (Sigma Aldrich) for 72 h. PC-3 and SW620 cells were cultured in a DMEM medium (Gibco, Carlsbad, CA, USA) containing 10% fetal bovine serum and 1% antibiotics. All cells were cultured at 37°C in a 5% CO[2] atmosphere. L-pampo (Lot#230,524–01 batch) was produced according to the manufacturing process SOP of the CHA Vaccine Institute. Each cell line was plated at 1 × 106/well in a 6-well plate and incubated in a 5% CO[2], 37°C for 24 h. After replacing the fresh complete media and incubating for 1h, the diluted L-pampo was treated to each cell line at a concentration of 1x (40 µg/ml) and cultured for 3 and 6 h, respectively. Cell viability test (WST assay) Each cancer cell line was plated at 5 × 105/well in a 24-well plate and incubated in a 5% CO[2], 37°C for 24 h. After replacing the fresh complete media and incubating for 1h, the diluted L-pampo was treated to each cell line and cultured for 24 h. The culture medium was removed and washed twice using 1 × DPBS. Then, the wells were incubated with a CCK-8 working solution (5µl per well) (Dojindo Laboratories Co., Kumamoto, Japan) for 1h. Converted orange formazan dyes from WST optical densities were measured at 450 nm with a microplate reader. Total RNA isolation The cells were washed twice using 1 × DPBS (WelGENE, Seoul, South Korea) to isolate the total RNA. After removing all the DPBS, 1ml of TRIzol was added (Invitrogen, Waltham, MA, USA) to each well. The TRIzol/cell lysate was collected in a 1.5 mL EP tube after briefly scraping each well and incubated at room temperature for 5 min. Then, 250 µl chloroform was added into the EP tube and shaken vigorously 10 times. The EP tube was centrifuged at 14,000 rpm at 4°C for 5 min. The cleared supernatant was transferred to a new EP tube, and 500 µl of isopropanol was added, mixing it carefully. The EP tube was centrifuged at 14,000 rpm at 4°C for 20 min. The isopropanol was discarded, and the pellets were allowed to air dry. For RNA-seq analysis, the RNA pellet was distributed to DEPC-DW at an appropriate concentration. RNA library preparation and RNA sequencing (RNA-seq) RNA quality verification, library preparation, and sequencing were performed at Macrogen (Seoul, Korea). RNA samples showed over 9.7 RNA integrity number (RIN), and libraries were generated using the Illumina TruSeq Stranded mRNA kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Barcoded samples were pooled and sequenced on NovaSeq6000 (Illumina, San Diego, CA, USA) to produce 100 bp paired end reads. This experiment was performed on 3 cell lines, 3 time points, and 3 technical replicates, resulting in 27 samples. RNA-seq data analysis In the first step, the sequencing quality of raw reads was checked by FASTQ^[106]41. After trimming adapter sequences using TrimGalore (Cutadapt version 0.6.10)^[107]42, resulting reads were aligned to the human reference genome (GRCh38) by STAR (version 2.7.10b)^[108]43 and quantification of transcript abundances was performed using RSEM (version 1.2.28)^[109]44. Finally, differentially expressed genes (DEGs) analyses for each cell line and time point were conducted by limma (R package, version 3.54.0)^[110]45, then we identified the genes with p-value < 0.05 and the absolute of log-fold-change > 1 as DEGs. Functional enrichment analysis Over-representation analysis (ORA) using DEGs or gene modules from WGCNA from each condition was performed for KEGG pathways using clusterProfiler (R package, version 4.6.2)^[111]46. The significance score of pathway enrichment was defined as [MATH: -log10(adj.Pvalue) :MATH] , and the adjusted p-value was applied by Benjamini-Hochberg (BH) procedure, also known as the False Discovery Rate (FDR) procedure. Influence analysis by network propagation Network propagation is a graph-based analysis approach, which was proposed in PropaNet.^[112]10. The main process is to propagate information of a node to adjacent nodes via the edges at each iteration. This process is repeated for a predetermined number of iterations or until convergence. As a result of the network propagation, the value of a node influences both the values of distant nodes and those of close neighbors. For a given network [MATH: N :MATH] , we can think that [MATH: I0(n) :MATH] is a degree of influence of a seed node [MATH: n :MATH] at the beginning of iteration 0. At each iteration [MATH: i :MATH] , the influence of the node can be estimated by the sum of those at the adjacent nodes in the network at iteration [MATH: i-1 :MATH] . Then, the following equation can be defined. [MATH: Iin=∑m∈NIi< mo>-1m :MATH] 1 When the random walk with restart (RWR) is applied to Eq. [113]1, [MATH: Ii=αI< /mi>0+(1-α)Ii-1 :MATH] 2 where the restart rate [MATH: α :MATH] is considered as the weight of prior information by restarting at seed nodes again. After iteration, the resulting [MATH: I :MATH] represents a diffusion result from the initial seed nodes at each node. Weighted gene co-expression network analysis (WGCNA) WGCNA (R package, version 1.72–5) was used in not only constructing the template network upon gene co-expression but also identifying gene modules^[114]47. To set an appropriate power value parameter that satisfies the scale-free network, the ‘pickHardThreshold’ and ‘pickSoftThreshold’ functions of the WGCNA R package were used. A hard threshold on the Pearson correlation coefficient (PCC) between gene expression was applied while constructing a template network. Meanwhile, the basic procedures of WGCNA were followed in the module identification step. A soft threshold on a topological overlap matrix (TOM) dissimilarity was determined and then gene modules were detected using the ‘blockwiseModules’ function, which enables hierarchical clustering analysis in one step. Subnetwork detection using TOPAS algorithm The TOP-down Attachment of Seeds (TOPAS) algorithm was proposed in a recent study^[115]23. The algorithm assumes that if a node in a network is found on the shortest path between seed nodes, it is considered as a potential seed connector. Given the template network, it extracts the largest seed connected subgraph with the shortest path. In the next step, RWR is applied in the same way as Eq. [116]2 to calculate the degree of influence from the seed nodes in the subgraph. In the last step, it prunes subgraph by iteratively deleting the nodes with the lowest influence and maximizing the length of seed nodes in the pruned subgraph. Submodule detection using Walktrap algorithm Walktrap algorithm was implemented^[117]21 for detecting submodules from large-sized modules using the igraph R package (version 1.3.1)^[118]48. The algorithm assumes that if two nodes are in the same community, their random walks to the nodes in other communities will be similar. In the first step, every node has an individual community and in the next step, the community organization is updated as the two closest communities are merged based on their distance. When we suppose that there are k communities ( [MATH: C1,⋯,Ck :MATH] ) in the network, the distance between two communities [MATH: i :MATH] and [MATH: j :MATH] is calculated as follows: [MATH: γij=∑k=1n(Pikt< mo>-Pjkt)2d(k) :MATH] 3 where [MATH: Pijt :MATH] is the probability of going from [MATH: i :MATH] to [MATH: j :MATH] through a random walk of length [MATH: t :MATH] ^[119]49. Quantification of network differences Jensen-Shannon Divergence is a statistical distance between two probability distributions similar to Kullback–Leibler (KL) divergence, but modified to make the KL divergence symmetric and finitely constrained^[120]24. The nJSD (network Jenson Shannon Divergence) was defined as the sum of entropy values measured at each of the genes in a protein interaction network^[121]50. To quantify the distance between networks under different conditions, the nJSD values were measured using the absolute values of gene expression change (log2-fold-change) as a probability distribution. Given two networks A and B consisting of [MATH: n :MATH] genes, a gene node [MATH: i :MATH] has neighbor genes [MATH: Ji :MATH] and [MATH: ej :MATH] is the expression value of the gene [MATH: j :MATH] . Then, the probability distribution [MATH: Pi :MATH] and the KL divergence ( [MATH: DKLi :MATH] ) of gene [MATH: i :MATH] are defined as follows: [MATH: pij=ej ∑l∈Jiel, whenj∈Ji0,other wise :MATH] 4 [MATH: PiA=pi1,< /mo>⋯,pij,⋯,pin :MATH] 5 [MATH: DKLi[PA|PB=∑l=1nPi(A)log Pi(A)Pi(B) :MATH] 6 Finally, the following equation is defined for the nJSD between network [MATH: A :MATH] and [MATH: B :MATH] . [MATH: nJSD=1n∑i=1n12DKLi[PA|PB+12DKLi[PB|PA :MATH] 7 Single sample gene set enrichment analysis (ssGSEA) and survival analysis We downloaded the RNA-seq expression data (TCGAbiolink^[122]51) and clinical data (UCSC Xena^[123]52, [124]https://xenabrowser.net/datapages/) of prostate cancer patients (TCGA-PRAD, n = 426). Single sample gene set enrichment analysis (ssGSEA) was run using ‘GSVA’ (R package, version 1.53.3)^[125]53 with the TPM (Transcripts Per Million) value of the RNA-seq data was used as input to measure the relative enrichment levels of OXPHOS and ROS genes in the WGCNA modules for TCGA-PRAD tumor samples. The biochemical recurrence (BCR) was used as a prognostic indicator for prostate cancer. Survival analysis was performed on BCR using ‘survminer’ (R package, version 0.4.9, [126]https://rpkgs.datanovia.com/survminer/index.html). For this, TCGA-PRAD patients were equally divided into four groups (Q1-Q4) based on their ssGSEA scores and the lowest quartile (Q1) and the highest quartile (Q4) were compared using Cox proportional hazard model. Supplementary Information [127]Supplementary Figures.^ (5.1MB, pdf) Acknowledgements