Impact of chemotactic factors and receptors on the cancer immune infiltrate: a bioinformatics study revealing homogeneity and heterogeneity among patient cohorts

Datasets and expression variability

We first assembled a series of public microarray data into one single database. Five different datasets were included for several major cancer types, in particular breast carcinoma (TCGA consortium),^31–34 colorectal carcinoma (http://www.intgen.org or TCGA consortium),^35,36 non-small cell lung carcinoma (TCGA consortium)^37–40 and melanoma,^41–45 to study the correlation between the expression level of metagenes indicating the presence of immune cell types and a variety of chemotactic factors and receptors. An extra dataset concerning breast carcinoma⁴⁶ was used for studying expression variability and treatment response. We also used two datasets for ovarian (TCGA consortium)⁴⁷ and head and neck carcinoma^48,49 each. The total number of tumors analyzed in this study amounts to 5953 (Table 1). We first examined the expression variability of each of the metagenes relevant to tumor-infiltrating immune cell types (determined by using the MCP counter method)^27,30 and mRNAs encoding annexins (ANXA), chemokines (CCL, CXCL, XCL), chemokine receptors (CCR, CXCR, XCR, CCRL), formyl peptide receptors (FPR), purinergic receptors (ADOR, P2RX, P2RY), compared to housekeeping genes (ACTB, GADPH, TUB), as this is typically done when protein expression is measured by immunoblot analysis. However the variability in the expression of housekeeping genes was not expected to be smaller than that of the genes of interest (including immune metagenes) due to the normalization methods of microarrays. Nevertheless, a systematic low variance of genes of interest (including immune metagenes) compared to housekeeping genes, might produce correlation coefficients (computed below) that are overestimated. Therefore, genes that have a systematic low variance should be removed from the analysis. We analyzed the relative variance of each gene of interest, expressed in a log₁₀ scale for all datasets (Figure 1). In addition, variance difference tests were performed, between housekeeping genes and genes of interest (Supplementary Figure 1). As a general tendency, the selected genes and metagenes corresponding to distinct immune types exhibited a normalized variance around 1 (i.e. 0 on a log₁₀ scale), with no major difference between immune-relevant and housekeeping genes. Nevertheless, some genes have a smaller variance compared to housekeeping genes, especially in the two breast cancer datasets (METABRIC and TCGA). As a result, we decided to include these close-to-invariant genes in the analyses, because the possible overestimation of correlation coefficients is not systematic.

General correlation patterns across cancer types

We computed the Spearman correlation coefficients between immune cell type-relevant metagenes and our selection of chemotaxis-relevant genes. We used Fisher’s method (or Fisher’s combined probability test, for multiple testing) to identify reproducible correlations for five cancer types, namely breast cancer (Figure 2), colorectal cancer (Figure 3a), non-small cell lung cancer (Figure 3b) and melanoma (Figure 3c), while performing a meta-analysis of five different data sets for each cancer (Table 1). This procedure was not applicable to ovarian and head & neck cancers because it was not possible to collect 2 × 5 sufficiently large datasets. Correlations that were reproducible across the five data sets were color-coded to indicate positive (red) or negative (green) associations between the expression levels of individual chemotactic factors/receptors and leukocyte subset-specific metagenes, while the absence of reproducible correlations were indicated by white boxes. As to be expected, we detected a mostly positive correlation between leukocyte subsets and chemotactic factors/receptors in all four cancer types (Figures 2 and 3). The reproducibility of correlations was less pronounced in melanoma (Figure 3C), probably because the datasets available for this cancer type were comparatively small (Table 1).

When all these figures were grouped together, (Supplementary Figure 2), one could observe that the variability seemed larger among cancer types than among distinct immune cell types. For this reason, we performed separate analyses of distinct cancer types in different (sub-)figures and then grouped immune cell types. We applied Fisher’s method to identify correlations that were reproducible across the aforementioned five cancer types (Figure 4). Unsupervised hierarchical clustering of these results identified two types of genes coding for chemotactic factors/receptors, namely (i) genes that exhibit reproducible positive correlations with most of the leukocyte subtypes (in the lower part of Figure 4) and (ii) genes that have almost no reproducible correlations with immune cell types (in the upper part of Figure 4). The interpretation of this separation in two groups is not obvious. The second group of chemotactic factors/receptors could be interpreted as a cancer type-specific immune system activation. To visualize this possibility, we removed the reproducible correlations across cancer types (Figure 4) from each organ-specific correlation pattern (Figures 2 and 3) with the scope of identifying cancer type-specific correlation patterns (Supplementary Figure 3–6). Again, this separation between fully reproducible and cancer specific-reproducible correlations supported the idea that part of the chemotaxis of different leukocyte subtypes may occur in an organ-specific fashion. We tested the robustness of this two-group separation of chemotactic factors/receptors, by changing the p-value threshold used for reproducibility tests from 0.05 to 0.01. Changing the threshold of the p-value did not affect the final result (Supplementary Fig. 7B) as in Figure 4.

Next, we investigated the reproducible correlation pattern identified in Figure 4 on two other cancer types, namely, ovarian cancer and head and neck carcinoma. Because the number of large datasets for these two cancers types is rather limited (two datasets for each of them, Table 1), it was not possible to determine reproducible correlation patterns within each of these two cancer types. For each dataset, we computed the correlation matrix between immune cell types and the selection of genes. Then, we removed the reproducible correlations identified in Figure 4. The remaining “specific” correlations are illustrated in Supplementary Fig. 8–11. The number of such specific correlations was rather small, exhibiting no obvious pattern and hence resembled those obtained for breast, colorectal, non-small cell lung cancer and melanoma (Supplementary Figure 3–6). Therefore, ovarian cancer and head and neck carcinoma exhibited a similar global correlation pattern between chemotactic factors/receptors and the density of the immune infiltrate, supporting the universal validity of the common correlation pattern (Figure 4) across distinct malignancies.

One question that can be raised concerns the cancer specificity of these reproducible correlations (fully reproducible or organ-specific reproducible) compared to normal tissue. Although a full analysis, comparing normal and cancer tissue is beyond the scope of this work, we produced a similar Figure than Figure 2, using normal breast samples available from the METABRIC dataset (Supplementary Fig. 12). The global pattern is not strongly different from Figure 2. Therefore, it seems that the mechanism involving chemotactic factors/receptors in immune infiltrate in cancer in not fundamentally different from a normal immune activation.

The reproducible correlation pattern across cancers (Figure 4) has been established without any knowledge about the abundance of the chemotactic factors/receptors in different immune cell types or in distinct cancer types. This information is not directly accessible from microarray data, because of normalization methods. Nevertheless, we explored the possibility that the abundance of chemotactic factors/receptors that have a reproducible correlation with immune cell types might be higher or lower than that of those failing to exhibit reproducible correlations, while investigating absolute expression in our microarray datasets and the mRNA expression data from the Human Protein Atlas.^50,51 Chemotactic factors/receptors were removed from the analyses if their expression was systemically low in the microarray data or in the spleen and lymph node (see Materials and Methods for the exact procedure). Among the 66 chemotactic factors/receptors that were classed as non-reproducible in their correlation with immune metagenes (in Figure 4), 12 exhibited a low expression. Among the remaining 39 chemotactic factors/receptors reproducibly correlating with immune metagenes (in Figure 4), only one exhibited a low expression. Hence, correlation reproducibility is anti-correlated with low expression (p-value = 0.029, exact Fisher’s test). Indeed, lowly expressed genes should not yield reproducible correlations with immune cell infiltrates. Because the choice for filtering thresholds was quite arbitrary (see Material and Methods), we applied the same procedure by a) doubling threshold on our mRNA datasets, b) doubling threshold on mRNA data from Human Protein Atlas. We still observed the same tendency; in the case a) Fishers’s p-value = 0.155 (filtering 2 genes from the reproducible list and 14 from the non-reproducible one), in case b), Fishers’s p-value = 0.028 (filtering 2 genes from the reproducible list and 14 from the non-reproducible one).

Figure 4 reveals 13 chemotactic receptors and chemokines that positively correlated with all identified cell subsets (irrespective of their stromal and leukocytic nature) contained in the tumor environment, regardless of the cancer type. This applies to CCR1,2,5; CXCR4; FPR1; ADORA2A; P2RY6,13,14; CCL18,19,21; and CXCL12. In other words, these factors appear intrinsically linked to the neoplastic process. Accordingly, CXCL12 and its receptor CXCR4 have been involved in promoting cancer proliferation, survival, invasion, metastasis, stemness and angiogenesis.^52–57 CCR1 and its ligand CCL5 have been associated with cancer cell invasion and metastasis.^58–60 Similar observations have been reported for the CCR2-CCL2, CCR5-CCL2,3,4,5, CCR7-CCL19,21 and CCR8-CCL18 axes.^61–69 The role of FPR1 in tumor progression remains controversial as it has been associated with tumor invasion in colorectal cancer but with tumor suppression in gastric cancer.^70,71 As far as P2RY6 is concerned, Placet et al. recently suggested its role in cancer cell survival.⁷² Overall, these 13 chemotactic receptors and chemokines may constitute universal targets for cancer therapy. Further consolidating the relevance of our results, additional correlations between tumor-infiltrating immune cell subtypes and chemokine/receptors across cancer types previously described in separate preclinical/clinical studies were found (Figure 4). Among these, we can mention the association between the infiltration of T cells and the requirement of the receptor ADORA2⁷³ or the role of CXCR3 and its ligands CXCL9 to 11 in the recruitment of CD8⁺ T cells, Th1 cells and NK cells.¹⁸ However, Figure 4 also revealed original interactions between immune cell subtypes. For instance, neutrophil infiltrates appeared negatively correlated with the chemokines CXCL9-11. Based on these results, the recruitment of T lymphocytes and neutrophils into the tumor microenvironment could be mutually exclusive. Indeed, a recent work demonstrating an inhibitory activity of CXCL9 peptides on neutrophil migration in murine models supporting such a hypothesis.⁷⁴

Impact of chemotactic factors and receptors on breast cancer prognosis

It is well established that breast cancer is under immunosurveillance and that neoadjuvant chemotherapy is particularly successful if the cancer is massively infiltrated by immune effectors including T lymphocytes.^75–80 Driven by the dichotomy between mRNAs coding for chemotactic factors/receptors that either (i) positively correlated with leukocyte infiltration across distinct cancer types or (ii) lacked a clear pattern of correlation (Figure 4), we investigated the impact of these two groups of genes on the survival of breast cancer patients. We subdivided the second group into two subgroups (see Supplementary Table 1): a) genes that positively correlate with leukocyte infiltration in breast cancer (but fail to do so across the entire spectrum of cancers analyzed in this study), b) remaining chemotactic factors/receptors (that fail to correlate with immune metagenes in breast cancer). We also considered the group of chemotactic factors/receptors that positively correlated with leukocyte infiltration across distinct cancer types, without the single genes that are characterized by low expression (see Supplementary Table 1).

We explored the link between chemotactic factors/receptors and response to treatment in the four out of the five Breast Carcinoma datasets for which the treatment response was annotated (as pathological complete response or simple response, see Table 2) after neo-adjuvant chemotherapy: Bonnefoi,³² Hatzis,³³ Tabchy,³⁴ and Korde.⁴⁶ We also included the METABRIC dataset⁸¹ in this analysis, because we had access to the information whether patients were alive or deceased, a status that is related to disease prognosis and treatment response. In four out of these five datasets (all except the Korde cohort, probably due to its small size), expression of chemotactic factors/receptors that positively correlated with leukocyte infiltration across distinct cancer types (as defined in Figure 4) was positively associated with treatment response (Bonnefoi, Hatzis and Tabchy datasets) or overall survival (METABRIC dataset) (Table 2). METABRIC results differed for non-reproducible chemotactic factors/receptors, probably because treatment response is measured indirectly. In contrast, no correlation was observed between the remaining chemotactic factors/receptors and prognostic features (Table 2). This is also illustrated by the Volcano Plots of different groups of chemotactic factors/receptors (Supplementary Figure 13). Because METABRIC dataset has the clinical information about survival, we attempted to build simple biomarkers (see Material and Methods) from the different groups of chemotactic factors/receptors (Supplementary Figure 14). Only the group of chemotactic factors/receptors reproducibly detected across different cancer types tended to separate the cohort into two groups based on the expression level of the biomarkers, with a difference between survival outcomes that happened to be almost significant (p = 0.081). Altogether, these data support the idea that chemotactic factors and receptors that positively and reproducibly correlate with the immune infiltrate of different cancer types play a positive role in local immunosurveillance, hence improving the probability of complete therapeutic responses in the context of neoadjuvant chemotherapy. The global function of this group of reproducible chemotactic factors/receptors is not well established, and therefore this group cannot be considered as a single metagene, like for immune cell types. Nevertheless it may be interesting to compare the differential abundance of these chemotactic factors/receptors with the differential expression of the immune cell metagenes (in Supplementary Table 2, one sided t-tests of differential expression of treatment response). Although the combined t-tests produce significant differential expression of immune cell metagenes, only one cohort (Bonnefoi) has significant p-values on its own (thresholds 0.05), which appears less significant than the tests for chemotactic factors/receptors of Table 2. It is always difficult to interpret such variation between cohorts. The main reason could reside in the treatment regimens and in the clinical definition of treatment response that were not the same for all cohorts.

Concluding remarks

The present bioinformatic study has been conducted with the explicit aim of identifying universal (rather than organ-specific) correlations between the expression level of mRNAs coding for chemotactic factors and receptors, and metagene signatures reflecting distinct leukocyte subtypes infiltrating a variety of cancer types. The specific strategy that we have chosen for this meta-analysis addresses important aspects with respect to the homogeneity and heterogeneity of patient cohorts. Undoubtedly, there are major interindividual variations in the density, composition, architecture and functional state of the tumor infiltrate.^9,11,82,83 Notwithstanding this fact, the variations in metagene expression across different cancer types were not higher than those observed for a well known ‘house-keeping’ gene (ACTB, GAPDH, TUB), suggesting that the heterogeneity of the immune infiltrate is not substantially higher than that of the expression of mRNAs coding for essential components of the cytoskeleton (ACTB) or glycolytic metabolism (GAPDH).

A second notion that emerges from this meta-analysis concerns the regulation of the immune infiltrate by chemotactic factors and receptors occurring in a universal (homogeneous) pattern across distinct cancer types, irrespective of their localization in the body, versus organ-specific (heterogeneous) patterns. Surprisingly, it appears that most of the (positive) correlations between the chemotactic factors and receptors analyzed herein and different infiltrating leukocyte subtypes are highly reproducible across different cancer types, while few of them are organ-specific. This points to the importance of common rules dictating the recruitment and permanence of immune cells into the tumor bed, irrespective of the precise cancer type that is concerned. Such universal rules appear to supersede in importance those that may influence organ-specific circuitries. Because this analysis is based on correlation, the causal interpretation is not obvious (are chemotactic factors/receptors attracted by cell types, or is it the contrary?). A better understanding of these circuitries would include a description/modeling of signaling networks at the microenvironment level, which is beyond the scope of the present work.

Most importantly, it appears that the mRNAs coding for chemotactic factors and receptors that reproducibly and positively correlate with the local presence of distinct immune cell types seem to have a positive impact on breast cancer prognosis. This contrasts with chemotactic factors and receptors that exhibit heterogeneous and organ-specific patterns of correlation with the immune infiltrate that have no such prognostic impact. On a positive note, this implies that, at least in the realm of chemotaxis, insights that have been gained in one cancer type can be extrapolated to other malignancies, irrespective of their cell of origin.

Datasets and immune infiltrate estimation

We explored public datasets of transcriptome microarrays. The complete list of these datasets is shown in Table 1. We use data that were already normalized, as provided from the different repository websites. The MCPcounter²⁷ R-package was used for estimating immune cell type activities. This method uses groups of genes to construct metagenes, whose expressions measure indirectly immune cell abundances. The identification of gene groups have been done in a strongly reproducible manner, learned from datasets where immune cell type abundances are known.

Test of variance difference

One sided F-test was used, to test if variance of genes of interest (including immune metagenes) is smaller than variance of housekeeping genes.

Patterns of reproducible correlation

For each dataset, we computed the correlation coefficients between immune cell activities and expression of individual mRNAs coding for chemotactic factors/receptors. We decided to use the Spearman’s correlation method, because the normal hypothesis necessary for using Pearson’s approach is not guaranteed for immune activities estimated by MCP counter method. For each cancer type, we used the Fisher’s method in the following way. We computed one-sided Spearman correlation test for both cases (H1: positive correlation and H1: negative correlation). For each case (H1 positive and H1 negative), we combined the p-value by using Fisher’s method (Chi-square test on sum of p-value logarithms). If both p-values were higher than 0.05, the correlation was considered as non-reproducible. In any other case, we considered that the correlation was reproducible. We then used the coefficient that has the smallest p-value and whose sign is compatible with the hypothesis H1 that had a significant combined p-value. These correlation coefficients are illustrated in Figure 2–3.

For the fully reproducible correlation coefficients listed in Figure 6, we applied the same procedure as above. Thus, we considered the correlations coefficients of Figure 2–3 and their associated p-values, as computed above (combined one-sided correlation tests by Fisher’s method). These p-values were one-sided ones, related to the sign of the associated correlation coefficient. We computed the p-values associated to one-sided correlation test of opposite sign (‘one sided p-value’ → (1-‘one sided p-value’). Therefore, we could apply Fisher’s method for both H1 hypothesis (positive and negative correlation) to a common set of cancer types. Again, if both p-values were higher than 0.05, the correlation was interpreted as non-reproducible. In the other case, we considered that correlation was reproducible; we then used the coefficient that has the smallest p-value and whose sign is compatible with the hypothesis H1 that has a significant combined p-value. These correlation coefficients are shown in Figure 4. Suppl. Fig 7 is produced by using the same analysis, replacing the p-value threshold of 0.05 (Figure 4) by 0.01 (Suppl. Fig. 7).