Unveiling the Link Between Inflammation and Adaptive Immunity in Breast Cancer

Inflammation has been recognized as an important driver in the development and growth of malignancies. Inflammatory signaling in cancer emerges from the combinatorial interaction of several deregulated pathways. Pathway deregulation is often driven by changes in the underlying gene regulatory networks. Confronted with such complex scenario, it can be argued that a closer analysis of the structure of such regulatory networks will shed some light on how gene deregulation led to sustained inflammation in cancer. Here, we inferred an inflammation-associated gene regulatory network from 641 breast cancer and 78 healthy samples. A modular structure analysis of the regulatory network was carried out, revealing a hierarchical modular structure. Modules show significant overrepresentation score p-values for biological processes unveiling a definite association between inflammatory processes and adaptive immunity. Other modules are enriched for T-cell activation, differentiation of CD8⁺ lymphocytes and immune cell migration, thus reinforcing the aforementioned association. These analyses suggest that in breast cancer tumors, the balance between antitumor response and immune tolerance involving CD8⁺ T cells is tipped in favor of the tumor. One possible mechanism is the induction of tolerance and anergization of these cells by persistent antigen exposure.

1. Introduction

Breast cancer is the most frequently diagnosed malignancy in women worldwide (1). Given its prevalence, a great effort has been made to understand the mechanisms that lead to its development including genetic analysis, mutation status for known oncogenes ad tumor suppresors, expression status of known associated receptors like ER and ERBB2 (HER2), as well as transcriptomic assays as exemplified in Cancer Genome Atlas Network (2). Breast cancer tumors appear to be highly heterogeneous in all these aspects, nevertheless a shared set of characteristics are recognizable during tumor development, which are widely known as the hallmarks of cancer (3, 4). In this view, cancerous phenotypes result from a complex of interacting biological processes or pathways that are subverted in favor of tumor survival, growth and invasion.

Inflammation is necessary to maintain homeostasis. Nontheless, alterations in the course of inflammatory response can lead to pathological states, specially when it cannot be resolved and becomes chronic. The effects of inflammation include changes in tissue properties like blood vessel permeability, alteration of extra cellular matrix and the initiation of immune responses that involve recruitment, proliferation and differentiation of innate as well as adaptive immune cells, actions that are coordinated in a combinatorial way by cytokine and chemokine production (5, 6). There is a strikingly similar scenario when we look at tumors.

The role of inflammation in tumor development is far from simple. It is known that acute inflammatory response is an important component in the spontaneous elimination of tumors (7). Also there are examples in which correlation exists between chronic inflammation and some types of cancer like colon, gastric, esophageal, thyroid and breast cancers (8–12). Although it is recognized that chronic inflammation has a role in the tumorigenic process, it is not well understood how or why inflammation is also necessary for tumor maintenance. In particular, breast cancer tumors show significant infiltrates of immune system cells (13). Infiltrates can vary in the amount and specific types of immune cells they contain. Some evidence points to a relation between mixtures of immune cells in breast cancer subtypes and patient outcome in relation to metastasis recurrence and survival time (5, 14). This suggests that inflammatory response in conjunction with the immune system plays an important role on breast cancer maintenance.

Considering the intricate set of relationships that inflammation and immunity play in a complex phenotype such as a breast cancer, an integrative approach capable of capturing at least part of the complexity of biological processes related to inflammation results appealing. In this context, the construction of gene regulatory networks via the association between gene expression levels, offers us a valuable view on how groups of genes are being collectively coordinated.

The modeling of the gene regulatory structure through the use of complex networks (15–17) opens up the use of techniques to interrogate network structure (15), as is the case with network modularity (18, 19). Resulting modules can then be associated biological functions (20). To find the modular structure in a network is a non-closed problem, which has been subject of intense research in network theory (15, 21). Furthermore, the applications of finding modular structure in real-world networks has also been matter of interest in several fields, including social (22) and biological systems (20, 23).

Previous work from our group has shown how network modularity structure is associated to biological functionality in transcription factor networks (20) and distinctive network structure for each of the major breast cancer molecular subtypes (24). Furthermore, by means of an automated analysis of network inference, module detection and enrichment analysis we have been capable to detect specific processes in breast cancer molecular subtypes.

In this work, we approach the entangled nature of the inflammatory process in the maintenance of breast cancer phenotype through the inference of an inflammation-related gene regulatory network and their associated genes (their immediate or first neighbors in network terminology), then finding the modular structure of the network and associated biological processes for each module.

We found a hierarchical modular structure in transcriptional networks associated to inflammatory response, where genes tend to be connected to others with similar differential expression patterns: overexpressed genes are more connected between them, as well as underexpressed ones. Modules of the network are mainly associated to immune system, extracellular matrix and cell adhession.

A comprehensive integrative framework allows us to observe a broader landscape of how inflammation may have influence in the establishment of pathological phenotypes. At the same time, this exploratory approach could help to direct research toward more specific questions about therapeutic options.

2. Materials and Methods

2.1. Data Acquisition

A graphical description of the followed pipeline is depicted in Figure 1. We used a curated database that consists of 719 microarray assays of the platform Affymetrix HGU 133 plus 2 from 641 untreated primary breast cancer samples and 78 healthy mammary tissue samples. The data was obtained from five separate data series deposited in the Gene Expression Omnibus https://www.ncbi.nlm.nih.gov/geo/(25), each data series contais case and control samples (Table 1). Raw data. CEL files were pre-processed and normalized as described in https://github.com/CSB-IG/rnw/tree/master/normalization-preprocessing. The resulting expression matrix was used as input to perform network inference. Differential expression (DE) analysis contrasting tumors to healthy mammary tissue was performed with the normalized expression matrix using Empirical Bayes statistic from Bioconductor's library Limma (30) for the R statistical software.

FIGURE 1

Figure 1. Graphical description of the methods followed in this work. (1) Gene expression datasets from GEO were selected. All datasets contained normal tissue samples and primary, untreated tumor samples. (2) Datasets were preprocessed, merged and normalized to obtain an expression matrix (Rows correspond to genes and columns correspond to samples). (3) Tumor samples expression matrix was used to calculate all statistical dependencies with the MI function between pairs of genes to construct the network. (4) The network was filtered to obtain the top 10,000 interactions ranked by MI value of genes in the inflammatory response process and first neighbors. (5) The inflammation network was analyzed to detect modules via the Infomap algorithm. (6) The resulting modules were tested for functional enrichment of Gene Ontology.

TABLE 1

Table 1. GEO identifier and references for the data used here.

2.2. Network Inference

Since we are interested in the integration of inflammatory response with other biological processes in the breast cancer phenotype, we used a gene regulatory network approach to obtain groups of associated genes. Although numerous algorithms and correlation measures can be used to infer networks based on expression data (31), it can be formally demonstrated that for datasets with more variables than samples and inherently noisy data, the best statistical dependency measure is Shannon's mutual information (MI) since it is capable to detect non-linear dependencies and is not affected by data transformations like normalization (it is reparameterization invariant) (32–34). A network in this context is a mathematical object composed of a colection of nodes that represent genes, and a colection of edges that represent the statistical dependency between pairs of genes.

We used MI calculations implemented in the ARACNe algorithm (34) to calculate pairwise statistical dependency for all gene pairs in the platform (MI threshold set with p ≤ 1) which gives us a completely connected graph (a network where all possible interactions exist) for the whole-genome. Our network contains only the top 10,000 most stringent interactions between genes from the inflammatory response and their first neighbors, which correspond to 0.005% of possible correlated gene pairs. The reasoning behind this is that the genes defining the phenotype must have the strongest statistical dependencies between them, hence are at the top of the ordered list, and also in order to minimize the effect of false positives. Since we are interested in inflammation, to start with a list of genes of interest can help us to discover other genes associated to them at the transcriptional level. This however does not guarantee that we recover all genes in the list but the ones within our imposed threshold.

The network was curated starting from the complete graph using a reference gene list of the inflammatory process obtained from Gene Ontology (Inflammatory response process GO:0006954) with the following procedure: First, we ranked all interactions in descending order (largest to smallest MI value). Second, we searched for the genes sharing the strongest interactions with those in the inflammatory response process. This was done by filtering those interactions where at least one gene is part of the inflammatory response. From these interactions we took the top 10,000 interactions and obtained the names of all the genes involved which comprises the list of inflammatory response and associated genes. Third, because we wanted to recover interactions between inflammation gene first neighbors, we extracted from the network the top 10,000 interactions between the genes of the list (inflammation and first neighbors).

2.3. Hierarchical Module Detection

Gene Regulatory Networks, such as those we inferred, contain associations in the scale of hundreds of genes connected by thousands of interactions. Nevertheless, the connectivity patterns are not uniform through the structure of the network. Highly interconnected groups of genes may be indicative of coordinated expression associated with biological functions, making relevant the identification of such groups. This approach has been previously proved by our group in the context of GRNs on breast cancer subtypes (23). We showed that gene modules are representative of distinct and meaningful biological functions.

Thus, in order to find to find connectivity patterns in our network, here we used Infomap (19), a well-known flow-based information clustering method to determine the modular structure of complex networks. Likewise, we use the expanded version of Infomap to find a finer modular structure over the two-level modules using the hierarchical version of the map equation (35).

Using the hierarchical map equation it possible to exploit the fact that the modules in a network are themselves organized into submodules and sub-submodules which can reveal a richer multilevel organization. This approach has been successfully applied as well in the case of GNRs associated with Her2+ cancer subtype (36).

2.4. Functional Analysis

The number of genes contained in each module can amount to several hundred. Additionally, we have to deal with the fact that individual genes can be annotated for more than one function or pathway. To obtain biological insights from gene sets like these, we use statistical over representation analysis to reduce such large sets of individual gene names to identifiable biological functions (37).

We tested or network modules and submodules for enrichment in Gene Ontology (GO) (38) Biological processes. We used GO because it offers a comprehensive annotation of molecules over a wide range of processes thus serving as a valuable first approach. Overrepresentation was calculated with WEBGESTALT (39). Statistical significance threshold was set at p ≤ 0.05 after Benjamini & Hochberg FDR correction. For each module, we performed Over-Representation Analysis (ORA) based on FDR-corrected hypergeometric test over a category of genes whose functions are annotated in the Gene Ontology Consortium database GO (38) with the R package HTSanalyzeR (40), choosing a significance p ≤ 1 × 10⁻⁵.

3. Results

3.1. Inflammation-Related Network Has a Characteristic Expression Pattern for Each Module

The Inflammation-associated Gene Regulatory Network (IGRN) contains the top 10,000 interactions ordered by MI value; this IGRN has 942 genes with three connected components that contain more than 10 genes (Figures 2A–C and Supplementary File 1). Information about the DE status for each node was mapped to the network. This revealed a definite composition of DE genes for each component.

FIGURE 2

Figure 2. Network obtained from the MI highest interactions between inflammatory process genes and their first neighbors. The network consists of 942 genes and 10,000 interactions (edges). Node color represent the differential expression status compared to healthy mammary tissue. Red: Overexpressed, Blue: Underexpressed, White or pale color means no differential expression (−1≥LogFch ≤ 1). This network consists of many connected components of which three of them (A–C) consists of more than 10 genes. Indicated with (D) are small components of less than 10 genes. By observing the way nodes aggregate in the largest component (A), it is evident that the network has an internal modular structure where differential expression levels seem to cluster with similar differential expression trends. We further explored how this network is organized and which known biological processes are being regulated.

The largest component contains 787 out of the 942 genes in the network. We explored the module structure of this component, which revealed a hierarchical structure of 4 first level and 28 second level modules or submodules (Table 2 and Supplementary File 2). Modules and submodules were labeled based on the highest page-ranked gene (Figure 3A). Hereon, modules and submodules are indicated with a subscript after the name of its highest ranked gene i.e. CCR5_m and CD2_sm.

TABLE 2

Table 2. Modules and submodules in the largest component of the network.

FIGURE 3

Figure 3. Modular structure of the largest component of the network. Twenty one submodules were identified. (A) infomap visualization. Nodes correspond to each detected module. Darkest colors and largest circles show those modules with more genes; link widths correspond to the flow between modules, proportional to the inter-module edges. Node labels refer to genes with highest page rank in the module. (B) inter-module edges preferentially link genes with similar expression patterns. Modules are visualized by differential expression. Red links are inter-module edges linking two overexpressed genes, blue links are analogous but regarding underexpressed genes, gray links show intra-module edges and inter-module edges between no-differentially expressed elements.

When we observe the DE status in individual submodules, it becomes evident how it tends to display characteristic patterns for each submodule. For instance, submodules such as LDB2_sm are integrated by genes with a tendency to underexpression, and IFI44L_sm integrated by genes tending to overexpression (Figure 3B). In the case of EXOC3L2_sm and ZFPM1_sm, all genes have a DE status of no change respect to normal mammary tissue. However, other submodules such as SPARC_sm and CCR5_sm show a mixed DE profile. Coordinated over and underexpression gene sets may account for modulation of opposing or conflicting pathways.

Now, based on the module detection method, a novel feature revealed is the information flow between the aforementioned modules, a measure of how much information is shared between them, which in turn is proportional to the number of inter-module edges, which is the number of links between genes belonging to different modules. The major contribution to this flow is by gene pairs with similar expression patterns: overexpressed genes join two different modules, and concomitantly underexpressed genes are responsible for the union between underexpressed modules. Interestingly enough, those modules with non-differentially expressed genes are connected between them (Figure 3).

The modular expression pattern of the IGRN suggests a coordinated activity of genes whose products (proteins) may be involved in the activation or inactivation of processes. It is not clear –although it may be possible– whether clusters of genes that are not differentially expressed could be related to processes of maintenance of mammary tissue, independent of the disease. The fact that genes included in the modules have a similar expression pattern concurs with the idea that network modules are not only a topological feature of the network, but also reflect a functional role in the definition of phenotypes (20, 23).

3.2. Network Modules Are Function-Specific

Modules are disjoint sets of nodes at the same hierarchical level. We performed overrepresentaion analysis for these gene sets, taking Gene Ontology:Biological Processes as the reference database. Statistically significant enrichment was found, but not for all modules and submodules.

Top enrichment scores for the CCR5_m include immune response, innate immune response, cytokine-mediated signaling pathway, adaptive immune response, also included is inflammatory response (Figure 4). SPARC_m top enrichments include extracellular matrix organization, cell adhesion, extracellular matrix disassembly, collagen catabolic process and angiogenesis. Since the inflammatory response process was enriched only in the CCR5_m we focused on the biological processes and pathways of its submodules.

FIGURE 4

Figure 4. Heatmap of GO Biological process enrichments for CCR5_m submodules. In this representation, each column corresponds to a detected module according to the network structure. Each row is the enriched biological process in the aforementioned modules. The color code represents the enrichment p-value: red color takes account for the most significant values. In this picture can be observed that, despite CCR5_sm, each module has a set of unique enriched processes, which could be related to the specificity of functions by each submodule.

Seven out of eleven CCR5_m submodules showed a statistically significant enrichment for GO Biological processes. Some of these enrichments were shared between submodules, with the largest category overlap between CCR5_sm and CD2_sm (Figure 4). These two submodules have numerous connections and are the two largest in terms of of gene and interaction counts (Figure 2).

3.3. CCR5 Submodule Is Responsible for Inflammatory Cell Recruitment

The largest submodule with 176 genes, is enriched for inflammatory response and other adaptive immune response processes. Genes coding for proteins involved in T lymphocyte signaling such as CD45, CD4 and CD28 which serve as coreceptors necessary for T cell activation. LYN, BTK, PIK3CD and VAV1 participate in the B cell signaling pathway, as well as in the Fcγ and Fcϵ signaling pathways through which receptor mediated endocytosis can be activated allowing antigen presenting cells to capture antigens bound to antibody molecules.

Said submodule contains genes for the receptors Fcγ RIIA and Fcϵ RIγ that participate in receptor mediated endocytosis signaling and recognize IgG and IgE class antibodies, respectively. Alternatively to the FcR pathways, phagocytic cells can produce molecules of the complement system that aid in the phagocytic process. C1q is a 18-mer integrated in a 1:1:1 proportion by the products of C1QA, C1QB and C1QC genes included in this submodule. Although components of the Complement System are produced in the liver, C1q is produced extrahepatically by macrophages and dendritic cells. Apart from its role in the classical complement pathway, C1q can bind to calreticulin exposed in membranes of apoptotic cells (41) and to antibodies that permeate and bind to intracellular antigens (42). This alternative functions give C1q an important role during apoptotic cell clearance by phagocytic cells (43) and even relates to autoimmunity development in certain contexts (44).

To induce effector function in CD4⁺ T cells it is necessary that antigens be presented by professional antigen presenting cells (APC) through MHC class II molecules. To achieve this, externally acquired antigens need to be processed in the phagosome and mounted in MHC-II molecules prior to exposure in the membrane of the APC. This submodule contains the gene HLA-DMB, responsible for the exchange of the clip peptide with antigen in MHC class II and shares numerous edges with the HLA−DRA_sm. This submodule also contains NCF1, NCF2, NCF4, and CYBB, whose products are members of the NADPH oxidase complex isoform NOX2. NOX2 catalyzes the production of reactive oxygen species and proton consumption in phagosomes but is also responsible for a slower acidification of the endosomes, maintaining a higher pH which delays protein degradation by phagosomal proteases (45). Slow protein degradation in endosomes is associated with antigen cross-presentation by MHC class I molecules (46), which allows APCs to activate CD8⁺ T cells who normally attack infected or transformed cells through presentation of internally-only produced antigens.

Genes of the signal transduction process include chemokine receptors and chemokines involved in immune cell recruitment. CCR5 and its ligands CCL4 (also known as Macrophage inflammatory protein-1β), CCL5 (also known as RANTES), and XCL2, are invloved in T cell and monocyte recruitment. CXCR3 ligands CXCL9 (MIG), CXCL10 (IP-10), and CXCL11 (I-TAC) are part of the neighboring submodules CD2_sm and CXCL10_sm. These cytokines are involved in leukocyte (CTLs, NKTs and macrophages) migration. CXCL10 and CXCL11 are also reported to be induced in response to interferon gamma (IFNγ) (47, 48) and are potent attractants for NK andactivated T Cells (49).

The molecules coordinated in this submodule suggest the coregulation of genes that elicit inflammatory cell recruitment mediated by cytokines, phagocytosis, antigen processing and presentation. These functions are at the beginning of adaptive immune response, and complemented by characteristic processes of effector immune cells, particularly cytotoxic CD8⁺ cells. This result will be discussed in the following section.

3.4. CD2 Submodule Contains Genes Characteristic of CD8+ T Cells

This submodule contains 87 genes, and shares most of its GO enriched processes with the CCR5_sm. CD2_sm is highly connected to the CCR5_sm sharing 831 edges. Here are expressed cytokines CXCL9, a chemotactic cytokine for lymphocytes recognized by CXCR3 and CCL8 also known as monocyte chemoattractant protein 2 (MCP-2), recognized by CCR5 which serves as a chemoattractant to monocytes, which later can differentiate to phagocytic cells, and to attract lymphocytes, including T cells.

In this module we observe genes for the components of the T cell receptor, including the α, β and ζ chains as well as the CD3 γ, δ and ϵ chains, which amounts for the entire functional TCR assembly. Coexpressed genes for co-stimulatory molecules CD8A and ICOS are also found. CD2_sm contains IFNG, perforin (PRF1), and granzymes A, B, K, and H (GZMA, GZMB, GZMK, GZMH) genes and EOMES. IFNG is the only member of the interferon family that is found in the network and is a mediator in antiviral, antimicrobial and antitumoral responses, while perforin and granzymes are molecules responsible for the effector cytotoxic function of CD8⁺ and NK cells. EOMES codes for a transcription factor paralogue of T-bet that is induced in CD8⁺ T cells and has been associated to the expression of IFNG, perforin and granzime B in this cells (50).

For CD8⁺ T lymphocytes, antigens must be presented in MHC-I molecules. MHC-I presented antigens generally come from cytosolic degradation of proteins in the cytosol by the proteasome. Normal cells present endogenous peptides this way, but CD8⁺ T cells need to be activated by APCs that present externally acquired antigens processed in a distinct compartment called the phagosome. These external antigens can reach MHC-I molecules, a process known as cross-presentation (51). In an inflammatory environment, stimulation with interferon gamma induces the expression of alternative components of the proteasome modifying the way proteins are cleaved and thus altering the repertoire of presented peptides (52).

As we mentioned in the previous section, in this submodule we observe a funtional relationship with CCR5_sm. This is backed by the fact that CCR5_sm triggers the adaptive immune response, meanwhile CD2_sm is related to effector functions, mainly involved in CD8⁺ cells response, which could imply a temporal relationship between both submodules, starting with CCR5_sm.

3.5. PSMB9 Submodule Is Associated to MHC Class I Antigen Presentation

This submodule contains TAP1 and TAP2, along with PSMB8 and PSMB9 genes. These are functionally related to antigen processing and presentation via MHC class I. Antigens presented by MHC class I are mainly peptides that result from protein degradation in the cytosol by the proteasome. When cells are stimulated by cytokines like IFNG, genes of alternative components of the proteasome PSMB8, PSMB9 and PSMB10 are induced and its protein products replace components of the catalytic center of the proteasome to form what is known as the immunoproteasome (52).

The immunoproteasome has a distinct pattern of protein cleavage and the resulting peptides are more efficiently transported by TAP to the inside of the ER and more efficiently loaded in to MHC class I molecules. This allows the presentation of a distinct repertoire of antigens to CD8⁺ T lymphocytes and is important for the recognition and elimination of infected or transformed cells via the Cytotoxic T Lymphocyte (CTL) induced cell death in which perforin and granzyme play effector roles.

Previously (23), we observed in a basal breast cancer network a module formed by PSMB9, TAP1, and UBE26 genes among others. Said module was the only one with a significant enrichment in processes related to apoptosis. It is worth mentioning that in this case (undifferentiated breast cancer), we are able to observe a similar pattern with a network obtained from the unique connectivity pattern of inflammation-related genes and their neighbor genes. This result may suggest that the inflammatory process could be involved in the apoptotic response more closely than it has been reported (53).

3.6. HLA-DRA Submodule Is Related to MHC Class II Antigen Presentation

This submodule contains genes such as HLA-DRA, HLA-DRB, HLA-DQA, HLA-DQB, and HLA-DMA, which are known members of the MHC class II, and whose products are crucial in the presentation of antigens acquired in an extracellular fashion (Figures 3A, 4). MHC class II genes are expressed in professional APCs and are necessary to elicit antigen-specific responses mediated by CD4⁺ T cells. Genes in the MHC class II complex are induced in monocytes after IFNγ signaling that activates transcription factor CIITA (54).

3.7. IFI44L Submodule Is Related to Antiviral Responses and Is the Most Significantly Enriched

This submodule contains genes induced by interferon signaling. The protein products of many of these genes are known for their participation in antiviral responses. For example, OAS3 catalyzes the formation of 2′-5′ oligomers from ATP leading to the activation of RNAse 2 which degrades endogenous and viral RNAs. IFIH1 and IFIT1 are are cytoplasmic sensors of viral RNAs. PARP14 which is also coexpressed in this submodule is related to downregulation of IFNγ induced cytokines and enhances the transcription of STAT6-dependent genes (Figure 4).

3.8. IGLC2 Submodule Contains Genes for Antibody Production

This submodule consists of immunoglobulin chain genes, including heavy and light chains with variable and constant regions. Among the immunoglobulin isotype defining chains, this submodule contains the IGHA1, IGHG1, and IGHM heavy chains characteristic of IgA, IgG, and IgM antibodies. IGH1A is overexpressed, which is suggestive of IgG production, although it needs to be demonstrated at the level of protein presence in the tissue. In contrast to what is shown in CD2_sm where the T cell receptor components are coregulated, genes of the components of the B cell receptor are not present except for the IGHM chain (Figure 3).

3.9. CXCL10 Submodule Takes Part in Lymphocyte Recruitment

This submodule contains genes that code cytokines CXCL10 and CXCL11 which are CXCR3 ligands involved in lymphocyte migration as chemoattractants and are induced by IFNγ signaling (49). CCL2 and CCL8, known as MCP1 and MCP2 also are chemoattractants for leukocytes. The coordinated expression of these genes suggests an active inflammatory process where leukocytes are recruited to the tumor.

3.10. Association Between Inflammation and Adaptive Immunity

In order to quantitatively evaluate the association between inflammation and adaptive response, we performed a pathway deregulation analysis of our data by implementing the ‘Pathifier” algorithm (55) to identify sets of samples significantly deregulated (this is, samples with a pathway deregulation score, PDS > 0.4) in the inflammation and adaptive response related pathways.

From these analyses we found that out of our 641 tumor samples, a total of 395 (61.6%) were significantly deregulated in both, inflammation and adaptive response, whereas 210 (32.7%) were strongly deregulated in adaptive immunity but not in inflammation. Also 22 (3.4%) tumor samples resulted in deregulated inflammation without significant adaptive immunity changes and just 14 tumors (2.18%) were not significantly affected in adaptive immunity nor in inflammation.

It is noticeable that most tumors in our study are significantly affected in their adaptive immune responses (94.3%). This group includes most of the inflammation affected tumors (395 out of 417, or 94.7%). In contrast there are very few tumors without significant adaptive immune deregulation (36 out of 641 or 5.6% overall). Also scarce resulted tumors with notorious inflammation but with no significant adaptive immunity changes (22 out of 417 or 5.3%).

In order to further validate the aforementioned results, we have performed a similar analysis on an independent dataset including RNASeq data from the TCGA collaboration in breast cancer. The proportions are as follows. Out of 993 tumor samples, 973 samples (97.99% were strongly deregulated in both inflammation and adaptive response. Only 12 samples resulted significantly deregulated in inflammation with no noticeable deregulation in adaptive immunity (1.21%) and just 8 samples with affected adaptive immunity and no significant inflammation (0.81%). No sample among the 993 TCGA tumors was found without significant deregulation of both processes. Indeed, these figures support our claims of an important role of adaptive immunity in the breast cancer phenotypes (Supplementary Files 3, 4).

3.11. Key Genes in the Association Between Inflammation and Adaptive Immunity

Inflammation is a highly complex process, intrinsically including immunosuppressive feed-back signaling to activated T cells. For instance, expression of IL-10 and PD-L1 genes by activated myeloid cells is supposed to inhibit adaptive T cell response, whereas that of Il-12, for instance, should improve it. Similarly, adaptive immune response comprises a variety of immunosuppressive events, eventually elicited by T cell subsets expressing specific markers, e.g. FOXP3 for regulatory T cells. These events may be strongly associated to other relevant molecules in the inflammatory context, such as interferon gamma or Granzyme B.

To answer the previous questions, we have performed a generalized linear model (a form of logistic multivariate regression) for the association (membership) of the four groups of samples with the following results:

• FOXP3 expression is significantly associated (p-value = 0.0133, postive association) to adaptive immune response in breast tumors.

• FOXP3 and GZMB expression are significantly associated (p-value = 2.57 E-5 for FOXP3 and 8.14E-11 for GZMB, both postive associations) to ‘adaptive immune response plus inflammation” in breast tumors.

• The group with ‘no adaptive immune response nor inflammation” showed a discrete (although significant) negative association with FOXP3 expresion (p-value = 0.027).

• Interestingly, in all tumors FOXP3 expression is highly significant, positively associated to GZMB expression levels (p-value = 2.0E-16)

Adaptive immune response strongly depends on the B cell response. The association of the expression of genes related to B cell activation/Ig production and those involved in inflammation and T cell adaptive response is matter of interest. Being both adaptive immunity and inflammation extremely complex processes, comprising the concerted expression of hundreds of genes, we decided to evaluate that relationship by performing unsupervised learning (clustering) in the genesets associated to these processes in the four defined groups of tumor samples previously mentioned. The results are presented in the form of gene signatures displayed as clustered expression heatmaps (Supplementary File 5).

4. Discussion

4.1. Co-regulated Processes Suggest a Type I Inflammatory Microenvironment

An inflammatory reaction is elicited in response of alterations in the tissue, recognized in the form of danger signals. An initial production of signaling molecules is followed by the recruitment of immune cells which in turn integrate signals and produce new molecules to adjust the response. This is indeed what our results suggest (see Figure 5). Our findings are also consistent with the presence of a type I immune response in breast cancer. Type I responses are characterized by the production of IFNγ and the recruitment of cytotoxic T lymphocytes, which includes CD8⁺ T lymphocytes (56). This type of response is associated to virus or intracellular parasites, as well as transformed cancerous cells and culminates in the elimination of targeted cells.

FIGURE 5

Figure 5. Modules in the network are enriched in processes of complementary biological functions. CCR5_sm contains genes involved in antigen acquisition via receptor mediated phagocytosis, endosome components, as well as a number of chemokines that mediate immune cell recruitment. PSMB9_sm contains genes of components of the immunoproteasome and MHC-I molecules involved in the presentation of internally produced antigens. The immunoproteasome is induced by IFNγ signaling and is associated with a response to intracellular infections and transformation in tumor cells. HLA−DRA_sm contains MHC-II genes involved in antigen presentation to CD4⁺ T lymphocytes. Other modules like IFI44L_sm also contain interferon-induced genes annotated with antiviral functions. CD2_sm has genes involved in T lymphocyte function, including T cell receptor components and co-receptor molecule CD8A as well as other genes related to T cell cytotoxic function like perforin and granzymes. Inside tumors, antitumor mechanisms coexist with tolerance mechanisms which impair tumor immune destruction. Our results suggest an active process of tolerance where cytotoxic CD8⁺ lymphocytes turn to anergic/memory cells that no longer fight the tumor. This is supported by the coregulation of CD8A and IL7R genes in CD2_sm but not of CD8B. Other modules like IGLC2_sm contain genes of immunoglobulin chains, including constant and variable chains that may be indicative of the presence of B lymphocyte infiltrates.

Our dataset included samples from patients with developing tumors that received no neoadjuvant treatment. In this context, where type I inflammatory response is associated to good prognosis and tumor elimination, it seems paradoxical that a favorable immune response has been related to tumor growth.

4.2. Inflammation and the Shaping of the Microenvironment

The cellular composition of breast cancer tissue includes cancerous cells as well as immune cell infiltrates. Among the latter, our subnetwork suggests the presence of phagocytic and antigen-presenting cells (i. e. by the presence of MHC class II in HLA−DRA_sm and Fc receptors in CCR5_sm) and T CD8 and CD4 cells (TCR components CD3D, CD3E, CD3G, TRAC, TRB-C1, TRB-C2, and CD8α correceptor in CD2_sm, as well as CD4 and CD28 correceptors in CXCR5_sm). These are indicative of antigen-specific immune responses.

Naïve T lymphocytes that have escaped deletion by central tolerance mechanisms do not readily migrate to tissues. Instead, they first circulate and home to lymph nodes where they are exposed to antigens previously captured and presented by professional antigen presenting cells (57). If an antigen is recognized by the TCR the cell goes through changes in gene expression, proliferates and becomes an activated lymphocyte (58) (Figure 5). This process is called priming. Primed cells are able to enter inflammed tissues and exert their effector functions when they encounter their specific antigen. In some circumstances effector functions can be minimized or suppressed in the presence of antigen, what is known as peripheral tolerance.

Tolerance to self-antigens is necessary to avoid autoimmune tissue destruction, while tolerance to antigens derived from food and microbiome components keeps inflammation at bay. However, in the context of infections and cancer, tolerance may be part of the contributing factors of the pathology of chronic disease.

The effector function of CD8⁺ T cells is to kill cells infected with intracellular parasites (virus or bacteria such as mycoplasma) or transformed cells. This is achieved by means of membrane disruption and activation of the apoptosis extrinsic pathway. The latter is done by directed vesicle secretion of perforin and granzymes (PRF1, GZMA, GZMB, GZMH, GZMK in CD2_sm) over the target cell. In order to achieve this, an activated CD8⁺ cell must recognize its antigen presented in MHC class I molecules expressed in the surface of the target cell (59). Given that CD8⁺ cells respond to antigens produced an processed inside presenting cells, priming of naïve CD8⁺ by APCs must be made by cross presentation of externally acquired antigens (57).

In a CD8⁺ antitumoral reponse, antigens can be acquired in the tumor by phagocytic cells, then transported to secondary lymphoid organs and cross presented to prime CD8⁺ T cells after which activated CD8⁺ cells migrate to inflammation sites and ideally kill cells that present the antigen. The activation of CD8⁺ cells results in a pattern of clonal expansion, followed by clonal reduction and survival of memory cells, all this in consonance with the apperance, presentation and eventual clearance of the antigen. However, if antigen is persistently presented and it is recognized with high affinity, clonal reduction and effector function are altered (60). Chronic exposure to antigen is associated to downregulation of TCR and CD8 signaling components, and the expression of anergy-specific genes (i.e. Itch, Cbl-b and GRAIL) (60).

Redmod and Sherman (60) proposed a model that explains how chronic antigen exposure can lead to persistence of anergic T CD8⁺ cells. In this model, strong signaling through the TCR eliminates many of the activated cells either throug IL7Rα receptor downregulation or a decrease in antiapoptotic molecules. Meanwhile, strong TCR signaling in the surviving cells produces an increase in free calcium levels triggering anergy and a reduction in the Ras-ERK proapoptotic signaling.

In our network CD8A and IL7R genes are coregulated in the CD2_sm. IL7Rα and CD8αα are associated to the survival and development of memory CD8⁺ T cells (61, 62). This is particularly intriguing, given that our expression matrix contains both CD8A and CD8B genes for the CD8α and CD8β subunits of the CD8 correceptor, but only CD8A is included in the network. This may point out to purely CD8αα signaling, which is hypothesized to exert a weaker TCR signaling response. Anergic CD8 T cells need to be exposed to antigen to maintain anergy, and when antigen is removed, cells recover its response capacity.

For the construction of this network, interactions and genes were obtained from a background that included the whole genome (gene set defined by the HGU133 plus2 platform) and the strength and ranking of the statistical dependencies between gene pairs was not known a priori. Only 76 out of a total of 942 genes in the network are annotated as part of the inflammatory response process, which represents 8% of the network gene count. Meanwhile 694 genes in the network are annotated as part of the immune response process, that is 74% of the genes in the network.

The whole network was constructed taking into account all inflammation-associated genes, the modular and submodular structure of this network groups specific functional processes for particular modules (Figure 4). This could reflect the compartmentalization of immune response, activating sets of genes depending on the type of response that should be triggered. In these terms, general modules could be active as a first response under cellular damage or an external influence. For example, CCR5 and SPARC modules are associated to two different general cell processes: immunity (CCR5) and extracellular remodeling (SPARC). Our interpretation here is that modules partially reflect cell types present in the tumor microenvironment and that are sampled as part of the biopsy. This is specially suggested in the CD2_sm where a number of genes characteristic of CD8⁺ T lymphocytes are grouped together.

Immune response is a keystone in cancer behavior. The fact that several immune-related processes are divided into coexpressed compartments could optimize the immune response depending on the specific type of stimulus. A clear example of the aforementioned is the distinction between MHC class I and II in two separated and barely communicated submodules: PSMB9_sm is associated to MHC class I and HLA−DRA_sm takes account for MHC class II. Along with this, IFI44L_sm is associated to viral response (Figure 4). This effect has been previously observed by our group (23, 63, 64). Again, the topology of modules in the network is strongly associated to specific and separated functions.

Concomitantly, the interconnection between submodules inside the modules could be associated to fluxes of information between processes or cooperativity between modules. This is the case of CCR5_sm and CD2_sm. In Figures 3A,B it is possible to observe that the strongest link between submodules occurs in this pair. Furthermore, the number of common enriched processes between modules is the highest. Processes related to T-cell activation and adaptive immune response are shared between these two submodules, indicating possible cooperation mechanisms under the need to stimulate T-cell activity.

Interestingly enough, ZFPM1 module, despite the fact that it contains 173 genes does not result statistically enriched for any biological process at our chosen significance threshold. This results remarkable since the belonging genes are not differentially expressed (white genes clustered in left side of Figure 2). These genes may not be part of the global response that cell performs under an external stimulus.

The results obtained from our network of inflammation-associated genes suggest that in primary, untreated breast cancer tumors, the balance between antitumor response and immune tolerance involving CD8⁺ T cells is tipped in favor of the tumor. A possible mechanism being the induction of tolerance and anergization of these cells by persistent antigen exposure. This hypothesis is supported by the presence of modules in the network with genes that enrich processes and pathways related to antigen acquisition, processing and presentation, as well as genes characteristic of CD8⁺ T effector cells.

5. Summary

In breast cancer, several kind of efforts must be achieved in order to understand the complex interplay between factors shaping cancerous phenotypes. In this sense, the modular analysis of the network structure observed here, provides us with an alternative framework of the aforementioned interplay between specific modules that take account of particular biological processes that are activated in response to external stimuli. Furthermore, the observed results highlight the compartmentalized structure of co-regulated genes that may act and behave together in a well-defined time and space.

The global genetic regulatory program observed in the form of a gene network may provide insights to understand the complex mechanisms behind relevant processes for cancer development such as inflammation and the adaptive immune response. To disentangle the association between these two events by using an information theoretical approach have allowed us to analyze at an extremely detailed level such regulatory programs. Said efforts will surely increase our knowledge on the complexity of this dismal disease and provide important tools for a personalized medicine (65) aimed to produce better diagnostic, prognostic and therapeutics for breast cancer.

Author Contributions

TV-C implemented methods, performed calculations, analyzed results, drafted the manuscript. SA-C performed calculations, and contributed to the analysis of results. JE-E contributed to the analysis of results, collaborated in the discussion, reviewed the manuscript. EH-L conceived the project, contributed to the analysis of results, collaborated in the discussion, reviewed the manuscript. All authors read and approved the final version of the manuscript.

Conflict of Interest Statement

The 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.

Acknowledgments

This work was supported by CONACYT (grant no. 285544/2016 Ciencia Básica, and grant no. 2115 Fronteras de la Ciencia), as well as by federal funding from the National Institute of Genomic Medicine (Mexico). Additional support has been granted by the National Laboratory of Complexity Sciences (grant no. 232647/2014 CONACYT). TV-C is a Ph.D. student from Programa de Maestría y Doctorado en Ciencias Bioquímicas, Universidad Nacional Autónoma de México (UNAM) and received a fellowship from CONACYT (grant no 742069/59691). EH-L acknowledges additional support from the 2016 Marcos Moshinsky Fellowship in the Physical Sciences. Authors want to thak José María Zamora-Fuentes for his technical support and comments on the final version of this manuscript.

Supplementary Material

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

Supplementary File S1. A text file containing the edge list or connectivity map for the gene regulatory network (a Cytoscape SIF file).

Supplementary File S2. A text file with the definition of the modules and submodules for the largest component of the network (an Infomap MAP file).

Supplementary File S3. A multi-page PDF file with heatmaps presenting gene signatures related to the expression of genes related to specific pathways in the adaptive immunity and inflammation groups of samples in our original dataset.

Supplementary File S4. A multi-page PDF file with heatmaps presenting gene signatures related to the expression of genes related to specific pathways in the adaptive immunity and inflammation groups of samples in the validation (TCGA) dataset.

Supplementary File S5. A multi-page PDF containing gene expression boxplots for relevant genes in the different groups of inflammatory response.

References

1. Ferlay J, Soerjomataram I, Ervik M, Dikshit R, Eser S, Mathers C, et al. GLOBOCAN 2012 v1. 0, Cancer Incidence and Mortality Worldwide: IARC CancerBase No. 11 [Internet]. Lyon: International Agency for Research on Cancer [accessed 17 October, 2013].

2. Cancer Genome Atlas Network. Comprehensive molecular portraits of human breast tumours. Nature (2012) 490:61–70. doi: 10.1038/nature11412

CrossRef Full Text | Google Scholar

3. Hanahan D, Weinberg RA. The hallmarks of cancer. Cell (2000) 100:57–70. doi: 10.1016/S0092-8674(00)81683-9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell (2011) 144:646–74. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

5. DeNardo DG, Coussens LM. Balancing immune response: crosstalk between adaptive and innate immune cells during breast cancer progression. Breast Cancer Res. (2007) 9:212. doi: 10.1186/bcr1746

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Balkwill F, Mantovani A. Inflammation and cancer: back to Virchow? lancet (2001) 357:539–45. doi: 10.1016/S0140-6736(00)04046-0

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Cann SH, Van Netten J, Van Netten C. Dr William Coley and tumour regression: a place in history or in the future. Postgrad Med J. (2003) 79:672–80.

Google Scholar

8. Guarino V, Castellone MD, Avilla E, Melillo RM. Thyroid cancer and inflammation. Mol Cell Endocrinol. (2010) 321:94–102. doi: 10.1016/j.mce.2009.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Lin EY, Nguyen AV, Russell RG, Pollard JW. Colony-stimulating factor 1 promotes progression of mammary tumors to malignancy. J Exp Med. (2001) 193:727–40. doi: 10.1084/jem.193.6.727

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Terzić J, Grivennikov S, Karin E, Karin M. Inflammation and colon cancer. Gastroenterology (2010) 138:2101–14. doi: 10.1053/j.gastro.2010.01.058

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ernst P. The role of inflammation in the pathogenesis of gastric cancer. Aliment Pharmacol Therapeut. (1999) 13:13–8. doi: 10.1046/j.1365-2036.1999.00003.x

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Souza RF, Shewmake K, Beer DG, Cryer B, Spechler SJ. Selective inhibition of cyclooxygenase-2 suppresses growth and induces apoptosis in human esophageal adenocarcinoma cells. Cancer Res. (2000) 60:5767–72.

PubMed Abstract | Google Scholar

13. Ruffell B, Au A, Rugo HS, Esserman LJ, Hwang ES, Coussens LM. Leukocyte composition of human breast cancer. Proc Natl Acad Sci USA (2012) 109:2796–801. doi: 10.1073/pnas.1104303108

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Dushyanthen S, Beavis PA, Savas P, Teo ZL, Zhou C, Mansour M, et al. Relevance of tumor-infiltrating lymphocytes in breast cancer. BMC Med. (2015) 13:202. doi: 10.1186/s12916-015-0431-3

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Newman ME. The structure and function of complex networks. SIAM Rev. (2003) 45:167–256. doi: 10.1137/S003614450342480

CrossRef Full Text | Google Scholar

16. Albert R. Scale-free networks in cell biology. J Cell Sci. (2005) 118:4947–57. doi: 10.1242/jcs.02714

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Barabasi AL, Oltvai ZN. Network biology: understanding the cell's functional organization. Nat Rev Genet. (2004) 5:101–13. doi: 10.1038/nrg1272

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Newman ME, Girvan M. Finding and evaluating community structure in networks. Phys Rev E (2004) 69:026113. doi: 10.1103/PhysRevE.69.026113

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Rosvall M, Bergstrom CT. Maps of random walks on complex networks reveal community structure. Proc Natl Acad Sci USA (2008) 105:1118–23. doi: 10.1073/pnas.0706851105

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Alcalá-Corona SA, Velázquez-Caldelas TE, Espinal-Enríquez J, Hernández-Lemus E. Community structure reveals biologically functional modules in MEF2C transcriptional regulatory network. Front Physiol. (2016) 7:184. doi: 10.3389/fphys.2016.00184

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Newman ME. Modularity and community structure in networks. Proc Natl Acad Sci USA (2006) 103:8577–82. doi: 10.1073/pnas.0601602103

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Espinal-Enrıquez J, Siqueiros-Garcıa JM, Garcıa-Herrera R, Alcalá SA. A Literature-Based Approach to a Narco-Network. Barcelona: Springer (2014).

Google Scholar

23. Alcalá-Corona SA, de Anda-Jáuregui G, Espinal-Enríquez J, Hernández-Lemus E. Network modularity in breast cancer molecular subtypes. Front Physiol. (2017) 8:915. doi: 10.3389/fphys.2017.00915

PubMed Abstract | CrossRef Full Text | Google Scholar

24. de Anda-Jáuregui G, Velázquez-Caldelas TE, Espinal-Enríquez J, Hernández-Lemus E. Transcriptional network architecture of breast cancer molecular subtypes. Front Physiol. (2016) 7:568. doi: 10.3389/fphys.2016.00568

CrossRef Full Text | Google Scholar

25. Edgar R, Domrachev M, Lash AE. Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. (2002) 30:207–10. doi: 10.1093/nar/30.1.207

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Clarke C, Madden SF, Doolan P, Aherne ST, Joyce H, O'Driscoll L, et al. Correlating transcriptional networks to breast cancer survival: a large-scale coexpression analysis. Carcinogenesis (2013) 34:2300–8. doi: 10.1093/carcin/bgt208

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Lisowska KM, Dudaladava V, Jarzab M, Huzarski T, Chmielik E, Stobiecka E, et al. BRCA1-related gene signature in breast cancer: the role of ER status and molecular type. Front Biosci. (Elite Ed). (2011) 3:125–36. doi: 10.2741/e227

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tan TZ, Miow QH, Miki Y, Noda T, Mori S, Huang RYJ, et al. Epithelial-mesenchymal transition spectrum quantification and its efficacy in deciphering survival and drug responses of cancer patients. EMBO Mol Med. (2014) 6:1279–93. doi: 10.15252/emmm.201404208

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Pedraza V, Gomez-Capilla JA, Escaramis G, Gomez C, Torné P, Rivera JM, et al. Gene expression signatures in breast cancer distinguish phenotype characteristics, histologic subtypes, and tumor invasiveness. Cancer (2010) 116:486–96. doi: 10.1002/cncr.24805

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Smyth GK, Ritchie M, Thorne N, Wettenhall J. LIMMA: linear models for microarray data. In: Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Statistics for Biology and Health (2005).

Google Scholar

31. Bansal M, Belcastro V, Ambesi-Impiombato A, Di Bernardo D. How to infer gene networks from expression profiles. Mol Syst Biol. (2007) 3:78. doi: 10.1038/msb4100158

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hernández-Lemus E, Rangel-Escareño C. The role of information theory in gene regulatory network inference. Inform Theory (2011) 109–44.

33. Hernández-Lemus E, Velázquez-Fernández D, Estrada-Gil JK, Silva-Zolezzi I, Herrera-Hernández MF, Jiménez-Sánchez G. Information theoretical methods to deconvolute genetic regulatory networks applied to thyroid neoplasms. Physica A (2009) 388:5057–69. doi: 10.1016/j.physa.2009.08.013

CrossRef Full Text | Google Scholar

34. Margolin AA, Nemenman I, Basso K, Wiggins C, Stolovitzky G, Favera RD, et al. aracne: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics (2006) 7(Suppl 1):S7. doi: 10.1186/1471-2105-7-S1-S7

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Rosvall M, Bergstrom CT. Multilevel compression of random walks on networks reveals hierarchical organization in large integrated systems. PLoS ONE (2011) 6:e18209. doi: 10.1371/journal.pone.0018209

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Alcalá-Corona SA, Espinal-Enríquez J, De Anda Jáuregui G, Hernandez-Lemus E. The hierarchical modular structure of HER2+ breast cancer network. Front Physiol. (2018) 9:1423. doi: 10.3389/fphys.2018.01423

PubMed Abstract | CrossRef Full Text | Google Scholar

37. García-Campos MA, Espinal-Enríquez J, Hernández-Lemus E. Pathway analysis: state of the art. Front Physiol. (2015) 6:383. doi: 10.3389/fphys.2015.00383

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang J, Duncan D, Shi Z, Zhang B. WEB-based GEne SeT anaLysis toolkit (WebGestalt): update 2013. Nucleic Acids Res. (2013) 41:W77–83. doi: 10.1093/nar/gkt439

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Wang X, Terfve C, Rose JC, Markowetz F. HTSanalyzeR: an R/Bioconductor package for integrated network analysis of high-throughput screens. Bioinformatics (2011) 27:879–80. doi: 10.1093/bioinformatics/btr028

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ogden CA, Hoffmann PR, Bratton D, Ghebrehiwet B, Fadok VA, Henson PM, et al. C1q and mannose binding lectin engagement of cell surface calreticulin and CD91 initiates macropinocytosis and uptake of apoptotic cells. J Exp Med. (2001) 194:781–96. doi: 10.1084/jem.194.6.781

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ogden CA, Kowalewski R, Peng Y, Montenegro V, Elkon KB. IGM is required for efficient complement mediated phagocytosis of apoptotic cells in vivo. Autoimmunity (2005) 38:259–64. doi: 10.1080/08916930500124452

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Chen G, Tan CS, Teh BK, Lu J. Molecular mechanisms for synchronized transcription of three complement C1q subunit genes in dendritic cells and macrophages. J Biol Chem. (2011) 286:34941–50. doi: 10.1074/jbc.M111.286427

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Korb LC, Ahearn JM. C1q binds directly and specifically to surface blebs of apoptotic human keratinocytes: complement deficiency and systemic lupus erythematosus revisited. J Immunol. (1997) 158:4525–8.

PubMed Abstract | Google Scholar

45. Mantegazza AR, Savina A, Vermeulen M, Pérez L, Geffner J, Hermine O, et al. NADPH oxidase controls phagosomal pH and antigen cross-presentation in human dendritic cells. Blood (2008) 112:4712–22. doi: 10.1182/blood-2008-01-134791

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Delamarre L, Pack M, Chang H, Mellman I, Trombetta ES. Differential lysosomal proteolysis in antigen-presenting cells determines antigen fate. Science (2005) 307:1630–4. doi: 10.1126/science.1108003

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Farber JM. A macrophage mRNA selectively induced by gamma-interferon encodes a member of the platelet factor 4 family of cytokines. Proc Natl Acad Sci USA (1990) 87:5238–42. doi: 10.1073/pnas.87.14.5238

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Luster AD, Unkeless JC, Ravetch JV. γ-Interferon transcriptionally regulates an early-response gene containing homology to platelet proteins. Nature (1985) 315:672. doi: 10.1038/315672a0

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Weng Y, Siciliano SJ, Waldburger KE, Sirotina-Meisher A, Staruch MJ, Daugherty BL, et al. Binding and functional properties of recombinant and endogenous CXCR3 chemokine receptors. J Biol Chem. (1998) 273:18288–91. doi: 10.1074/jbc.273.29.18288

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Pearce EL, Mullen AC, Martins GA, Krawczyk CM, Hutchins AS, Zediak VP, et al. Control of effector CD8+ T cell function by the transcription factor Eomesodermin. Science (2003) 302:1041–3. doi: 10.1126/science.1090148

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Fehres CM, Unger WW, Garcia-Vallejo JJ, van Kooyk Y. Understanding the biology of antigen cross-presentation for the design of vaccines against cancer. Front Immunol. (2014) 5:149. doi: 10.3389/fimmu.2014.00149

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Strehl B, Seifert U, Krüger E, Heink S, Kuckelkorn U, Kloetzel PM. Interferon-γ, the functional plasticity of the ubiquitin–proteasome system, and MHC class I antigen processing. Immunol Rev. (2005) 207:19–30. doi: 10.1111/j.0105-2896.2005.00308.x

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Alotaibi MR, Hassan ZK, Al-Rejaie SS, Alshammari MA, Almutairi MM, Alhoshani AR, et al. Characterization of apoptosis in a breast cancer cell line after IL-10 silencing. Asian Pacific J Cancer Prevent (2018) 19:777. doi: 10.22034/APJCP.2018.19.3.777

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Schroder K, Hertzog PJ, Ravasi T, Hume DA. Interferon-γ: an overview of signals, mechanisms and functions. J Leukocyte Biol. (2004) 75:163–89. doi: 10.1189/jlb.0603252

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Drier Y, Sheffer M, Domany E. Pathway-based personalized analysis of cancer. Proc Natl Acad Sci USA (2013) 110:6388–93. doi: 10.1073/pnas.1219651110

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Eberl G. Immunity by equilibrium. Nat Rev Immunol. (2016) 16:524. doi: 10.1038/nri.2016.75

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Kurts C. Cross-presentation: inducing CD8 T cell immunity and tolerance. J Mol Med. (2000) 78:326–32. doi: 10.1007/s001090000108

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Obst R. The timing of T cell priming and cycling. Front Immunol. (2015) 6:563. doi: 10.3389/fimmu.2015.00563

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Lieberman J. Cell death and immunity: the ABCs of granule-mediated cytotoxicity: new weapons in the arsenal. Nat Rev Immunol. (2003) 3:361. doi: 10.1038/nri1083

CrossRef Full Text | Google Scholar

60. Redmond WL, Sherman LA. Peripheral tolerance of CD8 T lymphocytes. Immunity (2005) 22:275–84. doi: 10.1016/j.immuni.2005.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Kaech SM, Tan JT, Wherry EJ, Konieczny BT, Surh CD, Ahmed R. Selective expression of the interleukin 7 receptor identifies effector CD8 T cells that give rise to long-lived memory cells. Nat Immunol. (2003) 4:1191. doi: 10.1038/ni1009

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Madakamutil LT, Christen U, Lena CJ, Wang-Zhu Y, Attinger A, Sundarrajan M, et al. CD8αα-mediated survival and differentiation of CD8 memory T cell precursors. Science (2004) 304:590–3. doi: 10.1126/science.1092316

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Alcalá-Corona SA, de Anda-Jáuregui G, Espinal-Enriquez J, Tovar H, Hernández-Lemus E. Network modularity and hierarchical structure in breast cancer molecular subtypes. In: International Conference on Complex Systems. Cambridge, MA: Springer (2018). p. 352–8.

Google Scholar

64. Espinal-Enriquez J, Fresno C, Anda-Jáuregui G, Hernandez-Lemus E. RNA-Seq based genome-wide analysis reveals loss of inter-chromosomal regulation in breast cancer. Sci Rep. (2017) 7:1760. doi: 10.1038/s41598-017-01314-1

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Mejía-Pedroza RA, Espinal-Enríquez J, Hernández-Lemus E. Pathway-based drug repositioning for breast cancer molecular subtypes. Front Pharmacol. (2018) 9:905. doi: 10.3389/fphar.2018.00905

PubMed Abstract | CrossRef Full Text | Google Scholar