The landscape of cancer rewired GPCR signaling axes

Abstract

We explored the dysregulation of GPCR ligand signaling systems in cancer transcriptomics datasets to uncover new therapeutics opportunities in oncology. We derived an interaction network of receptors with ligands and their biosynthetic enzymes, which revealed that multiple GPCRs are differentially regulated together with their upstream partners across cancer subtypes. We showed that biosynthetic pathway enrichment from enzyme expression recapitulated pathway activity signatures from metabolomics datasets, providing valuable surrogate information for GPCRs responding to organic ligands. We found that several GPCRs signaling components were significantly associated with patient survival in a cancer type-specific fashion. The expression of both receptor-ligand (or enzymes) partners improved patient stratification, suggesting a synergistic role for the activation of GPCR networks in modulating cancer phenotypes. Remarkably, we identified many such axes across several cancer molecular subtypes, including many pairs involving receptor-biosynthetic enzymes for neurotransmitters. We found that GPCRs from these actionable axes, including e.g., muscarinic, adenosine, 5-hydroxytryptamine and chemokine receptors, are the targets of multiple drugs displaying anti-growth effects in large-scale, cancer cell drug screens. We have made the results generated in this study freely available through a webapp (gpcrcanceraxes.bioinfolab.sns.it).

Figure 1. GPCR extracellular signaling network and overview of the datasets considered

a) workflow of the procedure to generate enzyme-receptor interacting pairs; b) Sunburst charts providing an overview of the number of liganded and orphan G-protein-coupled receptors (GPCRs) considered; Classification of the liganded GPCRs based on their ligand type, i.e., either peptide, organic, or a combination of both; number of enzymes that are either currently included or absent within the Reactome pathways. The included ligands are additionally subcategorized based on the frequency of their types and the frequency of the Reactome pathway domains with which they are linked. For enzymes, only the second layer of this information is shown; c) The network of GPCRs and cognate ligands and biosynthetic enzymes; d) Funnel plot exhibiting the frequency of differentially expressed (TCGA vs. GTEx, BH corrected Padj<0.01) GPCRs, ligands, and enzymes corresponding to each whole tissue. e) Stacked bar-plot displaying the number of significantly enriched ligand-associated pathways within Reactome sub-domains: Signaling (‘R-HSA-372790: signaling by GPCR’) and Metabolism (‘R-HSA-1430728:metabolism’). Each bar represents the relative proportions of different categories as percentages of the whole, with respect to each tissue. An enriched pathway refers to a pathway enriched in TCGA with a GSEA enrichment score>0 and BH corrected Padj < 0.01.

Figure 2. Peptide Ligand-Receptor co-regulation

a) Heatmap (center panel) displaying the co-differential regulation for Receptor-Ligand pairs across different TCGA subtypes (color-coded at the top row). Darker red represents both receptor and ligand significantly co-up regulated in TCGA (i.e., LFC>1 and Padj<0.01) and darker blue represents both receptor and ligand significantly co-down regulated in TCGA (i.e., LFC<1 and Padj<0.01). Paler colors represent either of the receptor-ligand as significantly DE (i.e., |LFC|>1 for both but Padj<0.01 for only one of these). White cells indicate anti-regulation or no fold change at all in at least one of them. Only those pairs which are affected in at least 25% of TCGA subtypes are displayed. A dashed line separating the two clusters, created using Hierarchical clustering, is shown in the middle. Hierarchical clustering was performed using the ‘ward’ method to identify two homogeneous clusters by minimizing within-cluster variance based on the sum of squared differences of feature values i.e. scores∈[−2,+2] assigned to each pair. Heatmap (left panel) uses color codes to display the ligand’s mechanism of action in ‘Action’, G Protein-coupling associated with the GPCRs and GPCR Family. Heatmap (right panel) displays the information contained in the left panel as concise bar plots for each of the two clusters. For each cluster barplots indicate G-protein coupling preferences as well as GPCR families represented as count of the number of recepors; b) The scatter plot displays a GPCR-centric view of the heatmap in (a) based on tissue-wise or ligand-wise similarity. Here, ligand-wise similarity for each GPCR is the average of the pairwise Euclidean distances between co-differentially expressed (co-DE) profiles across its ligands (i.e., vertically), while tissue-wise similarity is a similar metric across tissues (i.e., horizontally). The size of the markers is proportional to the number of ligands, and the colors correspond to the average difference in the number of co-up and co-down profiles. GPCRs located in the top-right corner are more diverse in terms of their ligand interactions; c) Diverse co-DE profiles are observed with respect to different ligands and tissues corresponding to the same GPCR. For example, CXCR3 and its corresponding ligands exhibit an almost exclusively co-upregulated profile irrespective of tissue type, while CALCRL and its ligands display varied co-regulation profiles, implying a varying interaction with its ligands in different tissue (or cancer) types.

Figure 3:. Biosynthetic pathway enrichment

(a) Heatmap of GSEA pathway enrichment analysis for differentially expressed (DE) genes in cancer (TCGA) over normal (GTEx) samples, considering Reactome pathways related to signaling and metabolism. Red cells indicate pathways enriched in TCGA, blue cells indicate pathways enriched in GTEx, and white cells indicate no enrichment. Cells marked with ‘*’ correspond to a statistically significant enrichment (Padj < 0.01). (b) Concordant pathway enrichment analysis considering both DE genes (Enrichment Score, ES) and metabolites (Pathway Abundance Score, PAS). A score of 1 (red) indicates co enrichment in cancer (ES & PAS > 0), and a score of −1 (blue) indicates co enrichment in normal samples (ES & PAS < 0). (c) Bar plot displaying numerical values for ES and PAS for the ‘Purine catabolism’ pathway in three different tissues: Breast, Bladder, and Prostate. (d) Heatmap of DE genes for the ‘Purine catabolism’ pathway in Breast cancer. Red cells indicate upregulated genes in TCGA, blue cells indicate downregulated genes in TCGA, and cells marked with ‘*’ correspond to a statistically significant differential expression (Padj < 0.01). (e) A functional interaction network between genes (ovals) and metabolites (diamonds) in the ‘Purine catabolism’ pathway in Breast cancer. Red nodes indicate upregulated components, blue nodes indicate downregulated components, and green nodes indicate no information available. The network shows that over-activation of the pathway is contributed by over-expression of both genes and metabolites; (f) Sunburst chart (left panel) displaying the distribution of co-differentially regulated Receptor-Enzyme pairs across two different clusters similar to receptor-ligand co-regulation analysis (Supplementary SX). Darker red represents the number of Receptor-Enzyme pairs significantly co-up regulated in TCGA (i.e., LFC>1 and Padj<0.01) and darker blue represents the number of Receptor-Enzyme pairs co-down regulated in TCGA (i.e., LFC<1 and Padj<0.01). Paler colors represent the number of pairs wherein either of the receptor-enzyme is significantly DE (i.e., |LFC|>1 for both but Padj<0.01 for only one of these). Grey indicates pairs with anti-regulation or no fold change at all in at least one of them. Only those pairs which are affected in at least 25% of TCGA subtypes are utilized for the clustering. For each cluster, barplots (center and right) indicate G-protein coupling preferences as well as GPCR families represented as count of the number of receptors

Figure 4. Association of GPCR-peptide ligands axes to survival

(a) Survival association of G protein-coupled receptor (GPCR) components in various cancer types. The funnel plot shows the number of significant instances of individual components, such as receptors, ligands, and enzymes, whose expression values are associated with patient survival across various cancer types. A total of 302 unique significant receptor instances (log rank p-value < 0.05, FDR<0.1) were identified in 11 cancer tissues, 103 ligands in 10 cancers, and 71 enzymes in 10 cancers. (b) Scatterplots for individual receptor instances, enzymes and ligands associated with patient survival in various cancer types. (logrank-p<0.05, FDR<0.1). x-axis displays the number of tissues in which the genes are associated with lower survival (inverted blue triangles) and, vice-versa, y-axis displays the number of tissues in which the genes are associated with higher survival (red triangles). Only the genes which are significantly associated with survival in at-least two tissues are shown. c) The bubble plot shows the correlation between the combined-expression levels of GPCR-ligand pairs and patient survival across various cancer types. The pairs with a log rank p-value <0.05 and also lower than log rank p-values for individual GPCR/ligand are displayed. Bubble color is proportional to HR: i.e., HR>1 : High expression is correlated with high survival (red); HR<1: High expression is correlated with poor survival (blue). Bubble diameters are proportional to the -log10(log-rank p-value). Green highlighted bubbles represent the most significant instances (sample sizes>5, FDR<0.1). (d) The expression values for the CALCR-CALCB axes in breast cancer (right) and corresponding Kaplan-Meier curve for survival analysis of patients stratified based on the individual as well as combined receptor/ligand expression;(e) Scatter-plot showing with differentially co-expressed GPCR-ligand pairs (upper triangle=co-up regulation; lower triangle=co-down regulation) that were significantly associated with patient survival (red=higher survival; blue=lower survival). A total of 24 receptor-ligand pairs from 15 cancer subtypes were identified.

Figure 5. Association of GPCR-enzymes axes to survival

(a) The bubble plot shows the correlation between the combined-expression levels of GPCR-enzyme pairs and patient survival across various cancer types. The pairs with a log rank p-value <0.05 and also lower than log rank p-values for individual GPCR/enzyme are displayed. Further, the pairs with rate-limiting enzymes having less than 100 PubMed references are not shown. Bubble color is proportional to HR: i.e., HR>1 : High expression is correlated with high survival (red); HR<1: High expression is correlated with poor survival (blue). Bubble diameters are proportional to the -log(log-rank p-value). Green highlighted bubbles represent the most significant instances (sample sizes>5, FDR<0.1). (b) The CHRM3-CHAT axis is recurrently correlated with lower survival rates in Esophageal cancer. The KM plots (left) display the risk enhancement achieved when CHRM3-CHAT combined stratification was applied as compared to individual CHRM3/CHAT as evident by log rank p-values. Similarly, the ADORA2B-PNP axis is also consistently associated with lower survival rates in Uterine cancer, as evident from the KM plot on the right.

Figure 6. GPCR ligands with cancer cell line growth-inhibitory capacity

(a) The scatterplot displays 52 GPCR ligands that were tested in 578 cell lines in PRISM drug repurposing resource.13 GPCR ligands (in red) were found to significantly inhibit the growth of cancer cell lines (correlation>0.2, bimodality coefficient>0.35). (b) Venn diagram showing the overlap of GPCRs in the target proteins of these 13 ligands with the pooled list of GPCRs that were found to mediate survival associated GPCR-ligand/GPCR-enzyme axes. ( c) The heatmap displays the 13 shortlisted ligands and their corresponding GPCR targets. The targets marked with ‘*’ are the overlapping GPCRs with significant axes. Heatmap also displays the overall action mode of ligand, more detailed mechanism of action (MOA) and G-protein coupling associated with the GPCR; d) stacked barplot of top 10 mechanisms of actions (MOA) involving the top 10 GPCR drug with the highest growth inhibition potential (ranked according to the negative of the log2fold change in each cell line). Each stack is proportional to the number of cell lines of a given tissue where a drug belonging to that MOA is ranked among the top 10 inhibiting drugs. Stack coloring is tissue specific; e) stacked barplot, normalized for the total cell line-drug pairs count, of top 10 drug targets hit by the top 10 GPCR drug with the highest growth inhibition potential (ranked according to the negative of the log2fold change in each cell line). Each stack is proportional to the number of cell lines of a given tissue where a drug hitting the specific target is ranked among the top 10 inhibiting drugs. Stack coloring is tissue specific.