Regulation networks of immunity single cell stats in melanoma
Tumor cell-intrinsic heterogeneity shapes the immune cell infiltration and influence the outcome of immunotherapy , and various metabolic pathways orchestrate the behavior of tumor-infiltrating immune cells, which related to enhancing of antitumor immunity and immunotherapy . Moreover, distinct CD45 + cell clusters revealed by single-cell RNA-seq associated with clinical outcome of ICB therapies and reflected by different identifiers . Thus, observation of pathway activities in multiple cell stats could offer both intra- and inter-cellular immunity regulation. Herein, we constructed immunity networks based on BioPAX pathways for previously defined CD45 + cell clusters  and combined them into regulation networks (Fig. 1, Figure S2, Method Details). We fitted power law models for degree distributions of cell cluster specific networks and regulation networks (Figure S3, Method Details), and the R squares greater than 0.88 all over these models. Immune networks and regulation networks constructed by our pipeline are scale free and similar to the real biological networks.
In sum, seven hundred and ninety genes and 3048 interactions, including 26 receptors, 26 ligands, and 47 TFs were recognized by the regulation networks. Notably, IL-2 receptor is T cell stimulating cytokine , and it not only drives the expansion of T cells and the contraction phase of immune response , but also has an effect on cancer stem cells [31, 32]. Most impotently, low dose IL2 combined with other immunotherapy demonstrates benefit in patients with metastatic melanoma . In C10, a memory T cell cluster, IL2 controls expression of KLRK1 which controls state change of ITGAM through CD247(in complex with HLA-C, HLA-E and HLA-G in all 11 cell clusters). IL2 controlling ITGAM, which in complex with CD14 (both are cell surface markers, indicted that our framework efficiently highlighted regulation flows of immunity system in distinct and common immune cell stats.
Divergence regulation of overlap genes were revealed by network comparative analysis among multiple networks
Since our study highlighted different interaction flows of cluster specific networks, we accessed expression of genes and interaction pairs in distinct immune cell stats. There are several genes (86) such as major histocompatibility complex (HLA-C, HLA-DMA, HLA-DMB, HLA-DRA, HLA-E, HLA-G and HLA-H), CD247, cyclin dependent kinase (CDK11B), casein kinase 2 (CSNK2A1 and CSNK2B), FGFR1 oncogene partner (FGFR1OP) etc. that were activated in all 11 predefined T cell stats (Fig. 2A). However, consensus interactions among those common genes are partly observed. On the contrary, obvious inconsistent interaction pattern among common genes were discovered in B cell clusters (C1, C2), myeloid clusters (C3, C4) and CD8 T cell clusters (C6, C9) (, Fig. 2B, Figure S4). Even though identifiers of distinct immune cell stat differ from the others, there are common genes, which exercise functions through different flows, that were activated in multi-networks. For instance, mechanistic target of rapamycin kinase (MTOR) state change was controlled by HLA-G in C9 but not in C6, and MTOR controls state change of HSPA1B in C9, TOP1 and PRNP in C6 (Fig. 2C). STAT2 interacts with different genes in C6 (HLA-G) and C9 (MTOR, PIK3CD). Moreover, a specific duplex interaction was observed in C9 (TRAC controls state change of CD247, CD247 in-complex-with TRAC, Fig. 2C). These results provide further evidence for our explanation that common genes from different cell stats can active variety ways in these cells.
Interleukin-2 (IL-2) antigen stimulate memory CD8(+) T cells production, and high relative IL-2 production in T cells of melanoma tend to perform memory CD8(+) T cells phenotype and superior proliferative capacity compared to cells with low IL-2 production . In our research, IL2 and LIME1 specifically activated in a memory T cell cluster (C10). We also calculated different expression of IL2 in C10. IL2 gene significantly expressed higher in cells from responders than cells from non-responders of on-treatment patients. Furthermore, several genes activated distinctly in different networks. For instance, PVRIG specialized in lymphocytes(C5), AGER, BHLHB9, CDK3, GNG8, IL11RA, PKIA, USP50 in regulatory T cells(C7). ARHGAP19, FNIP1 in exhausted/HS CD8 T cells (C9) and DECR2, SESTD1, ZNF775 in exhausted CD8 T cells (C6) etc. Our framework identified different performance of immune related genes in cell cluster specific networks, which may lead a functional nuance.
Consensus functions across immune cell regulation networks
To investigate functional relevance of cell cluster specific networks, we employed online pipeline metascape (http://metascape.org/gp/index.html) for enrichment analysis. We found that gene sets from plural networks enriched in consensus functions, especially in metabolism of RNA, antigen processing and presentation, cell cycle and herpes simplex infection, regardless the different activation flows they presented (Fig. 3). Significantly, all networks partly enriched in several sub-terms of vital functions like T cell activation, cell cycle, cellular responses to stress, apoptotic signaling pathway, cytokine-mediate signaling pathway, herpes simplex infection, metabolism of RNA and so on [36–40]. As for some functions, such as negative regulation of immune system process, are enriched by all networks but in different sub-terms. In addition, exhausted CD8 T cell cluster (C6) and lymphocytes exhausted/cell-cycle (C11) cluster significantly enriched in regulation of complement activation overall situation, yet regulation of complement activation was not enriched by monocytes/macrophage cluster(C3), cytotoxicity (lymphocytes) cluster (C8) and exhausted/HS CD8 T cell cluster (C9). Several sub-terms of negative regulation of immune system process enriched by different networks despite all eleven networks enriched in negative regulation of immune system process by a significant level.
A core pathway: antigen processing and presentation
We discovered that shared topology of 11 networks covered two major histocompatibility complex union (HLA-C, HLA-E, HLA-G; HLA-DMA, HLA-DRA, HLA-DMB), one ribosomal protein union (RPL10A, RPL9P7, RPL39, RPL12, RPL36A, RPS29, RPS9), one mitochondrially encoded NADH: ubiquinone oxidoreductase core subunit (MT-ND2, MT-ND4L, MT-ND5), one NME/NM23 nucleoside diphosphate kinase union (NME1, NME1-NME2) and one proteasome and proteasome activator subunit (PSME1, PSMB8, PSMB9) (Fig. 4A). Enrichment analysis of these shared topologic structure showed concurrent functions like regulation of expression of SLITs and ROBOs, antigen processing and presentation of exogenous peptide antigen, antigen processing and presentation peptide antigen assembly with MHC class II protein complex and neutrophil deregulation (Fig. 4B, 4C). Especially, subunits of the shared structures lie in two flows of antigen processing and presentation, MHC I and II pathways (, Fig. 4D), which play upstream roles of CD8 T cell killing target cells, regulation of NK cell activity (MHC I) and CD4 T cell cytokine production and activation of other immune cells (MHC II). According to Gene Ontology  enrichment analysis executed by metascape, three antigen related functions are connected from each other, and the most significant enrichment is antigen processing and presentation.
Key genes from networks were related with immunotherapy response at single cell level
Prior knowledge showed that a large number of genes appeared relevance with not only cancer occurrence but also treatment response, and these appearances came up in the single cell level as well [12, 21, 43–45]. Expectably, there are a batch of genes in our regulation networks that differently expressed between responders and non-responders in both period of treatment (pre- and on-treatment) regardless of immune cell stats (Fig. 5A). Functional analysis showed that genes expressed higher in non-responders from on-treatment samples were enriched in cytokine-mediated singling pathway, transmembrane receptor protein tyrosine kinase signaling pathway and so on, yet different expression genes form pre-treatment samples located in immune response-activating signal transduction, response to toxic substance etc. Interestingly, several ligands and receptors differently expressed between responders and non-responders in pre- or (and) on-treatment samples, especially for a higher-in-non-responder set (CCL7, CCL18, MRC2, FPR2, LILRA3, SIRPA and CD14). Moreover, we detected a ligand, SIRPA, significantly expressed higher in non-responders and in pre-treatment samples. Consequently, SIRPA may play a resistance role in immunotherapy [46–48], and a triple element, which consisted with CAV1, SIRPA and CD14 and expressed higher in non-responder, may potentially support drug antagonism([46–51], Fig. 5B-C).
Although we detected dysregulated genes in the single cell level, it still needs further inspection to uncover transcriptome changes in immune cell stats. Therefore, we next assessed differently expressed genes from cluster specific networks in corresponding CD45 + immune cell stats (Fig. 5D, Figure S5-6). Our research suggested that ligand integrin subunit alpha M (ITGAM) was significantly upregulated in C3 from non-responders of on-treatment patients, in C1 and C10 from non-responders of pre-treatment patients, as well as in C8 from responders of pre-treatment patients. Regulating cell fate, ITGAM was dysregulated in multiple cell type (T cells, B cells and so on), which may contribute to tumor cell survival in immunotherapy . ligand NTRK1, which expressed higher in pre-treatment non-responders, was dysregulated in C10, and IL2, a type I cytokine which can be associated with durable regression in metastatic melanoma and renal cell carcinoma , also was dysregulated in and only in C10 but expressed higher in on-treatment responders. Two major histocompatibility complex, HLA-G and HLA-H, expressed higher in C2 from pre-treatment non-responders. Another major histocompatibility complex, HLA-DMB, expressed higher in both C2 and C10 in non-responders from on-treatment samples. Other famous genes, such as tumor necrosis factor receptor superfamily member 9 (TNFRSF9) and CD247 ligands, was also upregulated in C3 from (non-)responders of on-treatment patients, and another tumor necrosis factor receptor superfamily member (TNFRSF12A) was dysregulated in both C3 and C9. Thus, we have reason to suppose that some ligands and receptors may propose a new insight to understand drug sensitivity and resistance, as well as immunotherapy response.
Ligands and receptors of regulation networks showed robust selective power in immunotherapy response
Some receptors can mediate functions of immune cells through distinct signaling pathways . Changing of these receptors and corresponding ligands may lead unexpected immunotherapy outcomes. Our research also observed massive change of gene expression of these proteins. Hence, we applied logistics regression model to select contribution features from ligands and receptors of regulation networks and access immunotherapy precision of featured gene sets. In the test data sets, medium AUC of 10000 random is 0.7143, and most of the test AUCs are among 0.65 and 0.8 (Fig. 6A). Eventually, we identified 17 genes sets that associated with immunotherapy response. In an independent validation, GSE35640, all 17 gene sets performed AUC greater than 0.7 (Fig. 6B), which suggests robust selective power of ligands and receptors of regulation networks for immunotherapy response. This completed key gene analysis by providing gene sets and their scores which able to construct classification of therapy response.
Several TNF receptor superfamily members are integral to the immune-response regulation by enhancing T-cell growth and dendritic-cell function. These proteins related to modulation of cellular functions, proliferation, survival or deaths . Moreover, TNF receptors controlling TNF receptor signaling, which play its role in inflammation and cell death, could determine the cellular fate . Especially, TNFRSF12A (also known as TWEAK receptor, Fn14, or CD266) correlated with integrin β3 expression, which drives Glut3 expression, associated with clinical outcome and tend to be responsible for inducing cachexia in tumors [57, 58]. We identified one TNF receptor leaded gene set, consisting with TNF receptor superfamily member 9/25/12A (TNFRSF9, TNFRSF25 and TNFRSF12A), sphingosine-1-phosphate receptor 1 (S1PR1), C-C motif chemokine ligand 3 (CCL3), caveolin 1 (CAV1) and major histocompatibility complex, class I, G (HLA-G), that predicted immunotherapy response with AUC up to 0.7433 in independent validation and 0.8824 in GSE12575 (Fig. 6C). Besides, CAV1 catalysis precedes of resistance related gene CD14 [51, 59], as well as supported a firm set predicting immunotherapy response (0.7045 in GSE35640, 0.871 in GSE120575) with CD247, HLA-C, HLA-E, ITGAM and CD14. Thus, TNFRSF12A, which works for 11 of 17 selected gene sets, cooperating with CAV1 and other ligands and receptors are key factors in immunity regulation and flexible immunotherapy response (Fig. 6D).
Additionally, other gene sets with major histocompatibility complex (HLA-C, HLA-E, HLA-G), TNF receptors, ITGAM, CD247 and CD14 performed acceptable precision with AUC around 0.72 in the independent datasets and above 0.85 in GSE12575 as well. Furthermore, ITGAM, which regulate cell fate , also works 11 of 17 selected genes sets and HLA-C works in 9 Fig. 6E). These genes not only activate in majority immune cells but also have a quality for prediction of immunotherapy response. Our results indicate that using ligands and receptors in the regulation networks to train decision models could provide a brand-new view for immunotherapy response, and these models could be potential guidance for precision medicine.