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

Multiple soluble factors including proteins (in particular chemokines), non-proteinaceous factors released by dead cells, as well as receptors for such factors (in particular chemokine receptors, formyl peptide receptors and purinergic receptors), influence the recruitment of distinct cell subsets into the tumor microenvironment. We performed an extensive bioinformatic analysis on tumor specimens from 5953 cancer patients to correlate the mRNA expression levels of chemotactic factors/receptors with the density of immune cell types infiltrating the malignant lesions. This meta-analysis, which included specimens from breast, colorectal, lung, ovary and head and neck carcinomas as well as melanomas, revealed that a subset of chemotactic factors/receptors exhibited a positive and reproducible correlation with several infiltrating cell types across various solid cancers, revealing a universal pattern of association. Hence, this meta-analysis distinguishes between homogeneous associations that occur across different cancer types and heterogeneous correlations, that are specific of one organ. Importantly, in four out of five breast cancer cohorts for which clinical data were available, the levels of expression of chemotactic factors/receptors that exhibited universal (rather than organ-specific) positive correlations with the immune infiltrate had a positive impact on the response to neoadjuvant chemotherapy. These results support the notion that general (rather than organ-specific) rules governing the recruitment of immune cells into the tumor bed are particularly important in determining local immunosurveillance and response to therapy.

Introduction

Although cancer has been viewed for a long time as a cell-autonomous, genetic and epigenetic disease, there is little doubt that non-transformed cells, in particular leukocytes have a major influence on the development, progression, therapeutic response and fatal relapse of malignant disease.[1-4] In essence, the relationship between cancer pathogenesis and leukocytes can be defined as ambiguous. On one hand, chronic subclinical or manifest inflammation can contribute to oncogenesis by favoring local generation of reactive oxygen species, inducing genomic instability, by increasing the turnover of parenchymal or stromal cells, increasing the probability of genetic and chromosomal aberrations, by disrupting homeostatic regulation in the affected tissue, favoring the emergence of more ‘competitive’ cells, and by inducing local immunosuppression, hence subverting immunosurveillance.[5,6] On the other hand, specific leukocyte subpopulations may contribute to the recognition of premalignant and fully transformed cells, causing their elimination (the best-case scenario) or engaging in a precarious equilibrium state between malignant cells and immune effectors. In most cases, it is only when cancer cells manage to suppress the anticancer immune response (or to hide from recognition) that the disease becomes clinically detectable.[2,4,7-11]

The attraction of leukocytes into premalignant and malignant lesions is dictated by multiple soluble factors including so-called danger-associated molecular patterns (DAMPs) that include extracellular nucleotides and nucleosides (with ATP as a prominent chemotactic factor for myeloid cells), proteins that are usually confined in intracellular compartments yet are released from stressed, dying and dead cells (such as annexin-1, ANXA1; calreticulin, CALR, F-actin, and high molecular group protein-1, HMGB1)[12-17] and proteins that are actively secreted by cancer or stromal cells, in particular chemokines (C-C motif chemokine ligand 1 to 9 and 11 to 28, CCL1-9, CCL11-28; X-C motif chemokine ligand 1 and 2, XCL1 and 2; C-X-C motif chemokine ligand 1 to 17, CXCL1-17; C-X3-C motif chemokine ligand 1, CX3CL1).[18] These factors act on a series of receptors, such as the adenosine receptor family AdoR (ADOR A1, A2A, A2B and A3), metabotropic purinergic receptors (P2RY1, 2, 4, 6, 8 and 10 to 14), ionotropic purinergic receptors (P2RX1 to 7), formylpeptide receptors (FPR1 to FPR3) and chemokine receptors (XCR1, CXCR1 to 7, CX2CR1, CCR1 to 10, CCRL1 and 2) to induce the chemotaxis of different leukocyte subtypes into the tumor.[18-21] The expression of such chemotactic factors and their receptors determines the composition of the tumor immune infiltrate and has a major impact on the therapeutic response and prognosis of cancer.[9,22-25]

Most of the knowledge about the immune infiltrate that can be retrieved from large data sets is based on microarray or RNAseq data obtained on whole tumor specimens and biopsies, requiring bioinformatic deconvolution of the data.[26] Methods for estimating the abundance of leukocyte subsets in tumor, based on gene expression data (microarray or RNAseq), have been developed.[27-30]

Here, we report a meta-analysis of the relationship between the relative abundance of distinct leukocyte subsets and the level of expression of individual chemotactic factors and their receptors in a large set of human tumors.

Results and discussion

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[46] was used for studying expression variability and treatment response. We also used two datasets for ovarian (TCGA consortium)[47] 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 log10 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 log10 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.

Table 1.

Overview on the cohorts included in this meta-analysis.

Cancer type	Cohort name	Number of samples	Characteristics of the cohort	Treatment & outcome	Reference	Platform
Melanoma	Xu	83	Primary and metastatic tumors		GSE8401	Affymetrix Human Genome U133A Array
Melanoma	Harlin	44	Metastatic tumors		GSE12627	Affymetrix Human Genome U133A Array
Melanoma	Bogunovic	44	Metastatic tumors		GSE19234	Affymetrix Human Genome U133 Plus 2.0 Array
Melanoma	RikerMel	56	Primary and metastatic tumors		GSE7553	Affymetrix Human Genome U133 Plus 2.0 Array
Melanoma	Talantov	45	Primary tumors		GSE3189	Affymetrix Human Genome U133A Array
Colon	BittColon	307	Various colon tumors		GSE2109	Affymetrix Human Genome U95 Version 2 Array
Colorectal	Smith	177	Various colorectal tumors		GSE17536	Affymetrix Human Genome U133 Plus 2.0 Array
Colon	Vilar1	155	Colon tumors		GSE26682, 1st set	Affymetrix Human Genome U133A Array
Colon	Vilar2	176	Colon tumors		GSE26682, 2nd set	Affymetrix Human Genome U133A Array
Colon	TCGA	174	Various colon tumors		TCGA consortium	Agilent 244K Custom Gene Expression G4502A-07–3
Breast	METABRIC	1781	Various breast tumors		METABRIC Study	Illumina HT-12 v3
Breast	TCGA	522	Various breast tumors		TCGA consortium	Agilent 244K Custom Gene Expression G4502A-07–1
Breast	Bonnefoi	161	Locally advance or large operable breast tumors, estrogen receptor negative	FEC or ET treatment. Pathological complete response (complete disappearance of the tumour with no more than a few scattered tumour cells) vs no pathological complete response	GSE6861	Affymetrix Human X3P Array
Breast	Hatzis	198	HER2 negative breast tumors	Taxane-anthracycline chemotherapy pre-operatively and endocrine therapy if ER-positive. Pathological complete response (no invasive or metastatic breast cancer identified) vs rapid development	GSE25065	Affymetrix Human Genome U133A Array
Breast	Tabchy	178	Various type of breast tumors before treatment	FEC or FAC neo-adjuvant chemotherapy. Pathological complete response vs residual disease (clinical or radiological progression)	GSE20271	Affymetrix Human Genome U133A Array
Breast	Korde	61	Various type of breast tumors, stage 2 or 3 breast cancer with tumor size ≥2cm at patients selection, prior to AC treatment	4 cycles of TX, 4 cycles of adriamycin, cyclophosphamide on day 1 and 21 (neoadjuvant) and AC (neo-adjuvant or adjuvant). Response vs no response (change in tumor size by clinical exam and pathological response).	GSE18728	Affymetrix Human Genome U133 Plus 2.0 Array
Lung	AdenoConsortium	462	Various type of Adenocarcinomas		Director’s Challenge Lung Study, National Cancer Institute (NHI)	Affymetrix Human Genome U133A Array
Lung	Lee	138	Adenocarcinoma and squamous cell carcinoma		GSE8894	Affymetrix Human Genome U133 Plus 2.0 Array
Lung	Okayama	226	Adenocarcinoma		GSE31210	Affymetrix Human Genome U133 Plus 2.0 Array
Lung	Raponi	130	Squamous cell carcinoma		GSE4573	Affymetrix Human Genome U133A Array
Lung	TCGA	134	squamous cell carcinoma		TCGA consortium	Affymetrix Human Genome U133A Array
Ovarian	Yoshihara	43	Ovarian carcinomas		GSE12470	Agilent-012097 Human 1A Microarray (V2)
Ovarian	TCGA	520	Ovarian carcinomas		TCGA consortium	Affymetrix Human Genome U133A Array
Head & Neck	Peng	57	Head and Neck carcinomas		GSE25099	Affymetrix Human Exon 1.0 ST Array
Head & Neck	Rickman	81	Head and Neck carcinomas		E-TABM-302	Affymetrix Human Genome U133 Plus 2.0 Array

Figure 1.

Heatmap representation of gene expression variances and immune cell type activity variances, in log10, for the different datasets. For each gene, in each dataset, gene variance was normalized by the variance all gene expressions pooled together. For each immune cell type activities, in each dataset, immune cell type activity variance is normalized by the variance all immune cell type activities pooled together. Housekeeping genes are highlighted.

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

Figure 2.

Heatmap representation of reproducible Spearman’s correlation coefficients, in breast carcinomas.

Figure 3.

Heatmap representation of reproducible Spearman’s correlation coefficients, in colorectal carcinomas (A), in non-small cell lung cancers (B) and in melanomas (C).

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.

Figure 4.

Heatmap representation of reproducible Spearman’s correlation coefficients, across different cancer types, integrating the results shown in Figure 2 to Fig. 5 on breast, colorectal, non-small cell lung cancer and melanoma.

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.[72] 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[73] or the role of CXCR3 and its ligands CXCL9 to 11 in the recruitment of CD8+ T cells, Th1 cells and NK cells.[18] 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.[74]

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,[32] Hatzis,[33] Tabchy,[34] and Korde.[46] We also included the METABRIC dataset[81] 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.

Table 2.

Combined p-values are obtained by applying Fisher’s method.

Breast cancer dataset	Criterion of positive prognosis	P-value of exact Fisher test, for positively correlated, reproducible mRNAs coding for chemotactic factors/receptors	P-value of exact Fisher test, for positively correlated, reproducible mRNAs coding for chemotactic factors/receptors,with no low expression	P-value of exact Fisher test, for positively correlated, breast cancer-specific mRNAs coding for chemotactic factors/receptors	P-value of exact Fisher test, for other chemotactic factors/receptors
Bonnefoi	Pathological complete response	5.65e-16	3.95e-5	0.0258	1
Hatzis	Pathological complete response	0.00118	0.00104	0.443	1
Tabchy	Pathological complete response	1.24e-6	1.08e-6	0.0200	1
Korde	Treatment response	1	1	1	0.227
METABRIC	Alive	0.0141	0.0137	1	0.00183
Combined		6.56e-21	3.61e-21	0.0797	0.247

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.

Material and methods

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[27] 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).

Filtering for low expression

The first filter is based on microarray expression among the five cancer types. For each dataset, we computed the average (among different patients) of the expression rank, for each chemotactic factors/receptors. Because no obvious threshold emerged from these mean rank expressions, we choose it arbitrarily (value = 0.2) to filter out a reasonable amount of these chemotactic factors/receptors that are expressed at a low level.

The second filter is based on mRNA expression from the Human Protein Atlas [50,51,84] from the three methods (HPA, GTEx, FANTOM5), in the spleen and lymph node. Again, no obvious threshold emerged from these data. Therefore we arbitrarily considered that a gene has a low expression if its mRNA abundance is lower than (TPM/RPKM = 1), as determined by all three methods (HPA, GTEx, FANTOM5) and in both spleen and lymph node.

We applied these two filters to the list of chemotactic factors/receptors. We also considered a threshold of 2 for (TPM/RKM), applying the same analysis.

Fisher’s exact test for treatment response

We considered genes that were overexpressed in responsive tumors (or alive patients in METABRIC dataset), using a p-value threshold of 5%, a logFC threshold of log2(1.2) for the METABRIC, Bonnefoi and Korde datasets, a logFC threshold of log2(1.5) for Hatzis’ and Tabchy’s datasets (these thresholds varied in order to obtain non-empty intersections with the reproducible correlation gene list of Figure 4), using limma R-package.[85] For each dataset, we computed the contingency tables regarding overexpressed genes and reproducible correlation gene list of Figure 4 and applied Fisher’s exact tests. Significant results are shown in Table 2.