iScience ll OPEN ACCESS Article Immune Cell Associations with Cancer Risk Luis Palomero, Ivan Galván- Femenı́a, Rafael de Cid, ..., Antonis C. Antoniou, Conxi Lázaro, Miquel Angel Pujana aca20@medschl.cam.ac.uk (A.C.A.) clazaro@iconcologia.net (C.L.) mapujana@iconcologia.net (M.A.P.) HIGHLIGHTS Cancer risk genetic variants linked to immune/ stromal cell tissue content SH2B3 associated with BRCA1/2 cancer risk and immune cell counts Peripheral immune cell types linked to breast cancer age at diagnosis Palomero et al., iScience July 24, 101296 23, 2020 ª 2020 The Author(s). https://doi.org/10.1016/ j.isci.2020.101296 iScience ll OPEN ACCESS Article Immune Cell Associations with Cancer Risk Luis Palomero,1 Ivan Galván-Femenı́a,2 Rafael de Cid,2 Roderic Espı́n,1 Daniel R. Barnes,3 CIMBA,3 Eline Blommaert,1 Miguel Gil-Gil,4 Catalina Falo,4 Agostina Stradella,4 Dan Ouchi,5,13 Albert Roso-Llorach,5,13 Concepció Violan,5,13 Marı́a Peña-Chilet,6 Joaquı́n Dopazo,6 Ana Isabel Extremera,1 Mar Garcı́a-Valero,1 Carmen Herranz,1 Francesca Mateo,1 Elisabetta Mereu,7 Jonathan Beesley,8 Georgia Chenevix-Trench,8 Cecilia Roux,9 Tak Mak,9 Joan Brunet,10 Razq Hakem,11 Chiara Gorrini,9 Antonis C. Antoniou,3,* Conxi Lázaro,12,* and Miquel Angel Pujana1,14,* SUMMARY 1ProCURE, Catalan Institute Proper immune system function hinders cancer development, but little is known of Oncology, Bellvitge about whether genetic variants linked to cancer risk alter immune cells. Here, we Institute for Biomedical Research (IDIBELL), report 57 cancer risk loci associated with differences in immune and/or stromal L’Hospitalet del Llobregat, cell contents in the corresponding tissue. Predicted target genes show expres- Barcelona, Catalonia 08908, Spain sion and regulatory associations with immune features. Polygenic risk scores 2GCAT-Genomes for Life, also reveal associations with immune and/or stromal cell contents, and breast Germans Trias i Pujol Health cancer scores show consistent results in normal and tumor tissue. SH2B3 links pe- Sciences Research Institute (IGTP), Program for Predictive ripheral alterations of several immune cell types to the risk of this malignancy. and Personalized Medicine of Pleiotropic SH2B3 variants are associated with breast cancer risk in BRCA1/2 mu- Cancer (IMPPC), Badalona, Catalonia 08916, Spain tation carriers. A retrospective case-cohort study indicates a positive association 3Centre for Cancer Genetic between blood counts of basophils, leukocytes, and monocytes and age at breast Epidemiology, Department of Public Health and Primary cancer diagnosis. These findings broaden our knowledge of the role of the im- Care, University of mune system in cancer and highlight promising prevention strategies for individ- Cambridge, Cambridge CB1 8RN, UK uals at high risk. 4Department of Medical Oncology, Catalan Institute INTRODUCTION of Oncology, Bellvitge Institute for Biomedical The immune system maintains organismal integrity and function by continuously protecting itself from Research (IDIBELL), exogenous and endogenous assaults. The concept of ‘‘immunological surveillance of cancer’’ was first pro- L’Hospitalet del Llobregat, Barcelona, Catalonia 08908, posed by Burnet in 1970 (Burnet, 1970) and developed by Thomas about a decade later (Thomas, 1982). In Spain this theory, the immune system inactivates or eliminates cancer-prone cells that are detected early in 5Jordi Gol University Institute normal tissue (Ribatti, 2017). Although this idea remains a matter of debate, it is clear that some immune for Research Primary Healthcare (IDIAP Jordi Gol), factors decisively influence cancer development and progression. Immunosuppression due to primary im- Barcelona, Catalonia 08007, munodeficiency or due to therapies administered to prevent organ transplant rejection and certain virus Spain infections are associated with an increased risk of some cancers (Mortaz et al., 2016). In parallel, studies 6ClinicalBioinformatics Area, of mouse models with defined genetic alterations have demonstrated the relevance of immunosurveil- Fundación Progreso y Salud (FPS), Bioinformatics in Rare lance; for example, loss of Nkg2d, which encodes the activating receptor of natural killer (NK) and Diseases (BiER), CIBERER, T cells, increases the risk of spontaneous neoplasms (Guerra et al., 2008). INB-ELIXIR-es, Hospital Virgen del Rocı́o, Seville 41013, Spain Results from genome-wide association studies (GWASs) have identified risk-associated variants in loci cod- 7CNAG-CRG, Centre for ing for immune regulatory factors, such as NKG2D for cervical cancer risk (Chen et al., 2013). Indeed, Genomic Regulation, pathway-based analyses of GWAS results have highlighted the involvement of immune-related processes Barcelona Institute of Science in susceptibility to certain cancer types (Michailidou et al., 2017). In parallel, many germline genetic variants and Technology, Barcelona, Catalonia 08003, Spain can influence immune cell infiltration in tumors (Lim et al., 2018). Therefore, immune-centered investiga- 8Cancer Division, QIMR tions of normal or precancerous tissue could yield fundamental evidence for improving cancer risk estima- Berghofer Medical Research tion and prevention (Spira et al., 2017). However, whether common genetic variants linked to cancer risk Institute, Brisbane, QLD 4006, alter immune cell contents in the corresponding cancer target tissue, and/or at the systemic level, remains Australia 9Princess Margaret Cancer largely undetermined. Centre, The Campbell Family Institute for Breast Cancer The balance between immunological surveillance and tolerance is determined from a complex interplay Research, Ontario Cancer between different types of immune cells and other classes of stromal cells (Vinay et al., 2015; Gonzalez et al., 2018). Here, we describe an integrative analysis of genetic and transcriptome data from tissue Continued iScience July 24, 101296, 23, 2020 ª 2020 The Author(s). 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ll iScience OPEN ACCESS Article A B Cancer risk variant associated TCGA studies with immune/stromal cell tissue content Glioblastoma multiforme (GBM) High Head and neck squamous cell carcinoma (HNSC) content Thyroid carcinoma (THCA) Esophageal carcinoma (ESCA) stromal cells Breast invasive carcinoma (BRCA) Immune/ Lung adenocarcinoma (LUAD) Lung squamous cell carcinoma (LUSC) Liver hepatocellular carcinoma (LIHC) Pancreatic adenocarcinoma (PAAD) Skin cutaneous melanoma (SKCM) Kidney renal clear cell carcinoma (KIRC) Colon adenocarcinoma (COAD) Low Bladder urothelial carcinoma (BLCA) content TT TG GG Prostate adenocarcinoma (PRAD) Genotypes Uterine corpus endometrial carcinoma (UCEC) Ovarian cancer (OV) Testicular germ cell tumor (TGCT) C PRS associated with immune/stromal Blood Normal tumor- Primary tumor cell tissue content adjacent tissue tissue High Genetic data content RNAs ta G ta en RNAseq da da GWAS-identified GWAS-identified et ic ic et eq da cancer risk variants cancer risk variants stromal cells da en Immune/ ta G G Germline Somatic G ta TGATGCC genotypes genotypes TGATGCC Inferred cell tissue content T cells Macrophages Low B cells content NK cells Fibroblasts Low High ... risk PRS risk Figure 1. Strategy for Assessing the Effect of Immune/Stromal Cell Tissue Content on Cancer Risk (A) TCGA cancer projects analyzed in this study and data analysis workflow. (B) Association between gene expression-inferred immune/stromal cell tissue content and GWAS-identified risk variant. Institute, University Health (C) Association between gene expression-inferred immune/stromal cell tissue content and cancer PRS. Network, Toronto, ON M5G 2M9, Canada 10Hereditary Cancer defined as normal and located adjacent to surgically removed tumors, and from primary tumors analyzed Program, Catalan Institute of Oncology, Biomedical by The Cancer Genome Atlas (TCGA) (Cancer Genome Atlas Research Network et al., 2013). This enables Research Institute of Girona us to identify immune and stromal (hereafter ‘‘immune/stromal’’) cell tissue content associations with the (IDIBGI), Girona, Catalonia 17190, Spain risk of several human cancer types. Beyond single variants, polygenic risk scores (PRSs) also show associ- 11Princess Margaret Cancer ations with differences in inferred immune/stromal cell tissue contents. Consistent associations among im- Centre, Department of mune cell signatures, PRSs, and age at diagnosis suggest that higher immune cell infiltration reduces the Medical Biophysics, risk of breast cancer. We identify the lymphocyte SH2B adaptor 3 (LNK/SH2B3) locus as linking immune cell University Health Network and University of Toronto, counts and breast cancer risk, including that from BRCA1/2 mutation carriers. To evaluate this connection Toronto, ON M5G 2C1, further, we assess associations between breast cancer age at diagnosis and immune cell counts measured Canada at diagnosis in routine clinical blood tests; the results further suggest that peripheral immune cell status 12Hereditary Cancer influences breast cancer risk. Collectively, these findings may broaden our current knowledge of the bio- Program, Catalan Institute of Oncology, Oncobell, logical basis of cancer risk and thereby suggest strategies for cancer prevention. Bellvitge Institute for Biomedical Research (IDIBELL), L’Hospitalet del RESULTS Llobregat, Barcelona, Strategy to Evaluate Immune/Stromal Cell Tissue Contents that Influence Cancer Risk Catalonia 08908, and Spanish Biomedical Research TCGA has greatly increased our knowledge of human cancer through multilayer biological analyses, which Network Centre in Oncology include genetic and gene expression profiling of tissue considered to be normal and situated adjacent to (CIBERONC), Instituto de Salud Carlos III, Madrid the cancer (hereafter referred to as ‘‘normal’’) and primary tumors (Liu et al., 2018). In parallel, many suc- 28029, Spain cessful GWASs have identified hundreds of germline genetic variants associated with the risk of common 13Autonomous University of cancer types (Torkamani et al., 2018). By compiling GWAS results, we assigned cancer risk variants to 17 Barcelona, Bellaterra, TCGA projects based on tissue of origin correspondences (Figure 1A lists the cancer study acronyms, Catalonia 08913, Spain 14Lead Contact and Table S1 lists the cancer risk variants). As deconvolution analyses of bulk gene expression enable robust inference of cell type content in heterogeneous samples (Avila Cobos et al., 2018), deduced cell *Correspondence: aca20@medschl.cam.ac.uk content in normal tissue and tumors can be assessed for associations with cancer risk variants in multivar- (A.C.A.), iate analyses (Figure 1B). Moreover, as differences in cancer risk are more accurately defined by combina- clazaro@iconcologia.net tions of key variants in PRSs (Torkamani et al., 2018), it might be possible to better define the relevance of (C.L.), mapujana@iconcologia.net the immune/stromal cells by analyzing associations between PRSs and their corresponding signatures (Fig- (M.A.P.) ure 1C; Table S2 summarizes the number of normal tissue and primary tumor samples available for each https://doi.org/10.1016/j.isci. subsequent analysis). 2020.101296 2 iScience July 24, 101296, 23, 2020 iScience ll Article OPEN ACCESS To infer immune/stromal cell contents in normal and primary tumor samples using bulk tissue RNA sequencing (RNA-seq) data from TCGA, we applied a consensus-signature approach benchmarked against other methods (ConsensusTME; Jiménez-Sánchez et al., 2019). Using this approach, the computed immune/stromal estimations in the 17 identified TCGA datasets were typically found to be positively correlated with two other methods (Figure S1), as well as with estimates of leukocyte content measured by a different method (Taylor et al., 2018) (Figure S2). In turn, the estimates were generally found to be negatively correlated with aneuploidy (Figure S3), as expected (Taylor et al., 2018). In addition, the im- mune/stromal cell TCGA estimates showed significant differences between primary tumors with low or high levels of CD274/PDL1 and CD279/PDCD1 expression (Figures S4 and S5). Applying the method to RNA-seq data from whole blood samples of healthy adults also revealed positive correlations with im- mune cell enumerations using fluorescence-activated cell sorting (Newman et al., 2019) (Figure S6). Moreover, the estimates from this method were found to be highly correlated (Spearman’s r > 0.75) with the numbers of immune/stromal cells used to generate 100 pseudo-bulk breast tumors (Figure S7 and Methods). To further assess the validity of the inferred immune cell contents in TCGA, the deduced scores for each setting were assessed for associations with defined immune benchmark genes (Methods). In most set- tings, each inferred immune cell type content was found to be positively correlated with the expression of the assigned benchmark; the average Pearson’s correlation coefficient values for all signature-bench- mark pairs were 0.52 and 0.60 in the normal tissue and primary tumor sets, respectively (Figure S8 and Table S3). To assess further the coherence of the inferred immune cell contents, the corresponding scores were tested for association with the activity status of immune-related signaling pathways (Cubuk et al., 2018). This analysis revealed coherent clustering of immune and stromal cell types in normal tissue (Figure S9). Identification of Cancer Risk Variants Linked to Differences in Immune/Stromal Cell Tissue Content A total of 1,453 cancer risk variants were compiled from various sources; 214 of these were directly genotyped in TCGA, and the rest were imputed. After applying quality controls and filtering criteria (Methods), 627 and 966 variants were analyzed as potential immune/stromal quantitative trait loci (isQTLs). The isQTLs were identified using multivariate regressions including covariates of gender (when informative), age at diagnosis, tumor stage, and histology. The significance of the associations in each setting was concluded from 1,000 permutations (Methods). Tumor data were also analyzed because germline risk alleles are frequently associated with defined cancer histopathological and biomarker features (Michailidou et al., 2017). To avoid redundant tests, only cell signatures with eigenvalues >1 were examined in each setting (Table S4). Through this methodology, 22 signif- icant isQTLs were identified. These comprised normal tissue corresponding to esophageal carcinoma (ESCA), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), and uterine corpus endometrial carci- noma (UCEC), and primary tumors of breast cancer (BRCA), head and neck squamous cell carcinoma (HNSC), LUSC, and ovarian serous cystadenocarcinoma (OV) (Figure 2A and Table S5). Several of the identified isQTLs involved differences in endothelial and fibroblast cell content (Table S5), and these signals may also indicate links with immune cell differences: for instance, rs4072037 is associated with endothelial cell content in normal esophageal tissue, and this variant corresponds to a cis-expression (cis-e) QTL in several normal tissue types (GTEx Consortium, 2013) of genes whose products are functionally relevant in the immune system and biology of endothelial cells (Stenina-Adognravi, 2014), including GBA, GBAP1, TSP3/ THSB3, and MUC1, which are locus-mapped genes. In addition, the cancer pleiotropy rs11168936 variant (Feh- ringer et al., 2016) is associated with differences in fibroblast content in normal lung tissue corresponding to LUSC, and this variant is a cis-eQTL for C1QL4 in several normal tissue types (GTEx Consortium, 2013). Intrigu- ingly, C1q is a regulator of dendritic cell maturation (van Kooten et al., 2008), and we found this variant also to be associated with dendritic cell content (Figure 2B). In normal lung tissue corresponding to LUAD, rs17078110 is associated with B cells, and this locus codes for SASH1, a regulator of TLR4 signaling and cytokine production (Dauphinee et al., 2013). Among the isQTLs identified in tumors, rs3764419 is associated with cytotoxic cell con- tent in OV. This variant is a cis-eQTL for ATAD5 (GTEx Consortium, 2013), whose product is essential for proper B cell biology and immunoglobulin production (Zanotti et al., 2015). Overall, these data suggest that some can- cer risk variants are associated with immune/stromal cell tissue content, and that this link is mediated by alter- ations in genes of functional importance to the immune system. iScience July 24, 101296, 23, 2020 3 ll iScience OPEN ACCESS Article A B isQTL immune cell tissue contents Normal lung tissue Chr 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 0.40 BRCA Dendritic cell score ESCA 6 HNSC 0.30 LUAD LUSC OV LOD score 0.20 4 UCEC Normal tissue Primary tumor 0.10 2 CC CT TT rs11168936 0 C isQTL immune system functions Chr 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X BLCA BRCA 6 COAD ESCA HNSC LIHC LOD score 4 LUAD LUSC PADD PRAD 2 OV TGCT Normal tissue Primary tumor 0 Figure 2. isQTLs Linked to Cancer Risk (A) Graph indicating the association values (LOD scores) and relative chromosome locations of 22 isQTLs identified in normal tissue and primary tumor analyses (Table S5). (B) isQTL rs11168936 for dendritic cell content in normal tissue corresponding to LUSC. (C) Graph showing the LOD scores and relative chromosome locations of isQTLs identified using five TCGA-defined immune signatures (Table S6). Identification of Cancer Risk Variants Associated with Immune System Functions in Target Tissue and Primary Tumors The TCGA consortium examined 160 immune system-related gene expression signatures across hundreds of tumors and identified five of them as being informative for cancer classification: IFN-g response, lympho- cyte infiltration score, macrophage regulation, TGF-b response, and wound healing (Thorsson et al., 2018). Therefore, we sought to expand on the aforementioned cell-type-based associations by analyzing these five additional signatures using the same method as introduced earlier: multivariate regression with signif- icance determined from 1,000 permutations. This study identified 75 isQTLs, of which 11 variants had been identified in the previous isQTL analyses, which represents a significant concordance (Fisher’s exact test, p < 0.0001; Figure 2C and Table S6). Taking both analyses into account suggests that the risk of 13 cancer types may be influenced by immune/stromal cell tissue content. Of the 57 unique variants identified from all isQTLs, five were linked to tumor suppressor genes with recog- nized roles in the immune system: CDKN2A/B, DCC, MUC1, and SASH1. In addition, genomic enhancers identified in T helper, regulatory, effector, memory, and mononuclear cells were significantly over-repre- sented in this unique variant set relative to all human variants: > 2-fold enrichments, binomial test p values <0.05 (Ward and Kellis, 2012). Consistent with this observation, 8 (14%) variants corresponded to expres- sion (e) QTLs from 18 immune-related genes in normal human tissue (GTEx Consortium, 2013) and 13 (25%) corresponded to eQTLs identified in CD4+ and/or CD8+ T cells (Kasela et al., 2017) (Tables S5 and S6). To evaluate the relevance of these observed percentages, we examined the expected proportions when considering all cancer risk variants studied; lower percentages were identified in both analyses, with expected proportions of 11% (115/1,079) for eQTLs of immune-related genes in normal human tissue (GTEx Consortium, 2013) and 14% (151/1,079) for eQTLs in CD4/8+ T cells (Kasela et al., 2017). We then examined whether the eQTL gene targets documented within the isQTLs were functionally coherent by determining the proportion of significant gene expression-immune/stromal cell signature correlations and comparing the results with those from equivalent 1,000 random gene sets. Both isQTL sets (Tables S5 and S6) included 4 iScience July 24, 101296, 23, 2020 iScience ll Article OPEN ACCESS a higher proportion of eQTL gene targets that were positively correlated with immune/stromal cell signa- tures than expected by chance (Figure S10). Finally, variants correlated (r2 > 0.8) with each isQTL were in- tersected with various functional genomic data from B cells, monocytes, and CD4+ and CD8+ T cells, and for potential effects on protein coding sequences (Methods). These analyses identified two additional candidate genes (LIF and OSM) with established functions in the immune system, being involved in cyto- kine signaling (Table S7). Together, these data indicate that a substantial proportion of the isQTLs identi- fied influence genes whose expression is associated with immune system functions. PRS Associations with Immune/Stromal Cell Tissue Content Highlight Breast Cancer Risk The effects of individual cancer risk variants are generally small, but their combinations within PRSs can potentially identify individuals who are at substantially higher risk than average for the population (Torka- mani et al., 2018). Therefore, reported PRSs were computed in the corresponding normal tissue and pri- mary tumor TCGA settings and evaluated for associations with immune/stromal cell contents using multi- variate analyses as described earlier. The study of normal tissue was limited to breast. Despite the valuable TCGA resource, the available sample size sets limited the detection of nominal significant associations to those with correlation coefficients of r > 0.3 in normal breast and of r > 0.12 in BRCA; higher correlations would be required for all other normal or tumor settings (Figure S11). In normal breast, most immune/stromal cell contents tended to be negatively correlated with the corre- sponding PRSs; the PRS cell signature correlation coefficients for overall and estrogen receptor (ER)-pos- itive breast cancer were significantly less than zero (p values <0.001; Figure 3A). The ER-negative PRS could not be computed because of the relatively low number of normal samples of this subtype and with com- plete data. Analogous limitations were encountered when attempting to analyze triple-negative breast cancer (TNBC) and human epidermal growth factor receptor 2 (HER2)-positive breast cancers, and there were no HER2-specific PRSs to analyze whatsoever. Potentially protective cell types (i.e., those exhibiting a nominally significant negative correlation between cell content and PRS) in the aforementioned two breast cancer settings included dendritic cells, eosinophils, macrophage M2, monocytes, neutrophils, and T cell terminal differentiation (Figure 3A). In addition to the breast cancer PRSs, eight other scores (Fritsche et al., 2018) were examined in their cor- responding primary tumor TCGA settings. The distribution of the correlation coefficients between im- mune/stromal cell tissue content and the PRS was again found to be less than zero not only in BRCA but also in in glioblastoma multiforme (GBM; with a major contribution for fibroblast content) and thyroid car- cinoma (THCA; Figure 3B). Conversely, positive correlations were detected in bladder urothelial carcinoma (BLCA), OV, prostate adenocarcinoma (PRAD), skin cutaneous melanoma (SKCM), and, principally, in LUAD and LUSC (Figure 3B). Conversely, positive correlations were detected in BLCA, serous OV, PRAD, SKCM, and, principally, in LUAD and LUSC (Figure 3B). Therefore, risk stratification based on PRSs may also be linked to differences in immune/stromal cell content in normal and/or tumor tissue. LUAD and LUSC PRSs shared positive correlations (p < 0.05) with cytotoxic and NK cell tissue contents; however, these as- sociations may be influenced by smoking status, because LUAD current smokers showed an opposite trend (Figure S12). Combined analyses of normal tissue and primary tumor data further suggested common protective effects for high immune cell content in breast and colorectal tissue, and also potentially in brain and a few other settings (Figure 3C). In contrast, high immune cell content might principally increase the risk of lung, bladder, and pancreatic cancer (Figure 3C), although, as already noted, smoking may influence these as- sociations. Then, analyses of COAD subtypes (Methods) suggested protective effects for high immune cell content in genomic stable tumors (Figure 3D, left panel), but this association might be biased due to PRS development in overall incident cases. When analyzing the COAD molecular subtypes, lower risk of CSM3 might also be associated with higher immune cell content (Figure 3D, right panel). The sample sets of these subtype analyses were relatively small to obtain robust conclusions, but, when compared with normal colorectal tissue, an opposite trend was observed (Figure 3D, left panel), which suggests that immune cell infiltration has different roles between normal tissue and tumors. As described earlier, the normal breast and BRCA settings both showed PRS-cell signature negative cor- relations. To assess these observations further, the correlation estimates were compared with those from similar analyses using age at diagnosis instead of the PRSs. In normal breast tissue, the iScience July 24, 101296, 23, 2020 5 ll iScience OPEN ACCESS Article A B PRSs - Cell signatures PRSs - Cell signatures Normal breast tissue Primary tumors CD8 T cell activation Fibroblasts Gamma-delta T cells Endothelial Fibroblasts B cells B cells Mast cells T regulatory cells Eosinophils Type I interferon response Neutrophils T-Reg CD8 T cells Type II interferon response CD4 T cells Plasma cells Macrophage M1 Macrophages Dendritic cells Eosinophils Anergy CD4 T cells Monocytes CD8 T cells Macrophages Anergy Macrophage M2 NK cells Cytotoxic cells Macrophage M2 polarization Plasma cells Dendritic cells T cell terminal differentiation Cytotoxic cells T regulatory cells Endothelial Gamma-delta T cells Anti-inflammatory Cytolytic effector pathway Cytolytic effector pathway NK cells Neutrophils CD8 T cell activation Mast cells Type I interferon response Macrophage M2 Type II interferon response Macrophage M1 polarization Macrophage M1 polarization Macrophage M1 Anti-inflammatory T cell terminal differentiation T-Reg Monocytes Macrophage M2 polarization E 78) CA 7) PR S-2 A o BM 03) 1 B RC eral (87) R+ 23) 2 S 1 PR ll (5 ) ou D ( 7) V 1) PR PRS BL M ( ) -2 CA 55) 1 S CA E D (1 ) ou + ( 6) V 8) 8) (2 ) ) 2 8 0 SC 78 09 t-score P < 0.05 3) ) ) 2) t-score P < 0.05 65 66 TH R- (7 PR PRS ove (41 1 s O 13 SK (14 PR 1 BR PR (16 R 4 s O 40 LU (14 4 (8 (8 4 ( G (2 RC A E l (5 LU D (2 ( ( 2 RC R- 2 ll R+ R+ ll C Not tested ra ra 1 B CA E 0 ra 0 A A A 1E 2E ve ve A -2 -2 B v -4 2o 1o -4 S- S- PR 2 BR A er er - PR PR S- S- C PR PR R S- S- 2B S- S- S- PR PR S- S- PR C PRSs - Cell signatures D Normal tissue and primary tumors COAD PRSs - Subtypes Fibroblasts Cytotoxic cells Endothelial B cells Anergy Fibroblasts Monocytes Endothelial CD8 T cells Mast cells Eosinophils T regulatory cells CD4 T cells NK cells Macrophages M2 Mast cells Macrophages Gamma-delta T cells Gamma-delta T cells NK cells Eosinophils Macrophages M1 Cytotoxic cells Dendritic cells CD8 T cells Dendritic cells Monocytes Macrophage M2 Plasma cells Neutrophils Macrophages Eosinophils Dendritic cells Fibrob lasts Neutrophils Macrophage M1 Macrophages M2 Cytotoxic cells CD8 T cell activation Type II interferon response Monocytes NK cells Plasma_cells T cell terminal differentiation Neutrophils Plasma_cells Cytolytic effector pathway Mast cells Macrophage M1 polarization Macrophages T-Reg Endothelial CD4 T cells Anti-inflammatory Macrophage M2 polarization Macrophages M1 T regulatory cells Type I interferon response B cells B cells ou AD 50) RC CA C ( 3) ov + ( ) 2 P CA 2) 1 S L AD 9) V 3) S- s O M (5 ) PR V (1 7) 1 B R KIR (14 ) PR BR THC ll (57 ) S- CA A 9) RC R+ 1) ) CA BM 8) RC CO - ) ov D ( 8) all 7) ) t-score P < 0.05 A ER 209 1 AD 51 2 ER 2 7 79 CD4 T cells CD8 T cells 4 R (16 s O (32 er SKC (15 er 45 3 A (45 G - (7 (8 A (7 er 19 4 1 BL C (2 1 B E (2 (5 ( 2 T regulatory cells Gamma-delta T cells R P < 0.01 E a 0 S LU U MS 15) CI (33) GS 2) 4) S4 ) S1 ) S3 ) 4) -2 CM (56 CM (54 CM (30 t-score P < 0.05 PR ou 2 (2 (2 -4 R e( (1 1 4 er A S- 2B I S2 su N 2 S- 1 B 2S PR CM S- 0 tis 2 2B S- S- PR RS- S- PR -2 al PR PR PR -4 S- rm P PR No E Age at diagnosis - Cell signatures F Normal tissue Normal breast tissue Anergy 1.0 PRS Signature - PRS correlation Anti-inflammatory B cells CD8 T cell activation Cytolytic effector pathway Overall Cytotoxic cells Dendritic cells 0 ER+ Endothelial Eosinophils Fibroblasts (t-score) Macrophage M1 polarization Macrophage M2 polarization Macrophages -1.0 Macrophage M1 Macrophage M2 Mast cells Monocytes Neutrophils NK cells Plasma cells -2.0 T cell terminal differentiation CD4 T cells CD8 T cells r = -0.30 Gamma-delta T cells T-Reg r = -0.45 T regulatory cells -3.0 Type I interferon response Type II interferon response -1.0 0 1.0 2.0 3.0 4.0 5.0 Signature - Age at diagnosis correlation RC ) AD ) 2) AD 2) CA ) AD ) SC ) LI (49) HN (42) CA 3) t-score FDR < 5% KI (111 PR (72 TH (58 CO (58 LU (40 (5 (2 (4 4 (t-score) EC HC SC 2 0 Not tested UC LU BR -2 -4 Figure 3. Associations between Immune/Stromal Cell Signatures and PRSs (A) Unsupervised clustering of the results of the regression analysis between cell signatures and PRSs in normal breast tissue. The y axis depicts the cell type signatures, and the x axis shows the PRSs. Sources #1 and #2 of the PRSs are detailed in Methods. ER+ and ER indicate estrogen receptor-positive and estrogen receptor-negative subsets, respectively. The maximum sample size used in each analysis is shown. The color scale (t-score) is calculated as the b estimate divided by the standard error. Nominally significant associations are indicated by black-outlined rectangles. (B) Unsupervised clustering of the coefficients of the regression of cell signature values in primary tumor TCGA studies and the corresponding PRSs. The gray-filled rectangles indicate ‘‘not tested’’ correlations because the corresponding cell signatures were only defined for breast cancer. (C) Unsupervised clustering of the coefficients of the regression of cell signature values in combined normal tissue and primary tumor datasets, and the corresponding PRSs. The regression p values <0.01 are also indicated as depicted in the inset. 6 iScience July 24, 101296, 23, 2020 iScience ll Article OPEN ACCESS Figure 3. Continued (D) Unsupervised clustering of the coefficients of the regression analysis between PRSs and cell signatures in combined normal tissue and primary tumors of the COAD study, divided by cancer subtypes. (E) Unsupervised clustering of the results of the regression analysis of cell signatures in normal tissue and age at diagnosis across TCGA studies. Associations significant at a false discovery rate (FDR) < 5% are indicated by black-outlined rectangles. (F) Negative correlations between the coefficients of regressions of immune/stromal cell contents and age at diagnosis or PRSs in normal breast tissue, for all cases and only ER-positive cases. The correlation coefficients are shown. immune/stromal cell contents tended to show positive correlations with age at diagnosis (p < 0.001; Fig- ure 3E). Consequently, negative correlations were detected between the estimates from the two parallel analyses, considering all cases or solely ER-positive cases (Figure 3F). Therefore, relatively higher im- mune/stromal cell content in normal breast might be a factor protecting against development of malignancy. SH2B3 Connects Immune Cell Tissue Content with Breast Cancer Risk The identified cancer risk isQTLs could be explained by peripheral alterations in immune cells. Examination of GWAS results for blood cell traits revealed that the tumor COAD isQTL rs12412391 in chromosome 10 (Table S6) is in linkage disequilibrium (r2 = 0.93) with rs11190133, which is associated with differences in platelets in the UK Biobank study (Astle et al., 2016). These variants constitute an eQTL of NKX2-3 (GTEx Consortium, 2013), and, remarkably, loss of the mouse ortholog causes developmental alterations in the spleen, colonic crypts, and lymphocyte tissue homing (Pabst et al., 1999). In addition to this locus, the tumor BRCA isQTL rs11065979 in chromosome 12 (Table S5) was associated with blood count differences in ba- sophils, erythrocytes, eosinophils, leukocytes, monocytes, and neutrophils in the UK Biobank study (Astle et al., 2016) (Table S8). The same study also indicated an association with breast cancer risk (p = 0.0003; Table S8). This variant has also been linked to cancer pleiotropy (Fehringer et al., 2016) and psoriasis (Tsoi et al., 2017), among other traits (GWAS Catalog). A variant in linkage disequilibrium, rs3184504 (r2 = 0.89), had also been associated with breast cancer risk (Fehringer et al., 2016), serum IgA levels (Jons- son et al., 2017), and various autoimmune diseases (Webb and Hirschfield, 2016), among other traits (GWAS Catalog). To investigate further the role of the isQTL identified in chromosome 12 and linked to breast cancer risk, we analyzed association results from BRCA1/2 mutation carriers. Both depicted variants showed nominal as- sociations with breast cancer risk in women carriers of germline BRCA1 or BRCA2 mutations: BRCA1 mu- tation carriers, rs11065979 hazard ratio (HR) = 0.96, 95% confidence interval (CI) 0.92–0.99, p = 0.018; rs3184504 HR = 0.95, 95% CI 0.92–0.99, p = 0.006; BRCA2 mutation carriers, rs11065979 HR = 0.94, 95% CI 0.90–0.99, p = 0.019; and rs3184504 HR = 0.93, 95% CI 0.89–0.98, p = 0.003. Then, wider examination of this region in chromosome 12 identified several genetic associations (p < 0.01) with breast and/or ovarian cancer risk in these women (Figure 4A and Table S9). The chromosome 12 locus identified here includes many eQTL signals for SH2B3 in EBV-transformed lym- phocytes and normal tissue (Figure 4A, bottom panel). Next, to evaluate potential causality linked to SH2B3, complementary gene expression analyses were performed using the normal breast tissue TCGA data. First, the expression of SH2B3 was found to be positively correlated with most of the immune cell/ stromal cell signatures (Figure 4B); second, SH2B3 expression was also found to be positively correlated with age at diagnosis, adjusted for tumor stage and regardless of cancer subtype (Figure 4C); third, an 84-gene signature corresponding to gene and protein functional relationships with mouse Sh2b3 and/or human SH2B3 (Huan et al., 2015) was also positively correlated with age at diagnosis (Figure 4D); and last, SH2B3 expression was positively correlated with the protein measures of CD26, cell surface glycopro- tein receptor important for T cell activation (Klemann et al., 2016), and TFCR, transferrin receptor required for erythropoiesis and immune system development (Jabara et al., 2016) (Figure 4E). In addition, the asso- ciation between SH2B3 expression in normal breast and age at diagnosis was replicated in an independent dataset (Terunuma et al., 2014): n = 47, r = 0.30, p = 0.039. Therefore, an identified isQTL may influence breast cancer risk through perturbation of SH2B3 expression, which is expected to be fundamental for ac- curate systemic development and function of immune cell populations (Li et al., 2000; Velazquez et al., 2002; Jabara et al., 2016). iScience July 24, 101296, 23, 2020 7 ll iScience OPEN ACCESS Article A B 3.5 BRCA1 breast cancer Normal breast tissue BRCA1 ovarian cancer SH2B3 - Cell signatures BRCA2 breast cancer Endothelial BRCA2 ovarian cancer Macrophages M1 Dendritic cells 3.0 M2 macrophage polarization Macrophages M2 -log10 p value Macrophages Monocytes 1.0 Mast cells Fibroblasts rs3184504 Neutrophils M1 macrophage polarization 2.5 T-Reg T cell terminal differentiation Eosinophils Cytolytics effector pathway T cells CD8 NK cells Cytotoxic cells 111,800 112,190 0.0 T cells gamma delta 2.0 T cells CD4 rs11065979 chr12 (hg19) Type II interferon response Type I Iinterferon response CD8 T cell activation PCC B cells 1.0 0.5 T regulatory cells Plasma cells PHETA1 SH2B3 BRAP 0.0 -0.5 Anergy Anti-inflammatory Genes -1.0 ATXN2 ACAD10 -1.0 TSS TES EBV-transformed eQTL SH2B3 lymphocytes Lung NES 1 Salivary grand 0.5 Pancreas 0.2 Spleen 0 -0.2 Whole blood -0.5 -1 Promoter 3’-UTR C D E Normal breast tissue Normal breast tissue Normal breast tissue CD26 1 TFRC P = 0.006 HER2 HER2 2 P = 0.001 P = 0.006 Luminal Luminal SH2B3-functionally 0.30 SH2B3 expression related gene set TNBC TNBC 1 RPPA normalized value All 0.20 All 0 0 0.1 -1 0 -1 -2 -0.10 P = 0.037 -3 40 60 80 40 60 80 -2 Age at diagnosis Age at diagnosis -3 -3 -2 -1 0 1 2 SH2B3 normalized expression Figure 4. The SH2B3 Locus Shows Associations with Breast Cancer Risk and Immune Cell Features (A) Graph showing the chromosome 12 association results (-log10 p value, y axis) with breast and ovarian cancer risk (as depicted in the inset) in women carriers of BRCA1/2 mutations. The rs3184504 and rs11065979 variants are indicated. (B) Rank of expression correlations (Pearson’s correlation coefficients [PCCs]) between SH2B3 and immune/cell signatures in normal breast. All PCCs were >0, but three of them did not reach nominal significance (marked gray). (C) Positive correlation between SH2B3 expression in normal breast and age of diagnosis of breast cancer. The trend lines for all cases and subtypes (luminal, HER2-positive, and triple-negative breast cancer [TNBC]) are shown. The correlation p value from the multivariate regression analysis is shown. (D) Positive correlation between SH2B3 functionally related gene set in normal breast and age of diagnosis of breast cancer. (E) Positive correlation between SH2B3 expression and CD26 and TFRC protein expression as measured by TCGA reverse-phase protein array (RPPA) assays. The correlation p value from the multivariate regression analysis is shown. Peripheral Immune Cell Counts Are Associated with Breast Cancer Risk To assess the proposed link between breast cancer risk and peripheral immune cell counts, which in turn might be influenced by specific genetic variants and gene candidates, a retrospective case-cohort study was performed. Data on age at diagnosis, tumor stage and subtype, and blood test results from 259 breast cancer cases were compiled in a tertiary referral hospital (Methods). The cases were randomly selected from clinical health records and showed an average age at diagnosis of 55.6 years, 95% CI 54.0–57.1 years. The blood test data were those collected on the date closest to diagnosis: 6 patients had the blood test on the same date as their diagnosis, 40 were earlier (on average 40 days before), and 182 were later (on average 45 days later): the average time between the blood test and disease diagnosis was 23.9 days, 95% CI 16.3–31.5 days. A multivariate regression analysis including tumor stage and subtype revealed three 8 iScience July 24, 101296, 23, 2020 iScience ll Article OPEN ACCESS Basophils Leukocytes Monocytes P = 0.047 P = 0.024 P = 0.009 Tumor stage 0-2 4.5 1.25 Tumor stage 3-4 Luminal Luminal+HER2 HER2 TNBC 4.0 All 0.2 1.00 Normalized cell counts (log2) 3.5 0.75 3.0 0.1 0.50 2.5 0.25 0.0 2.0 20 40 60 80 100 20 40 60 80 100 20 40 60 80 100 Age at diagnosis Age at diagnosis Age at diagnosis Figure 5. Positive Correlation between Peripheral Immune Cell Counts around Time of Diagnosis and Breast Cancer Age at Diagnosis Positive correlations between basophil, leukocyte, and monocyte blood counts and age of diagnosis of breast cancer in a retrospective case-cohort hospital-based study. The value of p shown here is that associated with the coefficient calculated as part of the multivariate regression analysis. The trends for tumor subgroups are shown as depicted in the inset. immune cell types to be significantly (p < 0.05) positively correlated with age at diagnosis: basophils, leu- kocytes, and monocytes (Figure 5). The trends were consistent for different tumor stages (0–2 and 3–4) and the major cancer subtype (i.e., luminal); the smaller patient sets of HER2-postive breast cancer (n = 18) and TNBC (n = 17) showed greater variability (Figure 5). The neutrophil to lymphocyte ratio, which is an estab- lished rate associated with breast cancer prognosis (Ethier et al., 2017), was not found to be associated with age at diagnosis in this study (p = 0.65). DISCUSSION The results of this study support the idea that the risk of certain cancers is influenced by the content of im- mune/stromal cells in the target tissue and/or by differences in peripheral immune cell counts. Of the 17 cancer settings analyzed, 57 risk loci comprising 13 cancer types were associated with differences in im- mune/stromal cell content with respect to the corresponding normal tissue and/or primary tumors. The gene candidates linked to these associations include several with key functions in the immune system, and they show significant enrichments in immune-related regulatory features and expression profiles. Detection of associations between immune/stromal cell signatures and PRSs provide further evidence that differences in these cell contents influence cancer risk. Nevertheless, the role of some cell types is multifaceted; for example, endothelial cells can regulate trafficking and activation of immune cells in a given tissue, but, critically, also determine angiogenesis (Hendry et al., 2016). Similarly, a given genetic variant may influence the expression of more than one gene target and/or indirectly alter the immune sys- tem by different mechanisms, such as by provoking oncogenic-induced signals. Unexpectedly, there appear to be opposing cancer risk effects for immune cell contents across cancer types. These might be due to differences in tissue microenvironment conditions, such as inflammation caused by smoking or other factors (Shalapour and Karin, 2015). However, the study was limited by the rela- tively low numbers of normal tissue samples available for analysis, and, potentially, by gene expression al- terations in normal tissue adjacent to neoplasms. This study had more power to detect significant results in normal breast tissue and BRCA and, consequently, the results prove to be more relevant and coherent in iScience July 24, 101296, 23, 2020 9 ll iScience OPEN ACCESS Article these settings. Carrying out similar analyses in other normal and cancer tissue contexts would appear to be worthwhile. Such additional studies would benefit from complementary molecular marker and signaling analyses, which would definitively establish the functional consequences of inferred cell alterations. At the same time as providing insight into the biological basis of cancer initiation, this study yields data that could be useful for analyses of cancer risk and prevention. Associations of PRSs with immune cell signatures could inform preventive strategies by modulating specific cell functions and/or their signaling molecules in individuals at high risk (Spira et al., 2017). This idea is particularly relevant in breast cancer. Our study shows consistent associations between immune/stromal cell signatures and breast cancer PRSs or age at diagnosis in normal tissue. A recognized risk locus connects differences in most peripheral immune cell types to breast cancer risk. This locus harbors the SH2B3 gene, which is altered in hematological neoplasms and autoim- mune diseases (Maslah et al., 2017, p. 3). Common genetic variation at this locus has been linked to cancer pleiotropy, including breast cancer susceptibility (Hung et al., 2015; Fehringer et al., 2016). We extend these observations by identifying potential associations with breast and ovarian cancer risk in BRCA1/2 mutation carriers. Our study shows consistent expression correlations of SH2B3 or SH2B3 functionally related genes with age at diagnosis using normal breast tissue data. Thus, pharmacological enhancement of SH2B3 func- tion might reduce cancer risk in individuals with high PRSs and/or carriers of BRCA1/2 mutations. However, the functional impact on SH2B3 remains to be established, and, therefore, prospective studies determining the expression and/or functional differences of SH2B3 among individuals with specific alleles in the corre- sponding locus, and their associations with peripheral immune cell counts and cancer risk, are needed. The effect of systemic differences of immune cell counts on breast cancer risk is further supported by un- expected associations between basophil, leukocyte, and monocyte blood counts and age at diagnosis from a retrospective case-cohort study. Relatively low monocyte counts collected over a 1-year period of disease diagnosis have recently been associated with increased breast cancer risk (Kresovich et al., 2020). However, high baseline leukocyte counts in a prospective study of postmenopausal women were found to be associated with increased breast cancer incidence (Margolis et al., 2007). In our study, we aimed to assess whether individuals’ status of having relatively low peripheral immune cell counts was asso- ciated with initial cancer development, hypothetically due to reduced immunosurveillance. Our results are consistent with this explanation, and among other factors, altered SH2B3 function might give rise to these observations. As a whole, the results of this study may be useful for improving cancer risk estimation, and for identifying preventive approaches. Limitations of the Study The present report identifies cancer-associated genetic variants and polygenic risk scores linked to the alteration of immune and/or stromal cell systemic and/or tissue contents. These links could explain the greater cancer risk. However, the study has several limitations that should be borne in mind. The cell con- tent inferences were based on gene expression profiles, and therefore, molecular and cellular analyses are required to corroborate them and accurately assess their functional consequences. The observed associ- ations could also be indirect in some instances. The study was also limited by the original sample sets, and observed associations could be confounded by other factors, such as the level of tissue inflammation, in- dividual hormonal status, and lifestyle aspects. The genetic basis of the proposed associations between blood cell count and age at breast cancer diagnosis in the studied cohort remains unknown, and it is un- clear whether similar associations exist in BRCA1/2 mutation carriers. Resource Availability Lead Contact Further information and requests for resources and data should be directed to and fulfilled by the Lead Contact Miquel Angel Pujana (mapujana@iconcologia.net). Materials Availability No materials were generated. Data and Code Availability TCGA data were obtained from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov) and from the corresponding consortium publications. Individual genetic data were obtained following 10 iScience July 24, 101296, 23, 2020 iScience ll Article OPEN ACCESS specific approval: dbGaP Data Access Committee project #11689. The R software algorithms developed by others and applied in this study are detailed in the Transparent Methods supplemental file. A complete pipeline was implemented and is available at https://github.com/pujana-lab/systematicQTL/. METHODS All methods can be found in the accompanying Transparent Methods supplemental file. SUPPLEMENTAL INFORMATION Supplemental Information can be found online at https://doi.org/10.1016/j.isci.2020.101296. ACKNOWLEDGMENTS The results presented here are partly based on data generated by the TCGA Research Network (https:// www.cancer.gov/tcga), and we would like to express our gratitude to the TCGA consortia and their coor- dinators for providing the data and clinical information. This study was supported by the following patient foundations in Catalonia: DACMA, GINKGO, Sosciathlon, and ‘‘Viladecans Contra el Cáncer’’. The work was also supported by grants from the Carlos III Institute of Health funded by FEDER funds – a way to build Europe – (ISCIII; Ministry of Science, Innovation and Universities; PI16/00563, PI18/01029, PI19/00553, and CIBERONC), the Generalitat de Catalunya (SGR 2017-449 and 2017-1282; and PERIS MedPerCan and URD- Cat), and the CERCA Program. AUTHOR CONTRIBUTIONS L.P. and M.A.P. designed and performed the experiments. I.G.-F. and R.d.C. performed variant imputation in TCGA. L.P., R.E., E.B, and M.G.-V. performed bioinformatic analyses of genetic variants, PRSs, gene expression, and cell signatures. D.R.B. and A.C.A. performed association studies in BRCA1/2 mutation car- riers. J. Beesley and G.C.-T. contributed to the association studies and performed candidate target gene analyses. M.G.-G., C.F., A.S., and A.I.E. contributed to the cohort study. D.O., A.R.-L., and C.V. contributed to the data normalization and multivariate regression designs. M.P.-C. and J.D. contributed to the pathway analyses. C.H., F.M., E.M., C.R., T.M., J. Brunet, R.H., and C.G. contributed data, evaluated the results, and critically revised the manuscript. C.L. contributed to the study design. M.A.P. wrote the paper. DECLARATION OF INTERESTS M.A.P. is recipient of an unrestricted research grant from Roche Pharma for the development of the Pro- CURE ICO research program. C.F. received support from Pfizer unrelated to this study. Received: March 18, 2020 Revised: May 23, 2020 Accepted: June 15, 2020 Published: July 24, 2020 REFERENCES Astle, W.J., Elding, H., Jiang, T., Allen, D., Ruklisa, I., Sander, C., and Stuart, J.M. (2013). The cancer Veenstra, T.D., and Karsan, A. (2013). SASH1 is a D., Mann, A.L., Mead, D., Bouman, H., Riveros- genome atlas pan-cancer analysis project. Nat. scaffold molecule in endothelial TLR4 signaling. Mckay, F., Kostadima, M.A., et al. (2016). The Genet. 45, 1113–1120. J. Immunol. 191, 892–901. allelic landscape of human blood cell Trait variation and links to common complex disease. Chen, D., Juko-Pecirep, I., Hammer, J., Ivansson, Ethier, J.-L., Desautels, D., Templeton, A., Shah, Cell 167, 1415–1429.e19. E., Enroth, S., Gustavsson, I., Feuk, L., P.S., and Amir, E. (2017). Prognostic role of Magnusson, P.K.E., McKay, J.D., Wilander, E., neutrophil-to-lymphocyte ratio in breast cancer: a Avila Cobos, F., Vandesompele, J., Mestdagh, P., and Gyllensten, U. (2013). Genome-wide systematic review and meta-analysis. Breast and De Preter, K. (2018). Computational association study of susceptibility loci for cervical Cancer Res. 19, 2. deconvolution of transcriptomics data from cancer. J. Natl. Cancer Inst. 105, 624–633. mixed cell populations. Bioinformatics 34, 1969– Fehringer, G., Kraft, P., Pharoah, P.D., Eeles, R.A., 1979. Cubuk, C., Hidalgo, M.R., Amadoz, A., Pujana, Chatterjee, N., Schumacher, F.R., Schildkraut, M.A., Mateo, F., Herranz, C., Carbonell- J.M., Lindström, S., Brennan, P., Bickeböller, H., Burnet, F.M. (1970). The concept of Caballero, J., and Dopazo, J. (2018). Gene et al. (2016). Cross-cancer genome-wide analysis immunological surveillance. Prog. Exp. Tumor expression integration into pathway modules of lung, ovary, breast, prostate, and colorectal Res. 13, 1–27. reveals a pan-cancer metabolic landscape. cancer reveals novel pleiotropic associations. Cancer Res. 78, 6059–6072. Cancer Res. 76, 5103–5114. Cancer Genome Atlas Research Network, Weinstein, J.N., Collisson, E.A., Mills, G.B., Shaw, Dauphinee, S.M., Clayton, A., Hussainkhel, A., Fritsche, L.G., Gruber, S.B., Wu, Z., Schmidt, E.M., K.R.M., Ozenberger, B.A., Ellrott, K., Shmulevich, Yang, C., Park, Y.-J., Fuller, M.E., Blonder, J., Zawistowski, M., Moser, S.E., Blanc, V.M., iScience July 24, 101296, 23, 2020 11 ll iScience OPEN ACCESS Article Brummett, C.M., Kheterpal, S., Abecasis, G.R., van Kooten, C., Fiore, N., Trouw, L.A., Csomor, E., Shalapour, S., and Karin, M. (2015). Immunity, and Mukherjee, B. (2018). Association of Xu, W., Castellano, G., Daha, M.R., and inflammation, and cancer: an eternal fight polygenic risk scores for multiple cancers in a Gelderman, K.A. (2008). Complement production between good and evil. J. Clin. Invest. 125, 3347– phenome-wide study: results from the Michigan and regulation by dendritic cells: molecular 3355. Genomics Initiative. Am. J. Hum. Genet. 102, switches between tolerance and immunity. Mol. 1048–1061. Immunol. 45, 4064–4072. Spira, A., Yurgelun, M.B., Alexandrov, L., Rao, A., Bejar, R., Polyak, K., Giannakis, M., Shilatifard, A., Gonzalez, H., Hagerling, C., and Werb, Z. (2018). Kresovich, J.K., O’Brien, K.M., Xu, Z., Weinberg, Finn, O.J., Dhodapkar, M., et al. (2017). Precancer Roles of the immune system in cancer: from tumor C.R., Sandler, D.P., and Taylor, J.A. (2020). atlas to drive precision prevention trials. Cancer initiation to metastatic progression. Genes Dev. Prediagnostic immune cell profiles and breast Res. 77, 1510–1541. 32, 1267–1284. cancer. JAMA Netw. Open 3, e1919536. Stenina-Adognravi, O. (2014). Invoking the power GTEx Consortium (2013). The genotype-tissue Li, Y., He, X., Schembri-King, J., Jakes, S., and of thrombospondins: regulation of expression (GTEx) project. Nat. Genet. 45, Hayashi, J. (2000). Cloning and characterization of thrombospondins expression. Matrix Biol. 37, 580–585. human Lnk, an adaptor protein with pleckstrin 69–82. homology and Src homology 2 domains that can inhibit T cell activation. J. Immunol. 164, 5199– Taylor, A.M., Shih, J., Ha, G., Gao, G.F., Zhang, X., Guerra, N., Tan, Y.X., Joncker, N.T., Choy, A., Berger, A.C., Schumacher, S.E., Wang, C., Hu, H., Gallardo, F., Xiong, N., Knoblaugh, S., Cado, D., 5206. Liu, J., et al. (2018). Genomic and functional Greenberg, N.M., Greenberg, N.R., and Raulet, approaches to understanding cancer aneuploidy. D.H. (2008). NKG2D-deficient mice are defective Lim, Y.W., Chen-Harris, H., Mayba, O., Lianoglou, Cancer Cell 33, 676–689.e3. in tumor surveillance in models of spontaneous S., Wuster, A., Bhangale, T., Khan, Z., malignancy. Immunity 28, 571–580. Mariathasan, S., Daemen, A., Reeder, J., et al. Terunuma, A., Putluri, N., Mishra, P., Mathé, E.A., (2018). Germline genetic polymorphisms Dorsey, T.H., Yi, M., Wallace, T.A., Issaq, H.J., Hendry, S.A., Farnsworth, R.H., Solomon, B., influence tumor gene expression and immune Zhou, M., Killian, J.K., et al. (2014). MYC-driven Achen, M.G., Stacker, S.A., and Fox, S.B. (2016). cell infiltration. Proc. Natl. Acad. Sci. U S A 115, accumulation of 2-hydroxyglutarate is associated The role of the tumor vasculature in the host E11701–E11710. with breast cancer prognosis. J. Clin. Invest. 124, immune response: implications for therapeutic 398–412. strategies targeting the tumor Liu, J., Lichtenberg, T., Hoadley, K.A., Poisson, microenvironment. Front. Immunol. 7, 621. L.M., Lazar, A.J., Cherniack, A.D., Kovatich, A.J., Thomas, L. (1982). On immunosurveillance in Benz, C.C., Levine, D.A., Lee, A.V., et al. (2018). An human cancer. Yale J. Biol. Med. 55, 329–333. Huan, T., Meng, Q., Saleh, M.A., Norlander, A.E., integrated TCGA pan-cancer clinical data Joehanes, R., Zhu, J., Chen, B.H., Zhang, B., resource to drive high-quality survival outcome Thorsson, V., Gibbs, D.L., Brown, S.D., Wolf, D., Johnson, A.D., Ying, S., et al. (2015). Integrative analytics. Cell 173, 400–416.e11. Bortone, D.S., Ou Yang, T.-H., Porta-Pardo, E., network analysis reveals molecular mechanisms Gao, G.F., Plaisier, C.L., Eddy, J.A., et al. (2018). of blood pressure regulation. Mol. Syst. Biol. 11, Margolis, K.L., Rodabough, R.J., Thomson, C.A., The immune landscape of cancer. Immunity 48, 799. Lopez, A.M., and McTiernan, A.; Women’s Health 812–830.e14. Initiative Research Group (2007). Prospective Hung, R.J., Ulrich, C.M., Goode, E.L., Brhane, Y., study of leukocyte count as a predictor of incident Torkamani, A., Wineinger, N.E., and Topol, E.J. Muir, K., Chan, A.T., Marchand, L.L., Schildkraut, breast, colorectal, endometrial, and lung cancer (2018). The personal and clinical utility of J., Witte, J.S., Eeles, R., et al. (2015). Cross cancer and mortality in postmenopausal women. Arch. polygenic risk scores. Nat. Rev. Genet. 19, genomic investigation of inflammation pathway Intern. Med. 167, 1837–1844. 581–590. for five common cancers: lung, ovary, prostate, breast, and colorectal cancer. J. Natl. Cancer Inst. Maslah, N., Cassinat, B., Verger, E., Kiladjian, Tsoi, L.C., Stuart, P.E., Tian, C., Gudjonsson, J.E., 107, djv246. J.-J., and Velazquez, L. (2017). The role of LNK/ Das, S., Zawistowski, M., Ellinghaus, E., Barker, SH2B3 genetic alterations in myeloproliferative J.N., Chandran, V., Dand, N., et al. (2017). Large Jabara, H.H., Boyden, S.E., Chou, J., Ramesh, N., neoplasms and other hematological disorders. scale meta-analysis characterizes genetic Massaad, M.J., Benson, H., Bainter, W., Fraulino, Leukemia 31, 1661–1670. architecture for common psoriasis associated D., Rahimov, F., Sieff, C., et al. (2016). A missense variants. Nat. Commun. 8, 15382. mutation in TFRC, encoding transferrin receptor Michailidou, K., Lindström, S., Dennis, J., Beesley, 1, causes combined immunodeficiency. Nat. J., Hui, S., Kar, S., Lemaçon, A., Soucy, P., Glubb, Velazquez, L., Cheng, A.M., Fleming, H.E., Genet. 48, 74–78. D., Rostamianfar, A., et al. (2017). Association Furlonger, C., Vesely, S., Bernstein, A., Paige, analysis identifies 65 new breast cancer risk loci. C.J., and Pawson, T. (2002). Cytokine signaling Nature 551, 92–94. and hematopoietic homeostasis are disrupted in Jiménez-Sánchez, A., Cast, O., and Miller, M.L. Lnk-deficient mice. J. Exp. Med. 195, 1599–1611. (2019). Comprehensive benchmarking and integration of tumor microenvironment cell Mortaz, E., Tabarsi, P., Mansouri, D., Khosravi, A., estimation methods. Cancer Res. 79, 6238–6246. Garssen, J., Velayati, A., and Adcock, I.M. (2016). Vinay, D.S., Ryan, E.P., Pawelec, G., Talib, W.H., Cancers related to immunodeficiencies: update Stagg, J., Elkord, E., Lichtor, T., Decker, W.K., and perspectives. Front. Immunol. 7, 365. Whelan, R.L., Kumara, H.M.C.S., et al. (2015). Jonsson, S., Sveinbjornsson, G., de Lapuente Immune evasion in cancer: mechanistic basis and Portilla, A.L., Swaminathan, B., Plomp, R., therapeutic strategies. Semin. Cancer Biol. 35 Dekkers, G., Ajore, R., Ali, M., Bentlage, A.E.H., Newman, A.M., Steen, C.B., Liu, C.L., Gentles, Suppl, S185–S198. Elmér, E., et al. (2017). Identification of sequence A.J., Chaudhuri, A.A., Scherer, F., Khodadoust, variants influencing immunoglobulin levels. Nat. M.S., Esfahani, M.S., Luca, B.A., Steiner, D., et al. Ward, L.D., and Kellis, M. (2012). HaploReg: a Genet. 49, 1182–1191. (2019). Determining cell type abundance and resource for exploring chromatin states, expression from bulk tissues with digital conservation, and regulatory motif alterations Kasela, S., Kisand, K., Tserel, L., Kaleviste, E., cytometry. Nat. Biotechnol. 37, 773–782. within sets of genetically linked variants. Nucleic Remm, A., Fischer, K., Esko, T., Westra, H.-J., Acids Res. 40, D930–D934. Fairfax, B.P., Makino, S., et al. (2017). Pathogenic Pabst, O., Zweigerdt, R., and Arnold, H.H. (1999). implications for autoimmune mechanisms Targeted disruption of the homeobox Webb, G.J., and Hirschfield, G.M. (2016). Using derived by comparative eQTL analysis of CD4+ transcription factor Nkx2-3 in mice results in GWAS to identify genetic predisposition in versus CD8+ T cells. PLoS Genet. 13, e1006643. postnatal lethality and abnormal development of hepatic autoimmunity. J. Autoimmun. 66, 25–39. small intestine and spleen. Development 126, Klemann, C., Wagner, L., Stephan, M., and von 2215–2225. Zanotti, K.J., Maul, R.W., Castiblanco, D.P., Yang, Hörsten, S. (2016). Cut to the chase: a review of W., Choi, Y.J., Fox, J.T., Myung, K., Saribasak, H., CD26/dipeptidyl peptidase-4’s (DPP4) Ribatti, D. (2017). The concept of immune and Gearhart, P.J. (2015). ATAD5 deficiency entanglement in the immune system. Clin. Exp. surveillance against tumors. The first theories. decreases B cell division and Igh recombination. Immunol. 185, 1–21. Oncotarget 8, 7175–7180. J. Immunol. 194, 35–42. 12 iScience July 24, 101296, 23, 2020 iScience, Volume July 24 Supplemental Information Immune Cell Associations with Cancer Risk Luis Palomero, Ivan Galván-Femenía, Rafael de Cid, Roderic Espín, Daniel R. Barnes, CIMBA, Eline Blommaert, Miguel Gil-Gil, Catalina Falo, Agostina Stradella, Dan Ouchi, Albert Roso-Llorach, Concepció Violan, María Peña- Chilet, Joaquín Dopazo, Ana Isabel Extremera, Mar García-Valero, Carmen Herranz, Francesca Mateo, Elisabetta Mereu, Jonathan Beesley, Georgia Chenevix- Trench, Cecilia Roux, Tak Mak, Joan Brunet, Razq Hakem, Chiara Gorrini, Antonis C. Antoniou, Conxi Lázaro, and Miquel Angel Pujana Supplementary figures legends Fig. S1. Evaluation of immune/stromal cell tissue content estimates in relation to two other methods. Related to Figure 1. (A) Heatmap showing the correlations (Spearman’s ρ) between ConsensusTME- based values and analogous TIMER cell type estimates. (B) Heatmap showing the correlations (Spearman’s ρ) between ConsensusTME- based values and analogous MCP-counter cell type estimates. Fig. S2. Evaluation of immune/stromal cell tissue content estimates in relation to independent leukocyte estimates. Related to Figure 1. Heatmap showing the correlations (Spearman’s ρ) between ConsensusTME-based values and independent leukocyte estimates using the approach of Taylor et al. (2018). Fig. S3. Evaluation of immune/stromal cell tissue content estimates in relation to aneuploidy scores. Related to Figure 1. (A) Heatmap showing the correlations (Spearman’s ρ) between ConsensusTME- based values and aneuploidy scores (Taylor et al., 2018) across major cancer types. (B) Heatmap showing the correlations (Spearman’s ρ) between ConsensusTME- based values and aneuploidy scores (Taylor et al., 2018) across in BRCA subtypes, which show positive correlations in claudin-low. Fig. S4. Differences of inferred immune/stromal cell content between primary tumors with low and high levels of CD274/PDL1 expression. Related to Figure 1. The graphs show the median cell content value in each group and the significance of the difference (Wilcoxon test P value). Fig. S5. Differences of inferred immune/stromal cell content between primary tumors with low and high levels of CD279/PDCD1 expression. Related to Figure 1. The graphs show the median cell content value in each group and the significance of the difference (Wilcoxon test P value). Fig. S6. Correlations between inferred blood immune cell contents and measures from fluorescence-activated cell sorting in blood samples. Related to Figure 1. Forest plot showing correlation estimates and 95% CIs of each inferred cell type (data from whole blood samples of healthy adults; n = 12, GEO GSE127813). Fig. S7. Correlations between inferred immune/stromal cell tissue contents and single cells used to generate pseudo-bulk breast tumors. Related to Figure 1. Each panel shows the correlation between immune cell signature scores (Y-axis) and the number of cells (X-axis) used to generate 100 pseudo-bulk breast tumors (data from Gene Expression Omnibus reference GSE75688). The trend lines, Spearman’s correlations and P values are shown. Fig. S8. Correlations between immune/stromal cell tissue contents and expression of immune benchmark genes. Related to Figure 1. Top panel, distribution of PCCs using data from normal TCGA tissue. Bottom panel, distribution of PCCs using data from primary tumors of TCGA. Mean PCCs and 95% CIs are shown. Fig. S9. Correlations between immune cell signatures and pathway signaling-inferred activities. Related to Figure 1. Unsupervised clustering of the correlation coefficients between inferred cell contents (Y-axis) and KEGG pathway activities (X-axis). Differentiated clusters in normal tissue are marked by red-outlined rectangles. Fig. S10. Gene targets of eQTL recognized in isQTLs are frequently correlated with the corresponding immune/stromal cell signatures. Related to Figure 2. Distributions of random gene sets (same gene set size and equivalent comparisons for each signature and TCGA setting) relative to the number of significant correlations between eQTL-target and immune/stromal signatures. Left- and right-hand panels show results for the first and second isQTL sets presented in the main text, respectively. Empirical test probabilities are shown. Fig. S11. Minimal correlation estimates to detect significant signature- PRS associations. Related to Figure 3. Left and right panels show the lowest correlations required in each normal and primary tumor setting, respectively, to detect nominal (P < 0.05) associations given the TCGA sample sizes. Fig. S12. LUAD and LUSC PRS correlations with NK cell content. Related to Figure 3. Top panels, positive correlations between NK cell content in primary tumors of LUAD and LUSC, and the corresponding PRSs. The adjusted-R2 and P values of the linear regression model are shown. Bottom panels, correlation trends of patients stratified by smoking status, as depicted in the insets. The estimate for LUAD cases classified as current smokers was found to be significantly less than zero (r = -0.12, P = 0.012). Figure S1 A TIMER ConsensusTME Neutrophils Neutrophils Dendritic cells Dendritic cells CD8+ T cells CD8+ T cells CD4+ T cells CD4+ T cells B cells B cells Macrophages Macrophages Spearman’s ρ T A C A BM D IHC C A C IRC D D D M V A GC −O RC CE SC TC PRA US HC TC HNS AA UA OA LC KC 1 −G −L −K GA −P −B −T −U −E −L −L −B −T −C −S 0.5 − − GA GA GA GA GA GA GA GA GA GA GA GA GA GA TC GA GA 0 TC TC TC TC TC TC TC TC TC TC TC TC TC TC -0.5 -1 B MCP-counter ConsensusTME Monocytic lineage Monocytes Fibroblasts Fibroblasts CD8+ T cells CD8+ T cells Cytotoxic lymphocytes Cytotoxic cells B lineage B cells NK cells NK cells Endothelial Endothelial Myeloid dendritic cells Dendritic cells Neutrophils Neutrophils C D D IHC D A CA SC M V AD IRC BM A Spearman’s ρ O US OA A RA RC HC KC N − U A L −G −L −K 1 GA −P −B −L −H −L −P −B −T −C −S GA GA GA GA GA GA 0.5 GA GA GA GA GA TC GA GA TC TC TC TC TC TC TC TC 0 TC TC TC TC TC -0.5 -1 TC GA TC −P AA GA D −S TC KC GA M −U TC CE Figure S2 GA C −E TC SC GA A − TC HNS GA C −L TC I GA HC TC −G GA BM −B TC LC GA A −P TC RA GA D −T TC HC GA A −L TC UAD GA TC −O GA V −C TC OA Leukocyte correlations GA D −T TC GC GA T −B TC RC GA A −L TC U GA SC −K IRC 0 1 -1 B cells 0.5 -0.5 Neutrophils CD4+ T cells CD8+ T cells Macrophages Dendritic cells Spearman’s ρ Figure S3 A Aneuploidy correlations B cells Neutrophils CD8+ T cells CD4+ T cells Dendritic cells Macrophages Spearman’s ρ T IHC D BM C GA CA GA IRC D C D D A M A V C A TC TGC −O TC COA CE AA NS UA RA HC SC US RC 1 KC L −G −L −K GA −P −B −U −H −E −L −L −P −T −B −S 0.5 − GA GA − GA GA GA GA GA GA GA GA GA GA TC GA GA 0 TC TC TC TC TC TC TC TC TC TC TC TC TC TC -0.5 -1 B Aneuploidy correlations B cells CD8+ T cells Dendritic cells Neutrophils CD4+ T cells Macrophages Spearman’s ρ w R2 BC lA lB -lo 1 na HE na TN in 0.5 mi mi ud Lu Lu 0 Cla -0.5 -1 Figure S4 B cells Cytotoxic cells Macrophages Mast cells P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 P < 2.2E-16 0.4 0.25 0.3 0.2 0.25 Score cells Score cells Score cells Score cells 0.00 0.0 0.079 0.0 −0.092 −0.072 −0.122 0.00 −0.249 -0.2 −0.085 −0.248 -0.25 −0.287 -0.3 -0.4 -0.25 -0.50 -0.6 -0.6 Low High Low High Low High Low High PDL1 expression PDL1 expression PDL1 expression PDL1 expression Dendritic cells Eosinophils NK cells Neutrophils P < 2.2E-16 P < 2.2E-16 P < 2.2E-16 0.4 P < 2.2E-16 0.2 0.25 0.25 0.2 Score cells Score cells Score cells Score cells 0.0 0.024 −0.008 −0.126 0.0 0.0 0.0 −0.116 −0.176 -0.2 −0.152 -0.2 −0.255 -0.25 −0.309 -0.25 -0.4 -0.4 -0.50 -0.6 -0.6 Low High Low High Low High Low High PDL1 expression PDL1 expression PDL1 expression PDL1 expression CD4+ T cells CD8+ T cells Macrophages M1 Macrophages M2 0.50 P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 0.50 P < 2.2E-16 0.25 0.25 0.25 0.25 Score cells Score cells Score cells Score cells 0.057 0.081 0.0 −0.042 0.0 −0.053 0.0 0.0 −0.135 −0.099 −0.216 −0.213 -0.25 -0.25 -0.25 -0.25 -0.50 -0.50 -0.50 Low High Low High Low High Low High PDL1 expression PDL1 expression PDL1 expression PDL1 expression T cells gamma-delta T-Reg cells Monocytes Plasma cells P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 P < 2.2E-16 0.2 0.3 0.25 0.25 Score cells Score cells Score cells Score cells 0.0 −0.093 0.0 −0.045 0.0 −0.016 0.0 −0.053 -0.2 −0.182 −0.283 −0.243 −0.167 -0.3 -0.25 -0.25 -0.4 -0.6 -0.6 -0.50 Low High Low High Low High Low High PDL1 expression PDL1 expression PDL1 expression PDL1 expression Endothelial Fibroblasts 0.6 P = 6E-14 P = 3E-12 0.3 0.3 Score cells Score cells −0.028 0.0 −0.068 0.0 −0.037 −0.088 -0.3 -0.3 -0.6 -0.6 Low High Low High PDL1 expression PDL1 expression Figure S5 B cells Cytotoxic cells Macrophages Mast cells P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 P < 2.2E-16 0.4 0.25 0.3 0.2 0.25 Score cells Score cells Score cells Score cells −0.044 0.095 0.00 0.0 0.0 −0.057 −0.091 0.00 -0.2 −0.088 −0.247 -0.25 −0.262 −0.32 -0.3 -0.4 -0.25 -0.50 -0.6 -0.6 Low High Low High Low High Low High PDCD1 expression PDCD1 expression PDCD1 expression PDCD1 expression Dendritic cells Eosinophils NK cells Neutrophils P < 2.2E-16 P < 2.2E-16 P < 2.2E-16 0.4 P < 2.2E-16 0.2 0.25 0.25 0.2 Score cells Score cells Score cells Score cells 0.045 0.01 0.0 −0.085 0.0 0.0 0.0 −0.109 −0.179 -0.2 −0.158 -0.2 -0.25 −0.336 −0.25 -0.25 -0.4 -0.4 -0.50 -0.6 -0.6 Low High Low High Low High Low High PDCD1 expression PDCD1 expression PDCD1 expression PDCD1 expression CD4+ T cells CD8+ T cells Macrophages M1 Macrophages M2 0.50 P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 0.50 P < 2.2E-16 0.25 0.25 0.25 0.25 Score cells Score cells Score cells Score cells 0.001 0.077 0.103 −0.011 0.0 0.0 0.0 0.0 −0.14 −0.104 −0.239 -0.25 −0.239 -0.25 -0.25 -0.25 -0.50 -0.50 -0.50 Low High Low High Low High Low High PDCD1 expression PDCD1 expression PDCD1 expression PDCD1 expression T cells gamma-delta T-Reg cells Monocytes Plasma cells P < 2.2E-16 P < 2.2E-16 0.50 P < 2.2E-16 P < 2.2E-16 0.2 0.3 0.25 0.25 Score cells Score cells Score cells Score cells 0.0 −0.054 0.011 0.003 −0.015 0.0 0.0 0.0 -0.2 −0.172 −0.199 −0.307 −0.277 -0.3 -0.25 -0.25 -0.4 -0.6 -0.6 -0.50 Low High Low High Low High Low High PDCD1 expression PDCD1 expression PDCD1 expression PDCD1 expression Endothelial Fibroblasts 0.6 P < 2.2E-16 P < 2.2E-16 0.3 0.3 Score cells Score cells −0.008 −0.02 0.0 −0.076 0.0 −0.094 -0.3 -0.3 -0.6 -0.6 Low High Low High PDCD1 expression PDCD1 expression Figure S6 0.40 Monocytes 0.80** CD8+ T cells 0.87*** CD4+ T cells 0.91*** Neutrophils 0.64 * NK cells 0.25 B cells -1.0 -0.5 0.0 0.5 1.0 Correlation (r) Figure S7 B cells vs B cells CD4+ T cells vs T cells ρ = 0.88 ρ = 0.84 1.1 CD4+ T cells (signature score) 1.0 P < 2.2E-16 P < 2.2E-16 ● ● 1.0 ● 0.9 ● ● ● B cells (signature score) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.9 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.7 ● ● ● ● ● ● ● ● ● ● ● ● 0.7 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● 0.6 0 5 10 15 20 25 0 10 20 30 B cells (single cells) T cells (single cells) CD8+ T cells vs T cells Macrophages vs Myeloid ρ = 0.81 ρ = 0.90 P < 2.2E-16 Macrophages (signature score) P < 2.2E-16 CD8+ T cells (signature score) 1.0 1.0 ● R = 0.9 , p < 2.2e−16 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.8 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● ● ● ● ● ● ● ● ● ● ● ● 0 10 20 30 0 10 20 T cells (single cells) Myeloid (single cells) Fibloblasts vs Stromal 0.9 ρ = 0.76 P < 2.2E-16 0.8 ● Fibroblasts (signature score) ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.7 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.6 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.5 ● ● ● ● ● ● 0.4 ● ● ● ● ● ● ● ● ● ● 0.3 ● 0 5 10 15 20 Stromal (single cells) Figure S8 Immune cell tissue content correlations with defined immune gene benchmarks Normal tissue 300 Mean = 0.52 250 95% CI: 0.49-0.54 200 Frequency 150 100 50 0 −1.0 −0.5 0.0 0.5 1.0 PCC Primary tumors 150 Mean = 0.60 95% CI: 0.58-0.62 Frequency 100 50 0 −1.0 −0.5 0.0 0.5 1.0 PCC Primary tumors Normal tissue Figure S9 Figure S10 eQTL-gene target correlations with immune/stromal cell signatures isQTL (1st set) isQTL (2nd set) Observed Observed P = 0.017 P = 0.013 Random 0.125 Random 0.20 0.100 0.15 0.075 Density Density 0.10 0.050 0.05 0.025 0.00 0.000 5 10 25 30 35 Significant correlations Significant correlations Figure S11 Minimal correlation value to detect a significant PRS-cell signature associations (P < 0.05) Normal tissue Primary tumors TCGA−PRAD.Effect_overall min(r):0.951 TCGA−SKCM.Effect_overall min(r):0.367 TCGA−BRCA.Effect−ER- min(r):0.311 0.8 TCGA−BLCA.Effect_overall min(r):0.743 0.3 Correlation (r) Correlation (r) TCGA−GBM.Effect_overall min(r):0.295 TCGA−BRCA.Effect−ER- min(r):0.743 TCGA−OV.Effect_OV−overall min(r):0.228 0.6 TCGA−PRAD.Effect_overall min(r):0.229 TCGA−KIRC.Effect_overall min(r):0.219 TCGA−OV.Effect_OV−serous min(r):0.228 TCGA−COAD.Effect_overall min(r):0.206 TCGA−THCA.Effect_overall min(r):0.454 0.2 TCGA−THCA.Effect_overall min(r):0.195 TCGA−LUSC.Effect_overall min(r):0.432 TCGA−BLCA.Effect_overall min(r):0.219 TCGA−LUSC.Effect_overall min(r):0.192 0.4 TCGA−LUAD.Effect_overall min(r):0.427 TCGA−LUAD.Effect_overall min(r):0.167 TCGA−BRCA.Effect−ER+ min(r):0.337 TCGA−BRCA.Effect−ER+ min(r):0.138 TCGA−KIRC.Effect_overall min(r):0.381 TCGA−BRCA.Effect_overall min(r):0.302 TCGA−BRCA.Effect_overall min(r):0.122 20 40 60 80 100 200 300 400 500 Number of samples Number of samples Figure S12 Tumor LUAD Tumor LUSC PRS - NK cells PRS - NK cells 0.2 Tumor stage I/II (normalized score) (normalized score) III/IV 0 0 NK cells NK cells -0.2 -0.2 -0.4 -0.4 P = 0.027 P = 0.045 Adjusted-R2 = 0.04 Adjusted-R2 = 0.02 -2.0 -1.0 -2.0 -1.0 0 PRS PRS Tumor LUAD Tumor LUSC PRS - NK cells PRS - NK cells 0.2 (normalized score) (normalized score) 0 0 NK cells NK cells -0.2 -0.2 -0.4 -0.4 -2.0 -1.0 -2.0 -1.0 0 PRS PRS Smoking status LUAD (n) Smoking status LUSC (n) Current smoker (63) Current smoker (56) Lifelong non−smoker (42) Lifelong non−smoker (4) Current reformed smoker for >15 years (67) Current reformed smoker for >15 years (41) Current reformed smoker for ≤ 15 years (102) Current reformed smoker for ≤ 15 years (105) Transparent Methods TCGA data Clinical and gene expression (RNA-seq fragments per kilobase of transcript per million mapped reads (FPKM) upper quartile normalized (UQ)) data from The Cancer Genome Atlas (TCGA) projects were obtained from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov) and from the corresponding publications. Genetic data at the individual level were obtained following approval by the dbGaP Data Access Committee (project #11689). Metastases and recurrent tumors were excluded from this study, making normal tissue (blood or solid tissue) and primary tumor samples the focus of the analyses. The cancer types are named using the corresponding TCGA study abbreviations (https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/tcga-study- abbreviations). For normal tissue, according to the TCGA protocols, these samples were collected > 2 cm from the tumor margin and/or did not contain tumor identified by histopathological review. The protein expression measures of CD26 and TFCR corresponded to those obtained by TCGA using reverse-phase protein arrays (RPPAs; level 4 data, https://tcpaportal.org/tcpa/). The COAD subtypes were defined based genomic/genetic alterations (chromosomal instability (CIN), genomic stable (GS), and microsatellite instability (MSI) tumors) and on molecular features (consensus molecular subtypes, CMS1-4) (Guinney et al., 2015). 1 Cancer risk variants The variants were compiled from the GWAS Central (Beck et al., 2014) and GWAS Catalog (Buniello et al., 2019) databases, and by literature searches using the PubMed MeSH terms “GWAS”, “association”, “cancer”, and “risk”. The variants are listed in Table S1. The UK Biobank GWAS results were taken from the public repository at http://www.nealelab.is/uk-biobank. Benchmark immune genes These genes were compiled from The Immunological Genome Project (ImmGen) (Shay and Kang, 2013) and CellMarker (Zhang et al., 2019) databases, and by a literature search using the MeSH terms corresponding to the specific immune cell types represented by the gene expression signatures. The benchmarks and their cell type assignments are included in Table S3. Genotype data and imputation Bulk genotyping data corresponding to the Affymetrix Genome-Wide Human SNP 6.0 Array were downloaded from the TCGA legacy archive (https://gdc- portal.nci.nih.gov/legacy-archive/). Of the initial normal tissue and primary tumor samples (n = 16,599), those corresponding to individuals of self-reported non-white origin (n = 4,770), and those of non-European origin based on principal component analysis using variants intersected in the 1000 Genome Project phase III (n = 2,598) were excluded from subsequent analyses; these filters were applied because summary 2 statistics of the GWASs used in this study are strongly biased towards populations of European origin. Normal and tumor samples were then examined separately for duplicates and up to third-degree relatives (kinship cutoff = 0.05), which resulted in the exclusion of an additional 672 samples. In the joint dataset, 765 samples were also excluded because they showed a gender mismatch in an analysis of pseudoautosomal genomic regions. Considering genetic variants, 108 samples that deviated by four or more standard deviations from the mean heterozygosity rate were also excluded. For imputation, variants were excluded if they fulfilled any of the following criteria: they mapped to chromosome Y, pseudoautosomal regions or the mitochondrial genome; they had a call rate < 100%; their minor allele frequency was < 0.01; they departed from Hardy–Weinberg equilibrium (P < 5x10-6); or they mapped to AT-CG sites. Finally, 7,686 samples (4,154 normal, comprising 3,287 blood-derived and 867 solid-tissue samples; and 3,532 primary tumors, of which 94.4% were paired) and 589,101 variants were retained for subsequent analyses. Imputation was performed using the Shape-IT V2 (Delaneau et al., 2008) and IMPUTE2 (Howie et al., 2009) algorithms, and the 1000 Genome Project Phase III panel as reference. Poorly imputed variants (accuracy score < 0.7) were excluded from subsequent analyses. A standard cutoff dose was applied to calculate genotypes using a hard-call threshold of 0.1 (i.e., 0 – 0.1, 0.9 – 1.1, 1.9 – 2.0 for reference homozygote, heterozygous and alternative homozygous genotypes, respectively). 3 Immune/stromal cell signatures Immune/stromal cell gene expression signatures for each TCGA cancer setting were computed using the ConsensusTME method (Jiménez-Sánchez et al., 2019), which was provided available as an R package (https://github.com/cansysbio/ConsensusTME). Ten single-cell breast cancer signatures (Azizi et al., 2018) were included in the TCGA BRCA analyses. Therefore, 18 signatures were examined in each normal tissue and primary tumor setting, except for normal breast and breast cancer tissue, for which a total of 28 signatures were analyzed. The signature scores were computed using the single-sample Gene Set Expression Analysis (ssGSEA) algorithm calculated within the Gene Set Variation Analysis (GSVA) software (Hänzelmann et al., 2013). These scores were calculated for normal tissue and primary tumors, but not for blood samples, since data from blood are limited to germline genotypes. Genes whose expression was uninformative in more than half the samples in a given setting were excluded from the signature calculations; otherwise, missing data were assigned the average value of the informative samples. Evaluation of signature scores computed by two different methods —ssGSEA and summing normalized gene expression Z-scores— revealed global coherence, whereby Pearson correlation coefficients (PCCs) were > 0.80 in 99% (571/578) of the score comparisons. To select independent signatures in each normal and cancer setting, we performed a principal component analysis using the prcomp function in R. Components with eigenvalues > 1 were retained to study quantitative trait loci (subsequent sections). Estimates of immune-related pathway activities were calculated using directed graphs from the Kyoto Encyclopedia of Genes and Genomes 4 (KEGG, https://www.genome.jp/kegg/). Briefly, gene expression profiles were converted into pathway module activity scores by taking into account the chain of reactions from a defined molecular input to a specific molecular output (Cubuk et al., 2018). The 84-gene signature linked to SH2B3 included the genes differentially expressed in Sh2b3-null cells and that participate in genetic and/or protein interactions to this gene/protein (Huan et al., 2015); SH2B3 was excluded from this signature for subsequent analyses. Pseudo-bulk breast tumors To generate 100 pseudo-bulk breast tumors, we used the single-cell RNA-seq data from the Gene Expression Omnibus (GEO) reference GSE75688 (Chung et al., 2017) and aggregated read counts using the aggregateData function in R (https://github.com/HelenaLC/muscat). Each simulated sample of 100 cells was forced to include >50% tumor cells (average 75.3%, 95% CI 72.53 – 77.93%). For non-tumoral cells, 10% of them were fixed as stromal (bulk average 7.22%, 95% CI 6.22 – 8.36%), while the other 90% were a random combination of B cells (average 5.16%, 95% CI 4.28 – 6.28%), T cells (average 6.21%, 95% CI 5.05 – 7.48%), and myeloid cells (average 6.11%, 95% CI 5.05 – 7.39%). Most of the myeloid cells were originally assigned to macrophages (Chung et al., 2017). Quantitative trait loci of immune/stromal cell tissue content The germline genetic calls corresponded to genotype data obtained from blood or normal tissue samples. For cases with both types of sample, the variants with discordant 5 calls were excluded from subsequent analyses. As specified above, the somatic genetic calls corresponded to primary tumors only. The immune/stromal cell-content quantitative trait loci (isQTL) were analyzed using the R/qtl2 package in R (Broman et al., 2019). These analyses included the covariates of gender (when informative), age at diagnosis, tumor stage and histology. The Haley–Knott regression method was used to compute the log odds (LOD) of the associations between genetic variants and immune/stromal cell scores. One thousand permutations were performed in each setting to obtain significance thresholds (Manichaikul et al., 2007) and the variant-signature associations with empirical values of P < 0.05 were considered significant isQTL. The gene targets were defined according to the genomic location of the identified variants. Additional targets were identified by analyzing all variants correlated (r2 > 0.8, 1000 Genomes phase 3, version 5) with each isQTL and intersect them with various functional genomic data, including promoter capture Hi-C (Javierre et al., 2016), annotated enhancers (Hnisz et al., 2013, p.), and eQTL (Schmiedel et al., 2018) from B cells, monocytes, and CD4+ and CD8+ T cells. In addition, correlated variants were queried using the Ensembl Variant Effect Predictor (McLaren et al., 2016) for potential effects on protein coding sequences. Computation of PRSs The PRSs were compiled from the literature and computed by summing the products of the per-allele LOD ratio assigned to each risk variant, and the corresponding allele dosage, for the total number of variants initially defined for each PRS. There was no 6 previous evidence of significant interactions or deviations from a log-additive model in BRCA PRSs (Mavaddat et al., 2019), but it is not known for other cancers. In the analyses of BRCA, OV (no normal tissue data available), and PRAD PRSs, two sets were analyzed, both based on GWAS-identified variants: set #1 (hereafter PRSs-1), which corresponded to scores derived from large collections of GWAS cohorts and validated in independent studies (Mavaddat et al., 2019); and set #2 (hereafter PRSs-2), which corresponded to scores derived from a phenome-wide longitudinal study using electronic health records collected by the Michigan Genomics Initiative (Fritsche et al., 2018). In both sets, PRSs were developed for all BRCA patients, and separately for the estrogen receptor (ER)-positive and ER-negative subtypes. The number of initial variants in these BRCA PRSs and those included in our study, based on available genotypes and obtained imputations were 307 and 185 for PRSs-1, and 3,820 and 3,629 for PRSs-2. As expected, the PRSs from the two sets were found to be positively correlated using germline or primary tumor data: BRCA PRSs PCCs = 0.60 – 0.66, P < 10-5; OV tumors PRSs PCC = 0.72, P < 10-25 (serous PCC = 0.72); and PRAD PRSs PCCs = 0.23 – 0.99, P < 0.01. The Michigan Genomics Initiative also provided PRSs for seven other cancer types, and the number of variants originally included and analyzed in this study were, respectively: 103 and 21 for PRAD; 42 and 41 for COAD; 16 and 16 for BLCA and SKCM; 15 and 15 for OV; 9 and 9 for GBM, LUAD and LUSC; 8 and 7 for THCA; and 7 and 6 for KIRC. 7 Cell signature associations with PRSs The bestNormalize package in R (https://github.com/petersonR/bestNormalize) was used to normalize the cell signature values. The transformation that produced the lowest value from the Pearson’s statistic divided by the degrees of freedom was taken to indicate the best function. The error distributions of the models and Q-Q plots were examined individually. The parameters of each signature transformation are provided in Table S10. Outliers were identified using the interquartile range rule and excluded from subsequent analyses; these were < 5% in all settings. Normalized signature values were used as dependent variables in a linear regression analysis relative to the PRSs. Stepwise analyses including covariates of gender, tumor stage and histology were performed, and the best model was selected based on the Akaike information criterion (AIC). For normal tissue, only those studies with at least 50 informative samples were analyzed. The small number of samples in each setting meant that these analyses could only detect significant (nominal P < 0.05) correlation estimates > 0.27 and > 0.09 in normal breast tissue and BRCA, and stronger correlations would be required in all other settings if nominal significance were to be reached (Fig. S11). The significance of the associations was corrected for multiple testing using the false-discovery rate (FDR) method. Cell signature associations with age at diagnosis The associations between the cell signature scores (dependent variables) and age at diagnosis were evaluated by multiple linear regression, including gender and tumor 8 stage as covariates, the best model being determined from an AIC-based stepwise selection algorithm. The statistical significances of the associations were corrected for multiple testing separately in normal tissue and primary tumor analyses (since the expected effects were the opposite of what they proved to be) using the FDR method. Breast cancer risk in BRCA1/2-mutation carriers Analyses were performed using data from the OncoArray and Collaborative Oncological Gene-environment Study (iCOGS) consortiums with the participation of the Consortium of Investigators of BRCA1/2 Modifiers (CIMBA). The OncoArray and iCOGS designs, quality controls, and statistical analyses have been described previously (Milne et al., 2017). Summary statistics from the retrospective likelihood method are reported. Analysis of blood cell parameters and age at diagnosis of breast cancer Clinical and histopathological data from breast cancer patients were compiled through manual curation of hospital records of the Catalan Institute of Oncology, L’Hospitalet del Llobregat (Barcelona, Catalonia, Spain). Patients were randomly selected from health records collected between 2009 and 2014. The compiled data included date of birth, age, gender (only women selected), date at diagnosis, tumor stage, subtype and/or ER status, and date at initial-diagnostic blood test. The blood test parameters analyzed were the normalized numbers (x109/L) of basophils, eosinophils, leucocytes, lymphocytes, monocytes, neutrophils, and platelets. Linear regressions of each of these parameters 9 on age at diagnosis, including tumor stage and subtype as covariates, were performed. The IDIBELL’s Research Ethics Committee approved this study (reference PR066/20). Supplemental References Azizi, E., Carr, A.J., Plitas, G., Cornish, A.E., Konopacki, C., Prabhakaran, S., Nainys, J., Wu, K., Kiseliovas, V., Setty, M., et al. (2018). Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell 174, 1293-1308.e36. Beck, T., Hastings, R.K., Gollapudi, S., Free, R.C., Brookes, A.J. (2014). GWAS Central: a comprehensive resource for the comparison and interrogation of genome-wide association studies. Eur. J. Hum. Genet. 22, 949–952. Broman, K.W., Gatti, D.M., Simecek, P., Furlotte, N.A., Prins, P., Sen, Ś., Yandell, B.S., Churchill, G.A. (2019). R/qtl2: software for mapping quantitative trait loci with high- dimensional data and multiparent populations. Genetics 211, 495–502. Buniello, A., MacArthur, J.A.L., Cerezo, M., Harris, L.W., Hayhurst, J., Malangone, C., McMahon, A., Morales, J., Mountjoy, E., Sollis, E., et al. (2019). The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res. 47, D1005–D1012. Chung, W., Eum, H.H., Lee, H.-O., Lee, K.-M., Lee, H.-B., Kim, K.-T., Ryu, H.S., et al. (2017). Single-cell RNA-seq enables comprehensive tumour and immune cell profiling in primary breast cancer. Nat. Commun. 8, 15081. Cubuk, C., Hidalgo, M.R., Amadoz, A., Pujana, M.A., Mateo, F., Herranz, C., Carbonell- Caballero, J., Dopazo, J. (2018). Gene expression integration into pathway modules reveals a pan-cancer metabolic landscape. Cancer Res. 78, 6059–6072. Delaneau, O., Coulonges, C., Zagury, J.-F. (2008). Shape-IT: new rapid and accurate algorithm for haplotype inference. BMC Bioinformatics 9, 540. Fritsche, L.G., Gruber, S.B., Wu, Z., Schmidt, E.M., Zawistowski, M., Moser, S.E., Blanc, V.M., Brummett, C.M., Kheterpal, S., Abecasis, G.R., Mukherjee, B. (2018). Association of polygenic risk scores for multiple cancers in a phenome-wide study: results from The Michigan Genomics Initiative. Am. J. Hum. Genet. 102, 1048–1061. Guinney, J., Dienstmann, R., Wang, X., de Reyniès, A., Schlicker, A., Soneson, C., Marisa, L., Roepman, P., Nyamundanda, G., Angelino, P., et al. (2015). The consensus molecular subtypes of colorectal cancer. Nat. Med. 21, 1350–1356. Hänzelmann, S., Castelo, R., Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7. Hnisz, D., Abraham, B.J., Lee, T.I., Lau, A., Saint-André, V., Sigova, A.A., Hoke, H.A., Young, R.A. (2013). Super-enhancers in the control of cell identity and disease. Cell 155, 934–947. Howie, B.N., Donnelly, P., Marchini, J. (2009). A flexible and accurate genotype 10 imputation method for the next generation of genome-wide association studies. PLoS Genet. 5, e1000529. Huan, T., Meng, Q., Saleh, M.A., Norlander, A.E., Joehanes, R., Zhu, J., Chen, B.H., Zhang, B., Johnson, A.D., Ying, S., et al. (2015). Integrative network analysis reveals molecular mechanisms of blood pressure regulation. Mol. Syst. Biol. 11, 799. Javierre, B.M., Burren, O.S., Wilder, S.P., Kreuzhuber, R., Hill, S.M., Sewitz, S., Cairns, J., Wingett, S.W., Várnai, C., Thiecke, M.J., et al. (2016). Lineage-specific genome architecture links enhancers and non-coding disease variants to target gene promoters. Cell 167, 1369-1384.e19. Jiménez-Sánchez, A., Cast, O., Miller, M.L. (2019). Comprehensive benchmarking and integration of tumor microenvironment cell estimation methods. Cancer Res. 79, 6238–6246. Manichaikul, A., Palmer, A.A., Sen, S., Broman, K.W. (2007). Significance thresholds for quantitative trait locus mapping under selective genotyping. Genetics 177, 1963– 1966. Mavaddat, N., Michailidou, K., Dennis, J., Lush, M., Fachal, L., Lee, A., Tyrer, J.P., Chen, T.-H., Wang, Q., Bolla, M.K., et al. (2019). Polygenic risk scores for prediction of breast cancer and breast cancer subtypes. Am. J. Hum. Genet. 104, 21–34. McLaren, W., Gil, L., Hunt, S.E., Riat, H.S., Ritchie, G.R.S., Thormann, A., Flicek, P., Cunningham, F. (2016). The Ensembl Variant Effect Predictor. Genome Biol. 17, 122. Milne, R.L., Kuchenbaecker, K.B., Michailidou, K., Beesley, J., Kar, S., Lindström, S., Hui, S., Lemaçon, A., Soucy, P., Dennis, J., et al. (2017). Identification of ten variants associated with risk of estrogen-receptor-negative breast cancer. Nat. Genet. 49, 1767–1778 Schmiedel, B.J., Singh, D., Madrigal, A., Valdovino-Gonzalez, A.G., White, B.M., Zapardiel-Gonzalo, J., Ha, B., Altay, G., Greenbaum, J.A., McVicker, G., et al. (2018). Impact of genetic polymorphisms on human immune cell gene expression. Cell 175, 1701-1715.e16. Shay, T., Kang, J. (2013). Immunological Genome Project and systems immunology. Trends Immunol. 34, 602–609. Zhang, X., Lan, Y., Xu, J., Quan, F., Zhao, E., Deng, C., Luo, T., Xu, L., Liao, G., Yan, M., et al. (2019). CellMarker: a manually curated resource of cell markers in human and mouse. Nucleic Acids Res. 47, D721–D728. 11