{'fraction_non_alphanumeric': 0.028923126433082584, 'fraction_numerical': 0.013234643388712337, 'mean_word_length': 5.519800681877786, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Immuno-Oncology Integrative Networks: Elucidating the Influences of Osteosarcoma Phenotypes

In vivo and in vitro functional phenotyping characterization was recently obtained with reference to an experimental pan-cancer study of 22 osteosarcoma (OS) cell lines. Here, differentially expressed gene (DEG) profiles were recomputed from the publicly available data to conduct network inference on the immune system regulatory activity across the characterized OS phenotypes. Based on such DEG profiles, and for each phenotype that was analyzed, we obtained coexpression networks and bio-annotations for them. Then, we described the immune-modulated influences in phenotype-specific networks’ integrating pathway, transcription factor, and microRNA regulations. Overall, this approach seems suitable for representing heterogeneity in OS tumorigenesis.


Introduction
A current research direction in cancer therapy involves the search for synergistic combinations of cancer immunotherapies aimed to modulate the immune system response. The rise of the immune-oncology field 1 finds supports in the opportunities offered by next-generation therapeutic solutions such as novel molecularly targeted agents 2 and multimodal drug combinations. These latter solutions include not only small molecules 3 but also epigenetic drug combinations and combined epigenetic and targeted therapies. 4 Intuitively, a combination of agents may be expected to yield superior effects than those obtained by each agent alone, unless negative effects occur. Although parameters measurable at a clinical level allow the possibility to control such outcomes, it is much harder to assess the presence and intensity, thus the influence, of both co-inhibitory and co-stimulatory effects due to synergistic therapies.
Osteosarcoma (OS) represents the most common primary bone tumor in children and young adults and the second highest cause of cancer-related death in this age group. Despite aggressive chemotherapy, disease-free survival has not improved significantly in the past 20 years, and 50% of the patients subsequently develop fatal pulmonary metastasis. 5 Therapies are needed, and a relevant role is played by the bone immunobiology, requiring integrative models that combine genetics, genomics and systems medicine approaches. 6 Here, we reused data from evidence generated by both in vivo and in vitro characterization experiments on a panel of 22 OS cell lines that were phenotypically differentiated. 7 In particular, we are interested in detecting phenotype-specific regulatory influences of immunosystem-related processes. To allow for a thorough assessment of possible differential features, we propose a network inference approach centered on gene coexpression dynamics associated with regulation drivers such as transcription factors and microR-NAs underlying the immune-related influences.

Experimental data
Publicly available OS evidence (GSE28425) were recently derived from 22  . They all were used to compute differential expression in multiple phenotypes. Differentially expressed genes (DEGs) were grouped according to the characteristics of the cell lines such as tumorigenic and nontumorigenic, colony and non-colony forming, invasive and noninvasive, and finally, proliferative and nonproliferative. The RPMI-1640 medium (Lonza, Basel, Switzerland) was used for culturing cell lines, whereas short tandem repeats DNA fingerprinting (Promega, Madison, WI, USA) were used for the verification of cell line identity, and then compared with the EuroBoNet 8 and ATCC profiles. In vivo tumorigenicity of the cell lines was checked by injecting subcutaneously into locally bred NOD/SCID IL2R-gamma-0 mice 1 × 10 6 cells in 100-μL serum-free RPMI-1640. Mice were killed when tumors reached a size of 1000 mm 3 . Tumor growth was assessed weekly, and only for six months for the lowly tumorigenic cell lines. Colony-forming assays were plated out in 12-well suspension plate of CELLSTAR (VWR International, West Chester, PA, USA) using standard protocol and GelCount system (Oxford Optronix, Oxford, England) 2 Cancer Informatics was used to quantify colonies larger than 50 μm. About 25 000 cells (in duplicates) in RPMI-1640 were plated, containing fetal bovine serum (FBS), into 24-well invasion chambers with uncoated or Matrigel-coated membranes, 8.0-μm pores (BD Biosciences, Franklin Lakes, NJ, USA). About 10% FBS in RPMI-1640 in lower compartment was used as chemoattractant. The numbers of cells were manually counted after taking 9 pictures of each well/membrane stained with Hemacolor (Merck KGaA, Darmstadt, Germany). Live-cell imaging of IncuCyte system of Essen Bioscience (Birmingham, UK) was used to analyze the proliferation of cells. The cells were seeded in quadruplets in 96-well plates, containing 2000 cells per well, and 2-phase contrast photographs per well were taken every second for 5 days. Cell confluence (in %) over time (in hours) presented with the proliferation rates.

Bio-annotations
We performed with BiNGO 9 the Gene Ontology (GO) analysis of molecular functions aimed to report overrepresentation across DEGs. The procedure of Benjamin and Hochberg was used for multiple testing corrections (false discovery rate). Immune-related processes for DEGs were performed using ClueGO 10 (the kappa score was tuned to high significance, and Bonferroni correction applied). Target gene-miRNA interactions were first extracted from the DEG lists using miRTarBase 11 and TargetScan. 12 Then the associated pathways were then extracted from NCI-PID (the Pathway Interaction Database) [13][14][15] and finally integrated in the delivered association map to create a reference annotation context for the interactions.

Profiling
We analyzed the profiles of DEG in the OS phenotypes. The significant detections for tumorigenic, invasive, colony-forming, and proliferative phenotypes were, respectively, 124, 102, 131, and 63 DEGs, with the log(FC) considered at a cutoff of ±2. Profiling was done with limma (https://bioconductor.org/ packages/release/bioc/html/limma.html, from Bioconductor in R). Overall, only 4 genes, ie, ACTG2, KISS1, NPPB, and KRT17, showed negative regulation in all phenotypes. Supplementary file S1 reports the entire profiles obtained for phenotype.
Among these 4 listed genes, only 2 share some functional characteristics. The first gene, ACTG2, is an actin involved in a mediator of cell motility and in the maintenance of the cytoskeleton. KISS1, a gene known to suppress metastasis in some cancers (melanoma, breast cancer), is also putatively involved in cytoskeletal reorganization and cell adhesion, inhibiting invasion. The role of the other 2 shared genes is less clear with specific reference to OS is less clear. Concerning NPPB, member of the natriuretic peptide family, note that mutations have been reported in association with osteoporosis, whereas KRT17 encodes the type I intermediate filament chain keratin 17 which is involved, among other functions, in tissue repair ( Figure 1).

Bio-annotations
The DEGs in all phenotypes showed common functional aspects related to various binding activities. A few primary annotations emerged from the GO-annotated molecular functions, here reported concisely (for details, see Supplementary file S2). The tumorigenic phenotype involved insulinlike growth factor binding and extracellular matrix (ECM) binding (significant P values). The invasive phenotype involved platelet-derived growth factor (PDGF) binding, plus other functions seen with the tumorigenic phenotype, such as structure molecule activity. The colony-forming phenotype involved both oxidoreductase and antioxidant activities. The proliferative phenotype involved interleukin 1 receptor binding and both oxidoreductase and kynureninase activities. Of interest is also the possible overlap of annotated terms across phenotypic profiles. For instance, tumorigenic, invasive, and colony-forming phenotypes shared molecular functions related to ECM structural constituents containing different members of the collagen family. The PDGF-binding molecular function was instead shared between colony-forming, invasive, and proliferative phenotypes. Among the major molecular functions shared between any pair of phenotypes, there were functions such as receptor binding, carbohydrate binding, kynureninase activity, and diuretic hormone activity. The absence of common molecular functions between proliferative and other phenotypes, except colony forming, indicates a certain specificity of mechanisms for proliferative phenotype ( Figure 2).

Immune-related processes
The 4 networks displayed in Figure 3 emphasize a variety of GO enrichments obtained by ClueGO. 10 These network configurations are obtained directly from the annotation tool and following the enabled features. The views shown here are called "gene distribution views," following the ClueGO-CluePedia networks of terms. Genes can be seen at nodes, in particular relatively small or big hubs for enriched pathway terms. The size of the nodes refers to significance through P values and the links indicate gene participation to the pathway terms. The latter are clustered and colored to emphasize similarly annotated genes. Regarding the selection criteria, all DEGs in each phenotype contributed to the enrichments and in particular for the immune-related terms. In tumorigenic phenotype, both negative regulation of erythrocyte differentiation and regulation of neutrophil migration appear as relatively large communities, ie, interconnected nodes (annotated terms) cross-linked by gene sets. A minor group is formed around the type 1 interferon Sharma and Capobianco 3 signaling pathway. The invasive phenotype presents 3 different communities. One is enriched for leukocyte migration involved in inflammatory response. Another is enriched for neutrophil chemotaxis, and a third one for neutrophil migration.
A relatively minor group consists of regulation of myeloid leukocyte differentiation with a prominent role played by CDK6 and RBP1. The former (http://www.genecards.org) is a member of the cyclin-dependent protein kinase family and has been shown to phosphorylate, thus regulating the activity of tumor suppressor protein Rb (the expression of this gene is upregulated in some types of cancers). The latter is a Runx2 coactivator in both U2OS OS cell models. 7 The colony-forming landscape is more interconnected than the previous ones. Regulation of myeloid leukocyte differentiation represents the largest community centered on MYC, well-known multifunctional, nuclear phosphoprotein playing a role in cell cycle progression, apoptosis, and cellular transformation. By functioning as a transcription factor, it contributes to the regulation of leukocyte differentiation induced by a large community centered on BMP4. This is a gene whose encoded protein is a secreted ligand of the transforming growth factor β superfamily of proteins possibly involved in human cancers. Other minor terms appear annotated, but with marginal impacts.
The proliferative phenotype appears with a sparse landscape of annotated terms, being these however, very important. For instance, osteoclast differentiation involves MAFB, a gene which acts as a transcriptional activator or repressor involved either as an oncogene or as a tumor suppressor; FAM20C, a gene encoding a member of the family of secreted protein kinases playing a role in cell migration and adhesion; NDRG1, whose stress-responsive protein is involved in hormone responses, cell growth, and differentiation and is necessary for p53-mediated caspase activation and apoptosis. At the right side of the network configuration, it appears as major hub IL1B: the protein encoded by this gene is a member of the interleukin 1 cytokine family produced by activated macrophages and important mediator of the inflammatory response, thus involved in a variety of cellular activities, including cell proliferation, differentiation, and apoptosis. The tumor microenvironment (TME) signs are thus quite visible, particularly in the proliferative phenotype.

Coregulation Networks
General pathway-driven influences Figure 4 shows network integrative configurations building phenotype-specific contexts for functional assessment. In particular, the structure of such networks consists in co-expressed DEGs, whereas the functional relevance appears from the association of regulators such as transcription factors and miRNAs when target protein-coding genes are considered among the displayed DEGs. Mutations have also been mapped onto such configurations; likewise, the pathway terms are shown according to the enrichment by the network genes. MiRNA-gene target interactions for DE miRNAs were extracted from miRTarBase 11 and TargetScan. 12 In particular, miRTarBase (http://mirtarbase. mbc.nctu.edu.tw/) contains only experimentally validated miRNA-gene target interactions using reporter assays, Western blots, and CLIP-seq, whereas TargetScan (http://www.targetscan.org/vert_71/) matches seed region of the miRNA with 8mer, 7mer, and 6mer sites for predicting gene-miRNA target interactions. When searching with a gene symbol, the results appear by the different transcripts, classified by prevalence, and for each transcript, the sites with higher/lower probability of targeting by miRNAs appear too. This probability is estimated by including all the algorithms and parameters for each miRNA candidate. For instance, context++ is a model to rank miRNA target predictions in a conservative way (not many alternative 3′-UTR untranslated region isoforms). Therefore, an unbiased confidence score is built to assess miRNA targeting efficiency based on 14 different features useful to reconstruct gene-miRNA target networks. The gene-pathway interactions which are displayed (regions marked in light blue) are associated exclusively with cancer annotations retrieved from the PID 13 (https://github. com/NCIP/pathway-interaction-database, now migrated to NDeX 14,15 ; http://www.home.ndexbio.org/about-ndex/). Such experimentally verified associations have been mapped onto the gene-miRNA target networks.
Overall, pathway landscapes revealed commonalities across OS phenotypes due to the sharing of 32 cancer-related pathways (see Supplementary file S3). These comprise DEGs in part exerting distinct influences and in part resulting share between phenotypes. The best enriched pathways involve β 1-integrin and β 3-integrin cell surface interactions, which are both syndecan-1-mediated signaling events and targets of C-MYC transcriptional repression that contain genes from the col and cyclin D families, MYC and PLAU. Notably, cyclindependent kinases (CDKs) influence cell cycle regulation,

6
Cancer Informatics which is notoriously aberrant in human cancers, thus leading to uncontrolled proliferation. Conversely, the consequently altered activation of CDKs offers through their inhibitors both reasons and opportunities for newly proposed anticancer drugs. On one hand, due to little specificity toward single CDKs, pan-CDK inhibitors have been introduced, such as alvocidib and seliciclib. 17 On the other hand, CDKs link cell cycle to cell metabolism, thus affecting cancer cells in multiple ways and offering an even stronger rationale for representing targets of therapeutic solutions. 18 These pathways also refer to DEGs with confirmed somatic mutations. Especially, pathways related to the integrin proteins are highly relevant in OS because of the role in controlling interaction of tumor cells with their microenvironment. In particular, integrins mediate between the ECM and the cytoskeleton and they function as detectors of environmental signals, being thus involved across cancer hallmarks, thus in cell differentiation, adhesion, migration, proliferation, and survival (a comprehensive review can be found in the study by Bianconi et al 19 ). Clearly enough, all such processes are useful   followed by invasive and colony-forming phenotypes, followed by proliferative phenotype (bottom network). The Weighted Gene Co-Expression Network Analysis (WGCNA) package was used. 16 Blue labels indicate pathway terms, circles indicate miRNAs, rhomboids and hexagons indicate mutations (missense and coding silent), and red links indicate the presence as a connector of a transcription factor.

9
for promoting invasiveness, acquiring drug resistance, and also shaping the TME in favor of cancer growth and metastasis. 20 Our evidence refer to integrin interactors in these pathways, which directly or indirectly mediate interactions between cells and the ECM. [21][22][23] Also, β 1-integrins are implicated in tumorigenesis and chemokine stimulation by the CXCR4 receptor mediating cellular migration in OS. 24 Note that the presence of DEGs referred to ITGA11 and RHO proteins involved in CXCR4-mediated signaling events that may indirectly influence chemokines, leading to irregular functioning and affecting cell surface interactions. Most of the signaling pathways that were mainly mediated by the Hedgehog (HIF-1-α) transcription factor, syndecans, interleukins, and mTOR, were also shared between phenotypes. In particular, it is considerable that the relevance of targeting the Hedgehog signaling pathway in OS is due to its promigratory effects observed in both OS cell lines and primary human specimens, suggesting the design of inhibitors to reduce cell proliferation and tumor growth. 25 Finally, DEG-driven protein-protein interactions were observed as involved in various cancer-related pathways. More specifically, both tumorigenic and invasive phenotypes shared 30 pathways mainly related to signaling, such as Aurora (closely linked to the mutated FBN1). These are kinases with a role in mitosis and cytokinesis. In particular, by phosphorylating and activating PLK1, it promotes CDK1 activation after DNA damage, and by stabilizing N-MYC, it prevents its proteasomal degradation. Several inhibitors targeting Aurora members are known, whereas others are under investigation. 17 Furthermore, the main pathway shared by the colonyforming phenotypes and the other phenotypes are the transcription factors E2F and p73, plus various signaling pathways. FBXO32 and SGK1 proteins showed involvement in the FoxO family (translating the effects of environmental stimuli in gene expression and acting as tumor suppressor in multiple cancers, particularly through the Wnt/β-catenin pathway 26 ), significantly in the tumorigenic phenotype. Instead, the invasive phenotypes contained genes involved in the PDGF receptor signaling network (see the study by Takagi et al 27 ) and circadian rhythm. Pathways specific to the colony-forming phenotype comprised signaling pathways related to IL2, syndecan-3, and CD40/CD40L, whereas the proliferative phenotype showed nonspecific pathways. The tumorigenic phenotype pathways contained in particular CDK4 regulated by multiple miRNAs, such as highly overexpressed hsa-miR-449a in validated targets of C-MYC transcriptional activation. Regulation of nuclear β-catenin signaling and target gene transcription involved in Notch and Wnt signaling was found regulated by lowly downexpressed hsa-miR-186, hsa-miR217 and hsa-miR-590-3p. Note that these 2 major pathways are involved in OS in a complex way, and especially with Notch; here, the restoration of its signaling activity is a potential therapeutic strategy designed to eliminate tumor-initiating cells by promoting their differentiation. 28,29 Furthermore, the overexpressed hsa-miR-193-3p was found to regulate the highly overexpressed PLAU, related to β-integrin cell surface interactions, fibroblast growth factor, amb2 integrin, osteopontin-mediated, and validated transcriptional targets of the AP1 family members Fra1 and Fra2. An example is the highly downexpressed hsa-miR-142-3p regulating IL1A involved in IL1-mediated signaling pathways.
The invasive phenotype-regulated genes were not directly implicated in cancer pathways except for CCND1, CDK6, and GAS1, involved in validated targets of C-MYC transcriptional repression. Signaling events were mediated by Hedgehog and p73. Validated targets of C-MYC transcriptional activation and repression were connected by CCND2 and CCND1 regulated by overexpressed hsa-miR-646. CCND1 showed involvement in 16 different pathways such as integrin-related kinases and Notch and were regulated by numerous miRNAs. Pathways specific to the colony-forming phenotype included MYC and COL1A2, also involved in validated targets of C-MYC transcriptional activation pathway. The downexpressed hsa-miR-516a-5p, hsa-miR-375 and hsa-mir-125a-5p regulate MYC along with other lowly expressed miRNAs, whereas the overexpressed hsa-miR342-3p and the down-expressed hsa-miR-569 regulate COL1A2. These genes appear to share integrins in angiogenesis, β-integrin-cell surface interactions, Notch signaling, and interleukin-mediated signaling events. In proliferative phenotype, pathways related to validated targets of MYC involve CCND2 and NDRG1, respectively, regulated by numerous DE miRNAs. High DE miRNA regulating these genes included downexpressed hsa-miR-153 and hsa-miR-139-5p in case of activation and overexpressed hsa-miR-182 for repression. Lowly expressed hsa-miR-186 regulates CCND2, TXNIP, CLDN1, and GFRA1 that were involved in many different pathways such as IL3 signaling mediated by STAT5. Signaling events by Ret tyrosine kinase and nectin adhesion pathways were regulated by numerous miRNAs, which indicates potential biomarkers in bone metastatic diseases. 30 Figure 4 reports the pathway landscapes just described.

Discussion
Due to the particularly important microenvironmental properties of bones affecting immunology, targeted immunomodulation therapies seem to represent promising paths toward OS treatment. This is especially true in the presence of genomic instability that might potentially affect the capacity to predict genes relevant to mestastasis. 31 Because epigenetic alterations may precede cancer transformation, network approaches have been adopted also in previous applications focused on epigenetic treatment in response to multidrug resistance. 32,33 Here, instead, we show that through networks evidence is emerging from common dynamics at the pathway level for the examined 10 Cancer Informatics OS phenotypes and from distinct effects. The data were obtained by previous extensive work with experimentally characterized cell lines and considering previous studies. 34,35 In particular, following research focused on miRNAs 36 it appears that influences from these regulators that were identified across OS phenotypes clearly indicated a variety of patterns centered on different hub genes enriching cancer-specific pathway terms.
On one hand, it is clear that from the panel of OS cell lines, a pathway-unifying role is played by integrin proteins controlling TME interactions through the exchange with ECM. Consequently, translational programs to clinic need to address such highly complex tumor-stroma circuit whose centrality for cancer hallmarks is clear, and whose disruption fosters dysregulation toward OS (similar to other cancers) progression. These aspects are also key for therapy calibration toward both signaling molecules and enzymes, namely, Shh and FAK inhibitors. 37 On the other hand, when distinct patterns are searched, of relevance becomes the role of targets of C-MYC transcriptional repression referring to genes from the Col and Cyclin D families, such as MYC and PLAU, which appear especially from the colony-forming and invasive phenotypes. Several regulations from many significant DE miRNAs and relatively few DE transcription factors were here assembled into network configurations, thus establishing their influences in the pathway landscapes of each OS phenotype. This is only a beginning step because it is clear that the need of deciphering the role of regulators such as miRNAs and possibly other noncoding RNAs when these are involved in the TME dynamics and cancer progression is influenced by their altered expression. 38 In conclusion, the OS context here examined suggests that all the changes occurring due to mutation and microenvironmental factors reflect genetic, physiological, and spatiotemporal heterogeneities that call for adaptive therapies. 39 Perhaps, and not surprisingly, the network-driven OS phenotype analyses that we have presented appear as enabling tools for the identification of novel targets for therapy. Assessing coregulation remains quite a difficult inference task, but networks approximate well the underlying dynamics and allow transcription factor and miRNA regulators to be cast within pathway contexts. In particular, networks represent probably the best possible syntheses of such highly complex integrative contexts.