Prognostic Pathways Guide Drug Indications in Pan-Cancers

Pathway-level analysis is a powerful approach enabling the interpretation of post-genomic data at a higher level than that of individual molecules. Molecular-targeted therapy focusing on cascade signaling pathways has become a new paradigm in anticancer therapy, instead of a single protein. However, the approaches to narrowing down the long list of biological pathways are limited. Here, we proposed a strategy for in silico Drug Prescription on biological pathways across pan-Cancers (CDP), by connecting drugs to candidate pathways. Applying on a list of 120 traditional Chinese medicines (TCM), we especially identified the “TCM–pathways–cancers” triplet and constructed it into a heterogeneous network across pan-cancers. Applying them into TCMs, the computational prescribing methods deepened the understanding of the efficacy of TCM at the molecular level. Further applying them into Western medicines, CDP could promote drug reposition avoiding time-consuming developments of new drugs.

Introduction

Advances in RNA sequencing technologies provide unprecedented capacity to comprehensively identify altered genes and pathways implicated in tumorigenic processes, raising the hope of extending targeted therapies (1–3). Monitoring the transcriptome of tumor samples mirrors voluminous altered pathways during cell programming and even heterogeneous cancer initiation and progression (4–6). Pathway-level transcriptomic analysis is a powerful approach to interpreting post-genomic data at a higher level than that of individual biomolecules (7, 8).

Tumor innate or acquired resistance leads to treatment failure (9). However, the development of new cancer therapeutics currently requires a long and protracted process of experimentation and testing (10). Recent research showed that the proportion of cancer patients who could be tractable by approved agents following clinical indications is merely 5.9% in a cohort of 4,068 tumors (11). Most of the single-agent drugs under development today target a specific protein (12). However, the formation and progression of a tumor usually occur through a range of defects in endothelial cells, immune cells, stromal cells, and cancerous cells (13). Defects in the whole signaling pathways allow cancer cells to alter their normal programmers of proliferation, transcription, growth, migration, differentiation, and death (14). Therefore, targeting multiple molecules, i.e., cascades of signaling pathways involved in both tumor and surrounding tissues, appears to be a new paradigm for anticancer therapy, rather than a single protein (15, 16).

The high failure rate of chemotherapy and the tread toward multitarget paradigms emphasize the importance of discovering novel effective therapies. Discovered pathways have been proposed to accelerate the development and approval of effective and targeting cancer therapies, such as the PI3K/AKT/mTOR pathway in bladder cancer (17). Apoptosis pathway (18), Src kinase pathway (19), Hippo signaling pathway (20), Wnt pathway (21), and DNA repair pathway (22) were all proposed as therapeutical targets in cancers. Besides in cancer stem cells, Ivy’s group reviewed that targeting Notch, Hedgehog, and Wnt pathway to control stem-cell replication, survival, and differentiation is a new strategy under development (23). As emerging anticancer drug targets, all these pathways could make contributions to the development of the clinical treatments for cancers.

Here, we developed a strategy for in silico Drug (including herbal medicines) Prescription on cellular pathways in Cancers (CDP). Our CDP connected drugs (agents or perturbations or small molecules) to pan-cancers on a pathway level to discover novel effective therapies. First, CDP established a pathway dysregulation profile and identified pathway modes of action. Second, CDP determined the magnitude of pathways across 20 cancers based on a generality score, and the tumor specificity of pathways based on prognostic evidence. Third, we linked pathways with 1,309 drugs and connected cancers to all potential beneficial agents. We proposed two applications of CDP in drug design. First is the reposition of well-known drugs for new indications, which can not only save time and cost in the process of research and development but also avoid the side effects. Second is the combination of traditional Chinese medicines (herbal components) to prescribe new potential formula targeting cellular pathways.

Results

Pan-Cancer Pathway Panels Reveal Tumor Connections

To portray the molecular pathway-level characterization of pan-cancer, we collected and integrated case–control studies of gene expression profiles and clinical outcome data from the Cancer Genome Atlas (TCGA) including 20 cancers. The paired pan-cancer cohorts contained 1,626 samples covering more than 10 distinct malignant histologies (Figure 1A, inlet). One tumor type corresponded to one case–control study, consisting of 6 to 224 normal and tumor samples (Figure 1B).

FIGURE 1

Figure 1 The brief of data source in this study. (A) The number of samples in the pan-cancer cohorts, one case–control study for each cancer type. (B) The tissues involved in pan-cancer, including livers, lungs, stomachs, thyroids, pancreas, breasts, kidneys, and rectal and colons. (C) The distribution pathway dysregulation reflecting by normalized enrichment scores (NES) in the form of boxplot. (D) Hierarchical clustering of various according to gene expression profiles of all detected genes, using the correlation and average linkage.

We constructed an activity profile of pathways based on those gene expression profiles of case–control studies across 20 cancers. The activity profile consists of normalized enrichment scores (NESs) for 178 pathways, representing the dysregulation (enrichment or depletion) score of each pathway in each tumor type. The cancer-related NES matrix was referred to as CNS.

The molecular characterization of pan-cancer at the pathway level revealed heterogeneity among different tumor tissues and consistency between tumor and supportive tissues. Hierarchical averaged-linkage clustering of 20 tumors showed 5 subcategories, with proximal to relevant tissues and organs or with similar physical functions: 1) READ, COAD, and COADREAD; 2) breast invasive carcinoma (BRCA), UCEC, LUAD, and LUSC; 3) cholangiocarcinoma (CHOL) and liver hepatocellular carcinoma (LIHC); 4) KIRP, KIPAN, and KIRC; and 5) esophageal carcinoma (ESCA) and stomach adenocarcinoma (STAD) (Figures 1C, D). PAAD, PRAD, and THCA seem to be quite different from other tumors across the deregulated features of KEGG pathways. Although these three cancers did not converge in one branch, evidence supported their underlying links (24). For most cancers, the corresponding distribution of NES indicated the overall level of pathway dysregulation within each subcategory (Figure 1C), consistent with the hierarchy shown in Figure 1D. For example, CHOL and LIHC are both hepatobiliary tumors. STAD and ESCA are both gastroenterological tumors. The abbreviations of all above cancers are shown in Table 1. Therefore, pathway panels in pan-cancers hinted the potential links underlying various tumors, suggesting subcategories among cancers.

TABLE 1

Table 1 The numbers of samples in each cancer type.

Pan-Cancer Pathway Panels Reveal Mode of Action (MoA) of Pathways

To investigate the pan-cancer pathway dysregulations across pan-cancers, the 178 KEGG pathways were partitioned into three clusters based on the CNS matrix (Figures 2A, B): loss of functions (LoF, Cluster 1), especially pathways harboring tumor-suppressive genes or downregulated genes in tumors, whose disruption facilitates tumorigenesis (Figure 2C, Figure S1); the LoF cluster which consists of 38 pathways including Hedgehog signaling, TGF-β, and WNT signaling; and switch of functions (SoF, Cluster 2) which defines the pathways whose favorable or adverse roles switched across different tumors (Figure 2D, Figure S2). A set of 62 cellular pathways were categorized into the SoF cluster, in which cellular signaling pathways accounted for a quarter. The well-known NOTCH, MAPK, mTOR, VEGF, JAK/STAT, and ERBB were grouped into the SoF cluster, indicating changes in their roles across various tumors. In addition, immune processes including B/T cell receptor signaling, chemokine signaling, and epithelial cell signaling in infection also play different roles in the tumor initiation and progression across various cancers (25). Gain of functions (GoF, Cluster 3) defines typically pathways consisting of oncogenic genes or unregulated genes in tumors, whose abnormal activation initiates tumorigenesis (Figure 2E, Figure S3). There were totally 72 cellular pathways categorized into the GoF cluster, including cell cycle, p53 signaling pathway, and basal cell carcinoma. Figures 2A, B shows the distribution of the dysregulations in each mode of action (MoA) in a boxplot and heatmap.

FIGURE 2

Figure 2 Profiles of pathway dysregulations on expression across pan-cancers. (A) The distribution pathway dysregulation reflecting by normalized enrichment scores (NES) in the form of boxplot. (B) The distribution pathway dysregulation reflecting by normalized enrichment scores (NES) in the form of heatmap. Three modes of action (MoA) of pathways: loss of functions (LoF) (C), switch of functions (SoF) (D), and gain of functions (GoF) (E). Hierarchical clustering activity profiles (NESs) of 178 biological pathways based on GSEA, using the Euclidean distance metric and average linkage.

To validate the oncogenic and tumor-suppressive characteristics of GoF and LoF pathways, we detected the differential expressed genes (DEG) in comparison of tumor and normal samples for each cancer type. Considering the robustness and the number of DEGs, we determined the optimal threshold (| log₂(fold change) | > 2, FDR <0.01) for the differentially expressed analysis. Figures S4, S5 show the variations of the number of upregulated and downregulated genes along with the different thresholds, respectively. For most cancer types, the significant enrichments of upregulated genes (tumor vs. normal) in the pathways of the GoF cluster were observed (Figure 3A, hypergeometric test, p < 0.05). The pathways whose members did not overlap with the DEG (upregulated or downregulated) are not shown in Figures 3A, B. Similarly, downregulated (tumor vs. normal) genes were significantly enriched in the pathways of the LoF cluster (Figure 3B, hypergeometric test, p < 0.05). The significant enrichments validated the aforementioned MoA hypothesis of GoF and LoF across pan-cancer. Identifying the MoA of cellular pathways was crucial to exploring the opportunities to therapeutically target pathways, because in principle, only pathways with activity (LoF, GoF, or SoF) could be directly targeted with molecules able to inhibit their dysfunctions.

FIGURE 3

Figure 3 Pathway dysregulations on functions across pan-cancers. (A) The enrichment on the set of somatic mutations genes (SMGs) of differential expression genes in each cancer type. (B) The enrichment on the set of drug targets (DTs) of differential expression genes in each cancer type. Hypergeometric tests were performed to obtain statistical significance, i.e., p value. The blank cells in the heatmap means that the p value was 1. The darker the green color, the smaller the significant p value of the hypergeometric test. (C) The overlap of SMG-enriched or DT-enriched pathways on upregulated or downregulated pathways across pan-cancers, which indicated the potential of a pathway to be a target based on its mutation burden or target enrichment. (D) The distribution of generality scores across three clusters of pathways according to three MoAs.

Characteristics of High Mutation and Target Enrichment Concentrate Seven Important Pathways Across Pan-Cancer

To detect the important KEGG pathways in each cancer, we proposed the following two questions: (1) whether significantly mutated genes (SMG) were enriched in the GoF pathways. The assumption proposed by Gonzalez implied that mutations in certain genes provided a selective advantage to tumor cells that followed by growing and proliferating faster (26). Tracing the signals left by the selection benefited to identify the gene driving tumorigenesis (27). Therefore, the enrichments of mutational drivers on GoF pathways suggested the pathway driving tumorigenesis. (2) Considering the pathway-targeted therapies for cancer, we asked the second question, whether drug targets (DT) were enriched in the GoF pathways. The enrichments of common DTs on GoF pathways indicated their potential of benefiting in cancer treatments.

The SMG list consisted of 585 driver genes collected from two pan-cancer studies (28). The list of 1,869 DTs was curated from DrugBank which currently consisted of ~10,000 drug entries (29, 30). The hypergeometric test was performed for three MoAs of pathways, and the statistically significant p values are as shown in Figures 3A, B. Using the threshold of p values above 0.05 (hypergeometric test), SMG were enriched in 47 pathways and DT were enriched in 15 pathways. GoF pathways contributed to 11 and 7, respectively. Intersecting the upregulated and downregulated significant pathways with DT-enriched pathways showed two and eight pathways, respectively. Intersecting the upregulated and downregulated significant pathways with SMG-enriched pathways showed six and five pathways, respectively (Figure 3C, left panel). There were 7 pathways that were enriched for both mutational and drug-targeted evidence. The enriched pathways are involved in five cancer-related diseases, a metabolic disease, and two signaling pathways (Figure 3C, right panel). The only one pathway in the GoF cluster which was both DT-enriched and SMG-enriched was “Type_II_Diabetes_Mellitus,” indicating that the metabolic syndrome increased cancer risk which is consistent with the perspectives in the previous studies (31, 32).

Quantifying the Magnitude of the Pathways Across 20 Cancers and Within Specific Cancer

Twenty-Two Generalized, Important Pathways Across 20 Cancer Types

To quantify the magnitude of biological pathways across 20 cancer types, we defined a generality score which assigned the functional inclination of enrichment or depletion to pathways. The “generality score” quantifies the consistency of deregulated pathways across different cancer types and is a measure of the extent to which the levels of the pathway were perturbed by various cancer types. The higher the absolute value of the generality score, the larger the magnitude of a gene set or a pathway is. The sign of generality means the enrichments or depletions which indicated the positive or negative functional inclination in most cancers.

The generality scores for each pathway in the three clusters of three MoAs are shown in Figure 3D. The top 11 generally enriched and depleted pathways across various cancers are shown in Table 2. The most common enriched pathways were implicated in the fundamental biological processes maintaining the cell survival, cell cycle, DNA replication, and homologous recombination. In contrast, majority of common depleted pathways were the metabolic processes including tryptophan, propanoate, fatty acid, and tyrosine. The distribution of all generality scores across 178 KEGG pathways is shown in Figure S6.

TABLE 2

Table 2 The top 11 generally enriched and depleted pathways across various cancers.

To verify whether the generality could represent the features of the common pathways across 20 cancers, we counted the numbers of cancer types in which the NES was above 1 for the above 22 generalized pathways (Table 2). All the generalized, depleted pathways were observed to have NES of less than -1 in more than half of 20 cancer types. 9/11 of the depleted pathways (butanoate metabolism, valine leucine, and isoleucine degradation) were observed to have NESs of less than -1 in the third quartile (top 75%, 15 cancer types) of NES across 20 cancer types. As for the generalized, enriched pathways, they were all observed to have NESs exceeding 1 in more than half of 20 cancer types. 1/11 of the enriched pathways (basal transcription factors) was observed to have NESs of more than -1 in the third quartile (top 75%, 15 cancer types) of NESs across 20 cancer types. Figure S7 shows the NESs of each pathway across 20 cancers.

According to the generality profiles in Figure 3D, we could reassign 4 pathways to the new MoA cluster more exactly (Figure 3D). In hierarchical clustering, we labeled O-glycan biosynthesis as the LoF cluster. However, the generality of 7.2 showed the proximity of cluster GoF. The “glyoxylate and dicarboxylate metabolism,” “alanine aspartate and glutamate metabolism,” and “selenoamino acid metabolism” were reassigned to the LoF cluster in spite of being possessed by the GoF cluster in the unsupervised clustering in the former.

Therefore, the generality we proposed was a novel measure to reflect the commonness of a gene set or a pathway among multiple cancer types. Generality benefits the development of pathway-targeted therapy in drug combination and prescription of herbal medicine.

Tumor-Specific, Important Pathways for Each Cancer

The pathway-targeted therapy required detections of not only the generalized or shared pathways but also tumor-specific pathways activated in specific cancer (33). Intratumor heterogeneity explained differences both in molecular phenotype and clinical outcomes (34). The follow-up data from TCGA showed KICH with the highest 5-year survival rate (up to 73%) and STAD with the lowest 5-year survival rate (merely 13%). The survival and matched transcriptomic data (RNA-seq) from TCGA enabled us to perform the pathway-level survival analysis to identify prognostic pathways for each of the 20 cancer types. The pathway-level survival analysis revealed the prognostic implication of pathways which reflect the molecular phenotype in specific tumor.

To quantify the activation of each cellular pathway in specific patients, we defined a pathway activation (PA) score by single sample extension of GSEA (ssGSEA). ssGSEA measures the overexpression of a pathway in each patient by comparing the rank of the genes in the pathway with all other genes. The patient cohort of each cancer was stratified into two groups according to a given threshold of the PA score. To obtain the optimal cutoff, we iteratively applied a series of thresholds to a log-rank test, examining the most significant differences in survival outcomes between two groups. The PA score with the lowest log-rank p value was considered as the optimal cutoff. In total, 3,560 Kaplan–Meier (KM) plots corresponded to all 178 KEGG pathways across 20 cancer types. Figure 4 shows the cellular pathways which were generalized enrichments or depletions in breast cancer according to the generality score.

FIGURE 4

Figure 4 The examples of prognostic pathways identified in survival analysis in BRCA. Adverse prognostic pathways: basal transcription factors, cell cycle; favorable prognostic pathways: fatty acid metabolism, adipocytokine signaling pathway. The adverse and favorable pathways are defined as follows: 1) favorable prognostic pathways, for which a higher PA score of a pathway was correlated with a good/longer survival outcome, and 2) adverse prognostic pathways, for which a higher PA score of a pathway was correlated with a poor survival outcome.

The cellular pathway of a cancer was defined as a prognostic pathway with PA scores higher or lower than the optimal cutoff in an individual patient which yields a significant (p < 0.001, log-rank test) difference in overall survival. Totally, 319 significant prognostic pathways were defined in 20 cancers. The number of unique prognostic pathways was 152. We detected two types of prognostic marker pathways in terms of clinical outcomes: (1) favorable prognostic pathways, for which a higher PA score of a pathway was correlated with a good/longer survival outcome, and (2) adverse prognostic pathways, for which a higher PA score of a pathway was correlated with a poor survival outcome. The proportion of favorable and adverse prognostic pathways varied in different cancer types.

By uncovering tumor-specific and generalized prognostic pathways in pan-cancers, exploring the magnitude of the cellular pathways helps to accelerate the discovery of targetable pathways in cancers and facilitate generations of novel hypotheses.

A Strategy for Computational A Strategy for Computational Drug Prescription in Cancers (CDP)

To connect therapeutic agents with pathways, we introduced the 1,309 perturbations (approved drugs, small-molecule compounds, and herbal medicines) obtained from the connectivity map (cMap build 02, https://clue.io/cmap) as a drug library. In each cancer type, we performed “cMap analysis” (35) to obtain the normalized enrichment scores (NESs) of all 1,309 perturbations, which suggests the possible treatment efficiency scores based on transcriptomic level. Across all the 20 cancer types, a matrix consisting of 1,309 × 20 was produced and defined as perturbation-related NES (PNS).

According to the same cancer type, we combined CNS that we have defined before with PNS, to connect all possible pairwise perturbations and pathways. We constructed a landscape of heterogeneous panels including perturbations and pathways associated with cancer types (Figure 5A). The heterogeneous panels might provide evidence of possible effects of perturbations and pathways on specific cancer type. The basic hypothesis is that there are oncogenic and prognostic pathways in which the gene expression level is significantly upregulated in cancers, and the potentially efficacious drugs are those perturbations which possess the capability of reversing the upregulated pathway signatures. This theory is similar to that of cMap, in which using a rank-based pattern-matching strategy is based on a small list of genes as a query signature (35, 36).

FIGURE 5

Figure 5 The connections of cancers–drugs–pathways. (A) The heterogeneous panels including perturbations (bottom panel, PNS matrix) and pathways (top panel, CNS matrix) associated with cancer types. A matrix consists of 1,309 drugs, and 20 cancers (1,309 × 20), representing the potential relationship of treatments, were produced and defined as perturbation-related NES (PNS). A matrix consists of 178 pathways, and 20 cancers (178 × 20), representing the dysregulation (enrichment or depletion) of each pathway in each tumor type, were defined as the cancer-related NES matrix (CNS). (B) The pan-cancer reversing network performed by Cytoscape, including 63 drug–pathway–cancer triplets.

To achieve the goal of prescription, i.e., screening the perturbations that are potentially useful for some cancer type on some pathway, we built a prescription score (PS) to measure the reversal potency of corresponding perturbations for a significant, prognostic pathway in some cancer type using the CNS and PNS matrix. The calculation was shown as follows: PS = absolute|CNS - PNS|. The construction of CNS and PNS and the calculation of PS consist of the strategy for computational drug prescription in cancers (CDP).

To extend the CDP strategy to application, we focused on traditional Chinese medicines (TCMs or herbs, herbal medicines). The first application of the CDP pipeline is on drug prescription for herbal medicines in traditional Chinese medicine. There were totally 120 herbal medicines in the cMap drug library. We identified prioritized herb–pathway–tumor triplets across pan-cancers and constructed a reversing network (Table S2, Figure 5B). The threshold to obtain the prioritized herb–pathway–tumor triplets consists of (1) a pathway that is significant, prognostic in survival analysis in at least one cancer type and (2) the absolute of the CNS score above 1.5 which means a pathway is enough upregulated or downregulated in the cancer. The network which is composed of prioritized cancer–drug–pathway triplets is illustrated in Figure 5B. The pan-cancer reversing network consisted of 63 drug–pathway–cancer triplets, involving 30 herbs and 42 significant, prognostic pathways across 13 cancer types.

Computational drug prescription in cancers (CDP) is a system-based approach to quantify the reversed effects of perturbations on tumors and predict the drugs showing high potency to reverse pathway activity in pan cancer.

Discussion

In this study, we firstly constructed a profile of molecular dysregulation at the pathway level. Then we portrayed the mode of action (MoA) influenced by compounds and herbal components on pathways and identified oncogenic, prognostic pathways in which the gene expression level is significantly unregulated in cancers and adverse to good prognosis of patients. The MoA comprehensively represented characteristics of distinct pathways, implying the capacity of becoming drug targets. LoF and GoF are commonly changed in some cancers, and drugs targeting these pathways can be shared among these cancers, while drugs targeting SoFs should be carefully investigated before further researches, since one pathway in the SoF cluster represents different profiles based on gene expressions. By reversing the expression of the upregulated, prognostic pathways, we developed a strategy for in silico drug prescription on cellular pathways in cancers (CDP) to achieve the goal of drug prescription, such as reposition of the well-known drugs for new indications and discovery of herbal components indicating for tumors. It is noted that the predicted results could only provide an alternative unless an indication was validated in vitro and in vivo.

Applying CDP to TCM, we constructed the pan-cancer reversing network which consists of herb–pathway–tumor triplets. To give evidence of the application, we searched well-described relationships of herbs–cancers in the published literatures. With the corresponding pathways, the molecular mechanism of the herbs could be further explored.

Taking BRCA as example, the significantly prognostic pathway signature consisted of 13 cellular pathways, of which the most upregulated, dysregulated pathway was “Basal_Transcription_ Factors.” According to the PS value calculated by the CDP strategy, we focused on the top 3 herb medicines (lycorine, podophyllotoxin, cephaeline) that potentially reversed the pathway of “Basal_Transcription_Factors” (Table 3). The published evidence indicating the treatment on BRCA is shown as follows.

TABLE 3

Table 3 The top PS value calculated by the CDP strategy.

Lycorine is a natural alkaloid extracted from the Amaryllidaceae genus. Ying et al. showed the suppression on cell growth, migration, and invasion both in multiple cell lines like MCF-7, T47D, MDA-MB-231, and 4T1 and in xenograft models (37). Besides the in vitro and in vivo studies, they also proposed that the mechanism of inhibition was that the blockage of Src/FAK (focal adhesion kinase), which is localized in nuclei, could play roles in regulating transcription factors (38). Similarly, Wang et al. also showed the findings about the inhibition of tumor growth and metastasis in breast cancer both in vitro and in vivo and investigated that suppressing STAT3 phosphorylation via upregulating SHP-1 expression was the cause of antitumor effects (39). Both the previous efforts provided the evidence that lycorine affected the regulations of transcriptional factors, which was represented as Basal_Transcription_Factors pathway in our CDP predictions.

Podophyllotoxin is an active component purified from Podophyllum peltatum (40). In terms of biological effects, podophyllotoxin is known to have immunosuppressive activity and antiviral effects against herpes, measles, influenza, and venereal warts (41). An in vitro study showed that podophyllotoxin presented the growth inhibitory activity for breast cancer cell lines and the inhibition extent up to 50% (42). Podophyllotoxin derivatives were investigated to inhibit the growth of both the breast cancer cell line MCF-7 and small-cell cancer cell line H69AR (43). The molecular mechanisms of podophyllotoxin in breast cancer were not investigated in current studies. However, in lung cancer, the anticancer effects of synthetic podophyllotoxin derivatives were verified in vitro and achieved by triggering mitotic arrest and ER stress (44).

Cephaeline was an analog of the natural product emetine, which was a kind of terpenoid tetrahydroisoquinoline alkaloids (45). Emetine is a protein synthesis/translation inhibitor that was evaluated in Phase II clinical studies as a potential chemotherapeutic agent for the treatment of solid tumors over 30 years ago (46). We did not find the literatures which reported the roles of cephaeline in BRCA, but Muhammad et al. indicated that cephaeline was cytotoxic to SK-MEL, KB, BT-549, and SK-OV-3 human cancer cells (45). Geir et al. proposed that cephaeline could be an inhibitor of tumors with necrosis, which was inferred from gene expression in aggressive endometrial cancers (47). In breast tumors, necrosis has been related to high-grade disease, increased tumor size, estrogen receptor negative status, high microvessel density, and macrophage infiltration (48). This evidence supports that our CDP strategy facilitates the drug prescription, especially for herbs, whose molecular mechanisms remain largely unknown.

There is a well-known computational method called cMap to connect diseases and small molecules. Using a list of differentially expressed genes (both upregulated and downregulated genes) in comparison of cases versus controls, cMap applies a computational method based on the Kolmogorov–Smirnov test to associated small molecules with the gene list, which reflected disease-related phenotypes. Our CDP has an advantage over cMap by taking biological pathways, especially prognostic pathways, into consideration, which approaches more to the nature of diseases.

Conclusion

In summary, we developed a score system for in silico drug prescription on cellular pathways in cancers (CDP) by constructing a profile of molecular dysregulation at the pathway level. Then, we portrayed the three modes of action (MoA), which were influenced by compounds and herbal components on pathways including LoF, GoF, and SoF. Application of this system to herbal components derived from TCM or herbs revealed potential indications targeting cellular pathways in pan-cancers, deepening our understanding of the molecular mechanisms of herbs. In addition, we could also repurpose the well-known (Western) drugs into new indications, which not only saves time and cost in the development process but also avoids side effects. The computational pipeline complements to the traditional target-based approaches at the pathway level, connecting cancers with potential effective drugs.

Materials and Methods

Data Source

Transcriptomic Data

Paired tissue samples (tumor and normal) across 1,626 individual cohorts were obtained from TCGA using the R package of “TCGA2STAT.” The pan-cancer-normalized gene-level RNA-Seq data for the TCGA cohorts were arranged, as shown in Table 1. With the transcriptomic data, the clinical data for all the TCGA cohorts were also obtained for survival analysis.

Agent-Perturbed, Transcriptomic Data

The CMap database (build 02): four types of human cancer cell lines (MCF7, PC3, HL60, SKMEL5) were treated with 1,309 distinct perturbations including US Food and Drug Administration (FDA)-approved drugs and uncharacterized bioactive compounds (called perturbagens by the authors of cMap, here simply referred to as drugs), for a total of 6,097 treatments. Gene (total RNA) expressions from the 6,097 “instances” (an instance is a cell line treated with a drug at a dosage, and its non-treated control) were recorded in two batches of microarrays: 671 HG-U133A (Affymetrix) chips (on 407 drugs) and 5,426 HT-HG-U133A chips (for a total of 6,097 chips on 1,309 drugs). Raw data were downloaded from the CMap website (http://www.broadinstitute.org/cmap/).

Cellular Pathways

Cellular pathways are a collection of manually drawn KEGG pathway maps representing experimental knowledge on metabolism and various other functions of the cell and the organism. Each pathway map contains a network of molecular interactions and reactions and is designed to link genes in the genome to gene products (mostly proteins) in the pathway.

Construction of Activity Profiles of Pathways

We employed Gene Set Enrichment Analysis (GSEA) to connect cancers to biological pathways from Kyoto Encyclopedia of Genes and Genomes (KEGG). We took the 186 KEGG pathways as the representatives of the general physiological processes. The pathway list was narrowed down into 178 pathways since those members of 8 pathways lacked gene expression levels in the paired pan-cancer cohorts. The NES (normalized enrichment scores) for 178 pathways were presented in a 178 × 20 table (Table S1), representing the dysregulation (enrichment or depletion) score of each pathway in each tumor type. This cancer-related NES matrix was referred to as CNS. On the basis of 1,000 trials, the empirical p values for the observed enrichment of 178 KEGG biological pathways were obtained. The false discovery rate (FDR) correction was used to deal with the accumulation of type I errors (false positives) in multiple comparisons.

Calculation of Generality Score

Generality is a measure that quantifies the deregulated commonness of a gene set or a pathway among multiple cancer types. Let n be the total number of cancer types in a study (20) and m the number of gene sets or pathways (178). For each pathway i in a collection of pathways, both a and b statistics are respectively computed corresponding to enriched pathways and depleted pathways, giving a_i and b_i. Actually, the generality score could be calculated for any gene sets or pathways. For simplicity and easy understanding, we only referred to pathways.

All cancer types in an ascending order were sorted by the normalized enrichment score of the gene set i, giving r_ik the rank of gene set i in the k^th cancer type, where k = 1, 2,…, n. A rank matrix (m rows × n columns) of the position (1… n) of each pathway in certain cancer was constructed. The ranks r_ik of n cancer types for m pathways in the current analysis are contained in the Supplementary File.

The a and b statistics were weighed using the proportion of significant enrichments and depletions in n cancer types of pathways i. Significant enrichments of pathway i require that the NES is positive and significance p is less than 0.05 in certain cancer type k. Significant depletions of pathway i require that the NES is negative and significance p is less than 0.05 in certain cancer type k. The number of cancer types for pathway i was counted according to the significant enrichments and deletions, giving w₁ and w₂.

The following two statistics corresponding to the enrichments and depletions of pathway i across n cancer types are computed:

$$a_i=\frac{w_1}{w_1+w_2}{\displaystyle \sum _{k=1}^n\left(\frac{r_{ik}}{n}\times 1\right)}$$

$$bi=\frac{w_2}{w_1+w_2}{\displaystyle \sum _{k=1}^N\left(\left(1-\frac{(r_{ik}-1)}{n}\right)\times (-1)\right)}$$

$$generalityscore{S}_i=\{\begin{array}{c}{a}_i(|a_i|>|b_i|)\\ b_i(|a_i|<|b_i|)\\ 0(|a_i|=|b_i|)\end{array}$$

Set S_i = a_i, if |a_i| > |b_i|, set S_i = b_i if |a_i| < |b_i|. The generality score is set to zero where a_i and b_i have the same value, which means that pathway i is not a generally deregulated pathway across these n cancer types.

Survival Analysis

All the patients along with the survival information were included in the Kaplan–Meier (KM) survival analysis. The patients who lacked survival information were excluded from survival analysis. The KM curves show the differences in overall survival for patients stratified by an ssGSEA score (49) of a certain pathway for each cancer type (see details in Supplementary Files). The KM curve separation was assessed by the log-rank test. The pathways which achieved p < 0.001 were considered as prognostic pathways. For each cancer and each pathway, the patient cohort was stratified into two groups with the higher and lower enrichments (ssGSEA), and an optimal split was selected by tracing each ssGSEA score to find the smallest p of the log-rank test. If the activations of a pathway were measured for each patient, his or her survival time and status (alive or dead) will be used to distinguish the prognostic implications.

Clustering Analysis

Pathway classification was identified using hierarchical clustering 178 pathways with the normalized enrichment score (NES) across all tumor types. The hierarchical clustering used the Spearman correlation metric and average linkage. The modes of action of pathways were identified using agglomerative hierarchical clustering of pathways with an SD > 1.0 on the log2 scale across all pathways. Consensus was commutated using 500 samplings of the pathways, 80% of the pathways per sample, agglomerative hierarchical clustering, Euclidean distance metric, and average linkage.

Differentially Expressed Analysis

To test for the differential expression across two groups (tumor and normal), we used the R package samr on gene expression profiles (FPKM value) from TCGA RNA-seq datasets. The p values were adjusted for multiple testing based on the false discovery rate (FDR). For comparison of two sample groups, we applied SAM (significance analysis of microarrays) which consists of a modification of the standard t-test by the inclusion of an extra variable s₀ called variance).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author Contributions

ZZ, YR and YL designed the work. KeZ, QX, RR,YP collected the data, FM, KZ, KenZ, CY, and YZ analyzed the data, YS, FM and ZZ interpreted results and wrote the manuscript, ZZ, YL, and HG reviewed the manuscript. FM participated in the preparation of the figures and tables. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the National Key R&D Program of China(No. 2017YFF0106006) and the National Natural Science Foundation of China (No. 82002647).

Conflict of Interest

Authors FM, YP, KZ, QX, RR, YZ, YS, YL, HG, and YR were employed by CapitalBio Technology.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank The Cancer Genome Atlas (TCGA) project and their contributors for these valuable public data sets.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.849552/full#supplementary-material

References

1. Reinartz S, Finkernagel F, Adhikary T, Rohnalter V, Schumann T, Schober Y, et al. A Transcriptome-Based Global Map of Signaling Pathways in the Ovarian Cancer Microenvironment Associated With Clinical Outcome. Genome Biol (2016) 17(1):108. doi: 10.1186/s13059-016-0956-6

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Zhao Z, Zhang KN, Wang Q, Li G, Zeng F, Zhang Y, et al. Chinese Glioma Genome Atlas (CGGA): A Comprehensive Resource With Functional Genomic Data From Chinese Glioma Patients. Genomics Proteomics Bioinf (2021) 19(1):1–12. doi: 10.1016/j.gpb.2020.10.005

CrossRef Full Text | Google Scholar

3. Zhao Z, Meng F, Wang W, Wang Z, Zhang C, Jiang T. Comprehensive RNA-Seq Transcriptomic Profiling in the Malignant Progression of Gliomas. Sci Data (2017) 4:170024. doi: 10.1038/sdata.2017.24

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Lim E, Wu D, Pal B, Bouras T, Asselin-Labat ML, Vaillant F, et al. Transcriptome Analyses of Mouse and Human Mammary Cell Subpopulations Reveal Multiple Conserved Genes and Pathways. Breast Cancer Res (2010) 12(2):R21. doi: 10.1186/bcr2560

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Yang L, Wang WH, Qiu WL, Guo Z, Bi E, Xu CR. A Single-Cell Transcriptomic Analysis Reveals Precise Pathways and Regulatory Mechanisms Underlying Hepatoblast Differentiation. Hepatology (2017) 66(5):1387–401. doi: 10.1002/hep.29353

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Li S, Li L, Meng X, Sun P, Liu Y, Song Y, et al. DREAM: A Database of Experimentally Supported Protein-Coding RNAs and Drug Associations in Human Cancer. Mol Cancer (2021) 20(1):148. doi: 10.1186/s12943-021-01436-1

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Draghici S, Khatri P, Tarca AL, Amin K, Done A, Voichita C, et al. A Systems Biology Approach for Pathway Level Analysis. Genome Res (2007) 17(10):1537–45. doi: 10.1101/gr.6202607

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Khatri P, Sirota M, Butte AJ. Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges. PloS Comput Biol (2012) 8(2):e1002375. doi: 10.1371/journal.pcbi.1002375

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sharma P, Hu-Lieskovan S, Wargo JA, Ribas A. Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell (2017) 168(4):707–23. doi: 10.1016/j.cell.2017.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Schein PS, Scheffler B. Barriers to Efficient Development of Cancer Therapeutics. Clin Cancer Res (2006) 12(11 Pt 1):3243–8. doi: 10.1158/1078-0432.CCR-06-0329

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Rubio-Perez C, Tamborero D, Schroeder MP, Antolin AA, Deu-Pons J, Perez-Llamas C, et al. In Silico Prescription of Anticancer Drugs to Cohorts of 28 Tumor Types Reveals Targeting Opportunities. Cancer Cell (2015) 27(3):382–96. doi: 10.1016/j.ccell.2015.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Warner KD, Hajdin CE, Weeks KM. Principles for Targeting RNA With Drug-Like Small Molecules. Nat Rev Drug Discov (2018) 17(8):547–58. doi: 10.1038/nrd.2018.93

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Swartz MA, Iida N, Roberts EW, Sangaletti S, Wong MH, Yull FE, et al. Tumor Microenvironment Complexity: Emerging Roles in Cancer Therapy. Cancer Res (2012) 72(10):2473–80. doi: 10.1158/0008-5472.CAN-12-0122

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Martin GS. Cell Signaling and Cancer. Cancer Cell (2003) 4(3):167–74. doi: 10.1016/S1535-6108(03)00216-2

PubMed Abstract | CrossRef Full Text | Google Scholar

15. McClellan M, Benner J, Schilsky R, Epstein D, Woosley R, Friend S, et al. An Accelerated Pathway for Targeted Cancer Therapies. Nat Rev Drug Discov (2011) 10(2):79–80. doi: 10.1038/nrd3360

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Faivre S, Djelloul S, Raymond E. New Paradigms in Anticancer Therapy: Targeting Multiple Signaling Pathways With Kinase Inhibitors. Semin Oncol (2006) 33(4):407–20. doi: 10.1053/j.seminoncol.2006.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sathe A, Nawroth R. Targeting the PI3K/AKT/mTOR Pathway in Bladder Cancer. Methods Mol Biol (2018) 1655:335–50. doi: 10.1007/978-1-4939-7234-0_23

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bai L, Wang S. Targeting Apoptosis Pathways for New Cancer Therapeutics. Annu Rev Med (2014) 65:139–55. doi: 10.1146/annurev-med-010713-141310

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Shukla D, Meng Y, Roux B, Pande VS. Activation Pathway of Src Kinase Reveals Intermediate States as Targets for Drug Design. Nat Commun (2014) 5:3397. doi: 10.1038/ncomms4397

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tao Y, Cai F, Shan L, Jiang H, Ma L, Yu Y. The Hippo Signaling Pathway: An Emerging Anti-Cancer Drug Target. Discov Med (2017) 24(130):7–18.

PubMed Abstract | Google Scholar

21. Anastas JN, Moon RT. WNT Signalling Pathways as Therapeutic Targets in Cancer. Nat Rev Cancer (2013) 13(1):11–26. doi: 10.1038/nrc3419

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Helleday T, Petermann E, Lundin C, Hodgson B, Sharma RA. DNA Repair Pathways as Targets for Cancer Therapy. Nat Rev Cancer (2008) 8(3):193–204. doi: 10.1038/nrc2342

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Takebe N, Miele L, Harris PJ, Jeong W, Bando H, Kahn M, et al. Targeting Notch, Hedgehog, and Wnt Pathways in Cancer Stem Cells: Clinical Update. Nat Rev Clin Oncol (2015) 12(8):445–64. doi: 10.1038/nrclinonc.2015.61

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Benedetti M, Zona A, Beccaloni E, Carere M, Comba P. Incidence of Breast, Prostate, Testicular, and Thyroid Cancer in Italian Contaminated Sites With Presence of Substances With Endocrine Disrupting Properties. Int J Environ Res Public Health (2017) 14(4):355–66. doi: 10.3390/ijerph14040355

CrossRef Full Text | Google Scholar

25. Niwa T, Tsukamoto T, Toyoda T, Mori A, Tanaka H, Maekita T, et al. Inflammatory Processes Triggered by Helicobacter Pylori Infection Cause Aberrant DNA Methylation in Gastric Epithelial Cells. Cancer Res (2010) 70(4):1430–40. doi: 10.1158/0008-5472.CAN-09-2755

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Gonzalez-Perez A, Lopez-Bigas N. Functional Impact Bias Reveals Cancer Drivers. Nucleic Acids Res (2012) 40(21):e169. doi: 10.1093/nar/gks743

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Gupta SC, Sung B, Prasad S, Webb LJ, Aggarwal BB. Cancer Drug Discovery by Repurposing: Teaching New Tricks to Old Dogs. Trends Pharmacol Sci (2013) 34(9):508–17. doi: 10.1016/j.tips.2013.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kandoth C, McLellan MD, Vandin F, Ye K, Niu B, Lu C, et al. Mutational Landscape and Significance Across 12 Major Cancer Types. Nature (2013) 502(7471):333–9. doi: 10.1038/nature12634

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Law V, Knox C, Djoumbou Y, Jewison T, Guo AC, Liu Y, et al. DrugBank 4.0: Shedding New Light on Drug Metabolism. Nucleic Acids Res (2014) 42(Database issue):D1091–7. doi: 10.1093/nar/gkt1068

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wishart DS, Knox C, Guo AC, Shrivastava S, Hassanali M, Stothard P, et al. DrugBank: A Comprehensive Resource for in Silico Drug Discovery and Exploration. Nucleic Acids Res (2006) 34(Database issue):D668–72. doi: 10.1093/nar/gkj067

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Azrad M, Turgeon C, Demark-Wahnefried W. Current Evidence Linking Polyunsaturated Fatty Acids With Cancer Risk and Progression. Front Oncol (2013) 3:224. doi: 10.3389/fonc.2013.00224

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Esposito K, Chiodini P, Colao A, Lenzi A, Giugliano D. Metabolic Syndrome and Risk of Cancer: A Systematic Review and Meta-Analysis. Diabetes Care (2012) 35(11):2402–11. doi: 10.2337/dc12-0336

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yauch RL, Settleman J. Recent Advances in Pathway-Targeted Cancer Drug Therapies Emerging From Cancer Genome Analysis. Curr Opin Genet Dev (2012) 22(1):45–9. doi: 10.1016/j.gde.2012.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Swanton C. Intratumor Heterogeneity: Evolution Through Space and Time. Cancer Res (2012) 72(19):4875–82. doi: 10.1158/0008-5472.CAN-12-2217

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Lamb J, Crawford ED, Peck D, Modell JW, Blat IC, Wrobel MJ, et al. The Connectivity Map: Using Gene-Expression Signatures to Connect Small Molecules, Genes, and Disease. Science (2006) 313(5795):1929–35. doi: 10.1126/science.1132939

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, et al. A Next Generation Connectivity Map: L1000 Platform and the First 1,000,000 Profiles. Cell (2017) 171(6):1437–52.e17. doi: 10.1016/j.cell.2017.10.049

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Ying X, Huang A, Xing Y, Lan L, Yi Z, He P. Lycorine Inhibits Breast Cancer Growth and Metastasis via Inducing Apoptosis and Blocking Src/FAK-Involved Pathway. Sci China Life Sci (2017) 60(4):417–28. doi: 10.1007/s11427-016-0368-y

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Lim ST. Nuclear FAK: A New Mode of Gene Regulation From Cellular Adhesions. Mol Cells (2013) 36(1):1–6. doi: 10.1007/s10059-013-0139-1

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang J, Xu J, Xing G. Lycorine Inhibits the Growth and Metastasis of Breast Cancer Through the Blockage of STAT3 Signaling Pathway. Acta Biochim Biophys Sin (Shanghai) (2017) 49(9):771–9. doi: 10.1093/abbs/gmx076

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Ardalani H, Avan A, Ghayour-Mobarhan M. Podophyllotoxin: A Novel Potential Natural Anticancer Agent. Avicenna J Phytomed (2017) 7(4):285–94.

PubMed Abstract | Google Scholar

41. Gordaliza M, et al. Podophyllotoxin: Distribution, Sources, Applications and New Cytotoxic Derivatives. Toxicon (2004) 44(4):441–59. doi: 10.1016/j.toxicon.2004.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chattopadhyay S, Bisaria VS, Panda AK, Srivastava AK. Cytotoxicity of In Vitro Produced Podophyllotoxin From Podophyllum Hexandrum on Human Cancer Cell Line. Nat Prod Res (2004) 18(1):51–7. doi: 10.1080/1057563031000122095

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Merzouki A, Buschmann MD, Jean M, Young RS, Liao S, Gal S, et al. Adva-27a, a Novel Podophyllotoxin Derivative Found to be Effective Against Multidrug Resistant Human Cancer Cells. Anticancer Res (2012) 32(10):4423–32.

PubMed Abstract | Google Scholar

44. Chen JY, Tang YA, Li WS, Chiou YC, Shieh JM, Wang YC. A Synthetic Podophyllotoxin Derivative Exerts Anti-Cancer Effects by Inducing Mitotic Arrest and Pro-Apoptotic ER Stress in Lung Cancer Preclinical Models. PloS One (2013) 8(4):e62082. doi: 10.1371/journal.pone.0062082

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Muhammad I, Dunbar DC, Khan SI, Tekwani BL, Bedir E, Takamatsu S, et al. Antiparasitic Alkaloids From Psychotria Klugii. J Nat Prod (2003) 66(7):962–7. doi: 10.1021/np030086k

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Siddiqui S, Firat D, Olshin S. Phase II Study of Emetine (NSC-33669) in the Treatment of Solid Tumors. Cancer Chemother Rep (1973) 57(4):423–8.

PubMed Abstract | Google Scholar

47. Bredholt G, Mannelqvist M, Stefansson IM, Birkeland E, Bo TH, Oyan AM, et al. Tumor Necrosis Is an Important Hallmark of Aggressive Endometrial Cancer and Associates With Hypoxia, Angiogenesis and Inflammation Responses. Oncotarget (2015) 6(37):39676–91. doi: 10.18632/oncotarget.5344

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Zhou YD, Im YP, Mohammed KA, Jones DK, Muhammad I, Dunbar DC, et al. Terpenoid Tetrahydroisoquinoline Alkaloids Emetine, Klugine, and Isocephaeline Inhibit the Activation of Hypoxia-Inducible Factor-1 in Breast Tumor Cells. J Nat Prod (2005) 68(6):947–50. doi: 10.1021/np050029m

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Barbie DA, Amayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA Interference Reveals That Oncogenic KRAS-Driven Cancers Require TBK1. Nature (2009) 462(7269):108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar