The prognostic landscape of genes and infiltrating immune cells across human cancers

A cancer-wide atlas of prognostic genes

We assembled, curated, and integrated cancer gene expression and clinical outcome data from the public domain into a new resource for PREdiction of Clinical Outcomes from Genomic profiles (PRECOG, Fig. 1a, Supplementary Table 1, available at http://precog.stanford.edu). PRECOG encompasses nearly 30,000 expression profiles from 166 cancer expression datasets covering 39 distinct malignant histologies (Fig. 1b). Importantly, all data were manually curated with respect to relevant clinical parameters, including disease status, stage, and histology, and key results from selected studies (e.g., survival plots) were verified (Methods). To avoid ambiguities in outcome annotation, we restricted the results presented herein to the ~18,000 cases with overall survival data. We further confirmed consistency and quality of integrated datasets and curated annotations by assessing concordance between molecularly-inferred versus clinically-reported patient gender. In all, 98% of tumors were gender-concordant, with less than 2% of tumors affected by likely annotation errors (such as “off-by-one” mismatches) (Supplementary Fig. 1a–c). All microarray studies in PRECOG were consistently normalized and pre-processed (Methods) and associations between gene expression and survival were assessed by univariate Cox regression.

To compare prognostic associations across independent datasets, and to minimize the confounding influence of batch effects, the statistical associations between genes and clinical outcomes were assessed by z-scores. Z-scores are directly related to P values, but conveniently encode both the directionality and robustness of statistical associations (Supplementary Fig. 1d). Moreover, z-scores have a straightforward interpretation; they represent the number of standard deviations from the mean of a normal distribution. For example, |z| > 1.96 is equivalent to a two-sided P < 0.05 (Supplementary Fig. 1d). Unlike parameters such as hazard ratios, z-scores are independent of different time-scales measuring survival follow-up times, and of the range/scale of predictor variables, permitting direct comparison across studies and platforms. To facilitate cross-cancer analyses, z-scores for individual studies were combined to yield “meta-z scores” for the prognostic significance of each gene in each cancer type (Methods; Supplementary Table 1). We observed high concordance between meta-z scores and z-scores, where the latter were obtained by first merging expression data from multiple studies of the same cancer (e.g., lung adenocarcinoma, Spearman's R = 0.9, P < 2.2×10⁻¹⁶; Methods). To further evaluate the robustness of the meta-z metric, we calculated a global meta-z score for each gene across all cancers, and compared PRECOG to a validation set of 9 independent studies that were held-out (Supplementary Table 1). Globally prognostic genes were significantly correlated between PRECOG and the validation set (R = 0.55, P < 2.2 × 10⁻¹⁶; Supplementary Fig. 2a,b). In addition, pan-cancer prognostic genes were significantly concordant between PRECOG and another validation set comprised of studies profiled by RNA-seq from The Cancer Genome Atlas (TCGA) (R = 0.52, P < 2.2 × 10⁻¹⁶; Supplementary Fig. 2a,b). We also evaluated the influence of batch effects²¹ on z-score values. Notably, only modest differences in z-scores were observed following batch effect removal (e.g., for samples run on different dates) (Supplementary Fig. 2c–e).

Pan-cancer prognostic genes

PRECOG provides an unprecedented opportunity to quantify commonalities in prognostic genes across a large number of human malignancies. We found that prognostic genes (filtered at |meta-z| > 3.09, or nominal one-sided P < 0.001) are significantly more likely to be shared by distinct tumor types than expected by random chance (Fig. 1c, Supplementary Table 2). This result was reproducible across a broad range of statistical thresholds (Supplementary Fig. 3a,b), and is reminiscent of the high cancer-wide concordance reported among somatic aberrations influencing genome-wide copy number²². Conversely, cancer-specific prognostic genes are less frequent than expected by random chance (Fig. 1c, Supplementary Fig. 3a,b), and predominantly reflect tissues of origin (Supplementary Fig. 3c, Supplementary Table 2).

To obtain a global map of prognostic patterns, we clustered survival-associated z-scores across all 166 PRECOG datasets using AutoSOME, an unsupervised method that is robust against outliers and does not require pre-specification of the number of clusters²³ (Fig. 1d, Supplementary Table 3). Prognostic clusters include genes involved in cell adhesion and epithelial-mesenchymal transitions, vascularization, and immunological and proliferative processes (Supplementary Table 3). When clusters were ordered by a metric that integrates gene-level meta-z scores and cluster size, the two largest clusters were most highly ranked (Fig. 1d, left; Methods). One of these two clusters is broadly associated with inferior outcomes, and is functionally linked to cell proliferation and cell cycle phase (Fig. 1d, right). While this cluster is prognostic in many solid tumors, such as breast and lung adenocarcinoma, proliferation genes were not adversely prognostic in some cancers, including colon cancer and AML (Supplementary Table 1), two malignancies for which the clinical relevance of generally quiescent cancer stem cells has been demonstrated^24,25. The other large cluster is associated with favorable survival and is highly enriched in immunological processes and immune response genes (Fig. 1d, right; Supplementary Table 3), suggesting that the immune microenvironment is a key factor contributing to favorable survival across cancers.

To further explore cancer-wide prognostic signatures, we used PRECOG to define robust pan-cancer survival models. First, we determined the number of histologies needed to identify genes with maximal prognostic power. Using a cross-validation approach to avoid outliers, we observed quantitative improvements in the significance of pan-cancer prognostic genes until ~30 distinct histologies were sampled, after which additional gains were marginal (Fig. 2a, left). With this framework, we identified the top ten adverse and top ten favorable prognostic genes, largely consisting of genes regulating or oscillating with the mitotic cell cycle and markers of distinct lymphocyte subsets (Fig. 2a, right). Expression of the proto-oncogene FOXM1 (encoding forkhead box M1)²⁶ was most frequently associated with adverse risk, andoutperformed MKI67, which encodes a protein (Ki-67) used clinically as a marker of proliferation²⁷ (Supplementary Fig. 4a). Expression of KLRB1 (encoding CD161, a surface marker on several T cell subsets²⁸ and NK cells; Supplementary Fig. 4b) was most frequently associated with favorable outcomes. Notably, FOXM1 and KLRB1 were among the top pan-cancer genes in validation sets, including TCGA RNA-seq data (Fig. 2b and Supplementary Fig. 2). Bivariate models based on these two genes outperformed either gene alone (Supplementary Table 4). Moreover, when weighted by their Cox regression coefficients in a training set (Supplementary Table 4), a FOXM1-KLRB1 composite score stratified patient survival in diverse internal validation datasets, including carcinomas, brain tumors, and hematopoietic neoplasms (Fig. 2c), and was prognostic in multivariate models integrating common clinical indices (Supplementary Table 4). Furthermore, the FOXM1-KLRB1 composite score remained significant in validation datasets held-out from PRECOG, including TCGA RNA-seq datasets (Supplementary Fig. 4c,d).

We next asked whether integrating cancers in a meta-analysis might illuminate additional functional elements of biological programs impacting survival. To address this, we evaluated connectivity among top prognostic genes within protein-protein association networks. For adversely prognostic genes, a considerably higher connectivity was achieved when considering all cancer histologies together as opposed to separately (Fig. 2d, bottom; Supplementary Table 5), suggesting that individual malignancies in PRECOG contribute complementary information to a pervasive proliferation program (Fig. 2d, middle). We further characterized this network by integrating PRECOG with data from ENCODE²⁰, ChEA²⁹, and mSigDB³⁰ for 1,006 gene sets defined by shared transcription factor binding sites. Adversely prognostic genes in this network and across PRECOG were most significantly enriched in FOXM1 ChIP-Seq binding targets³¹ (Supplementary Fig. 4e, Supplementary Table 5), which together with FOXM1 itself, may be drivers of inferior survival. In contrast, we observed no difference in protein-protein connectivity among favorably prognostic genes, whether identified from pooled cancers or individual cancers (Fig. 2d, bottom; Supplementary Table 5). We hypothesized that transcriptome heterogeneity among distinct TALs, which account for the majority of favorable pan-cancer prognostic associations, may not be well captured by aggregated protein-protein association studies. In addition, several highly prognostic genes, including KLRB1, are expressed by multiple immune subsets (Supplementary Fig. 4b). Therefore we performed an “in silico dissection” of PRECOG to infer the specific leukocyte subsets associated with survival across cancers.

Leukocyte composition in bulk tumors

Infiltration of tumors by specific leukocyte cell subsets such as CD8⁺ and CD45RO⁺ memory T-lymphocytes has been largely linked with favorable outcomes in different cancers²,³², while others such as regulatory T-cells and macrophages can confer good or poor prognosis depending on context^33-36. To systematically and comprehensively map compositional differences in TALs and their relationships to survival, we applied a new machine-learning framework for Cell-type Identification By Estimating Relative Subsets Of known RNA Transcripts (or CIBERSORT)¹⁶. CIBERSORT outperforms previous deconvolution methods with respect to noise, unknown mixture content, and closely related cell types, in statistically estimating relative proportions of cell subsets from expression profiles of complex tissues (e.g., bulk tumors)¹⁶. As input, we used purified expression profiles for 22 distinct leukocyte subsets, and defined “barcodes” of gene expression signatures that robustly distinguish these cell types without requiring cell type-specific marker genes¹⁶. At a |meta-z score| > 3.3 (corresponding to two-sided P < 0.001), 28% of these barcode genes (152 of 547) are individually significant in PRECOG, out of 2,851 total pan-cancer prognostic genes at the same significance threshold. This is higher than expected by random chance (P < 0.001, Chi-squared test). Whether directly or indirectly compared against flow cytometry and immunohistochemistry, CIBERSORT exhibited robust performance on solid tumors, accurately estimating relative fractions of leukocyte subsets in colorectal cancer and lung adenocarcinoma (Fig. 3a), and follicular lymphoma¹⁶.

Applied to PRECOG, CIBERSORT revealed striking differences in relative leukocyte composition between hematopoietic neoplasms, brain cancers, and non-brain solid tumors (Fig. 3b, Supplementary Table 6). Variation in TAL content was also consistent and reproducible across independent studies of the same cancer type, including solid tumors (Supplementary Fig. 5a). Of note, while the majority of tumors profiled within PRECOG were unpurified and uncontrolled with respect to tumor content (Supplementary Table 1), CIBERSORT correctly inferred high fractions of plasma cells in multiple myeloma-enriched specimens (Fig. 3b). Furthermore, as expected, B-cell signatures were found to predominate in B-cell malignancies (Fig. 3b), suggesting that CIBERSORT has general utility for discerning cell of origin in diverse cancers.

Prognostic associations of TALs

To complement our gene-centric survival analysis, we assembled a global map of prognostic associations for 22 immune populations across human malignancies (Supplementary Fig. 6a). We observed considerable variation between cell subsets and cancer-specific outcomes, and many of these associations are statistically significant (Supplementary Fig. 6b–d). Pooling cancers yielded significant global leukocyte prognostic patterns, in which higher levels of estimated T-cell fractions were found to generally correlate with superior survival while increasing levels of myeloid populations primarily correlated with poorer survival. Intra-tumoral γδ T-cell³⁷,³⁸ and polymorphonuclear (PMN)³⁹,⁴⁰ signatures emerged as the most significant favorable and adverse cancer-wide prognostic populations, respectively (Fig. 3c, left). Moreover, when inferred leukocyte fractions were compared with KLRB1 expression across cancers, γδ T-cell and CD8 T-cell signatures were most highly correlated (Supplementary Fig. 5b), suggesting a link to the prognostic significance of this gene. We found no relationship between estimated PMN levels in datasets with annotated necrotic tissue content (Methods), suggesting that intra-tumoral PMNs are not simply a correlate of tissue necrosis⁴¹. Furthermore, consistent with previous reports³³,⁴², signatures of tumor-associated M2 macrophages were found to predict worse outcomes than pro-inflammatory M1 macrophages, and anti-CD3/anti-CD28-costimulated, but not resting, CD45RO⁺ memory helper T-cells were correlated with superior outcomes.

Prognostic TALs in solid tumors

By comparing leukocyte survival signatures in breast and lung cancer—two of the most highly profiled cancers in PRECOG—we identified two populations, PMNs and plasma cells (PCs), with unexpectedly strong yet reciprocal relationships to survival (Fig. 3d). PC signatures are significant predictors of favorable survival across solid tumors in general (Fig. 3c, right), and were the most inversely correlated prognostic population to PMNs (Fig. 4a) when assessed globally in a cross-correlation analysis between human cancers (Supplementary Fig. 5c). Estimated PC levels were not correlated with tumor stage (Supplementary Fig. 7a). Since PC signatures were found to be higher in tumors than in adjacent normal tissues (Supplementary Fig. 7b), the prognostic value of tumor infiltrating PCs is unlikely a proxy for general immunological health, supporting a role for antigen-driven processes required for their clonal expansion and emergent humoral immune responses⁴³. Furthermore, a simple ratio of estimated PMN to PC levels was found to be significantly prognostic in diverse solid tumors (Fig. 4b).

To experimentally evaluate the reciprocal survival associations of PMN and PC signatures, we assessed their infiltration of 187 lung adenocarcinomas using tissue microarray (TMA) analysis (Supplementary Table 7). Characteristics of both cell types were observed by H&E staining of tissue sections (Supplementary Fig. 7c,d), and the presence of tumor infiltrating plasmacytic cells (i.e., plasmablasts or plasma cells) was confirmed in fresh tumor specimens using both flow cytometry (Supplementary Fig. 7e) and morphological assessment (Supplementary Fig. 7f). Moreover, we confirmed an elevated presence of plasmacytic cells in non-small cell lung cancer (NSCLC) tumors, as compared to normal adjacent tissues (Supplementary Fig. 7g,h). In serial lung adenocarcinoma tissue sections, we stained for the presence of MPO (myeloperoxidase) and IGKC (Immunoglobulin kappa constant), markers of PMNs and PCs, respectively (Supplementary Fig. 8a). Since B-cells express varying levels of IGKC, we also tested for CD20, a surface marker of mature B-cells but not PCs (Supplementary Fig. 7e). We found <10% overlap with CD20, indicating the high specificity of IGKC for PCs (Supplementary Fig. 8b; Methods). Next, we quantitated the staining area for each marker in the tissue array (Methods; Supplementary Fig. 8c,d). While operating on differing scales and measured on independent tumor specimens, fractional levels of these three markers measured in situ on TMAs were comparable to relative infiltrate levels inferred by CIBERSORT (Fig. 4c). Moreover, in both continuous and binary models, we found a strong relationship between inferior survival and a higher ratio of PMN to PC levels in lung adenocarcinoma, whether measured in PRECOG (Fig. 4d), in held-out microarray validation datasets (Supplementary Fig. 8e), or by surrogate markers in tissue microarray specimens (Fig. 4e). Furthermore, TMA results remained significant in multivariate models incorporating relevant clinical parameters (Supplementary Table 7). Together, these data validate our computational approach, and demonstrate that tumor-associated PMNs and PCs exhibit opposite associations with overall survival.

Circulating leukocytes, including PMNs and B-lymphocytes contribute to the tumor microenvironment^44-46, and leukocyte frequencies of innate and adaptive effectors in peripheral blood can have prognostic value⁴³,⁴⁷. Therefore, we examined a subset of NSCLC patients from the TMA with available peri-operative complete blood counts to assess the concordance between levels of circulating leukocytes and TALs. While intra-tumoral PMN to PC ratios remained significantly prognostic within this subset, we found no significant correlation between circulating and infiltrating compartments, and no prognostic value from circulating leukocyte levels (Supplementary Table 7).

PRECOG assembly and quality control

To identify cancer gene expression datasets with corresponding patient outcome data we queried NCBI Gene Expression Omnibus (GEO), EBI ArrayExpress, NCI caArray, and Stanford Microarray Database for the terms, survival, prognosis, prognostic, or outcome. Perl scripts were implemented to download processed and raw data, and associated annotation. For data within NCBI, the array platform was determined from the SOFT format file, and the corresponding annotation file was retrieved from GEO. From these, the Probe ID, Genbank accession, HUGO gene symbol and gene description were extracted based on the internal headers of the SOFT annotation file. The desired fields were specified manually if this automated procedure failed. For older platforms, such as cDNA microarrays, where annotations had not been recently updated, we re-mapped the probe sequences to HUGO gene symbols via the Genbank or Refseq accession number through the NCBI Entrez gene identifier. In cases without available accessions, but with the DNA sequence of the probe, we performed the mapping using BLAT to compare probes to a Refseq reference and look for unique highest-scoring hits.

Scripts were written to extract sample annotation information from GEO SOFT format files and parse them into tables. Since the contents of annotation fields are not semantically enforced, sample data can be contained in various fields, including Sample_title, Sample_characteristics, Sample_description, and Sample_source. Moreover, not all fields are specified for every sample. To parse this information into tabular format, we attempted to estimate the correct variable name (column header) by searching for common substrings across samples. In some cases, a dataset clearly had survival information, but was not deposited with the genomic data. In such cases, we first searched supplementary information of corresponding literature for the missing information. Failing this, we contacted corresponding and first authors, of which roughly half supplied the requested data.

All tabulations of clinical annotations were further checked and manually curated. This process included verification of results in selected studies by direct comparison of Kaplan-Meier plots and time scales with those in the corresponding primary publications, as well as consistency of prognostic genes across studies. Separately, errors due to technical issues or the curation process were estimated by comparing annotated gender to the ratio of RPS4Y1 to XIST (male:female) expression levels⁵⁶ after microarray normalization, as detailed below (Supplementary Fig. 1a–c). Furthermore, identical samples present in more than one dataset were identified using MD5 checksums for Affymetrix data, and by cross-correlation analysis of expression vectors, and redundant samples were accordingly eliminated.

We applied the following gene expression normalization strategy to allow unification of data from diverse microarray platforms within PRECOG. For Affymetrix GeneChip data, raw CEL files were obtained when possible, and were normalized with the MAS5 algorithm (affy package v. 1.26 of Bioconductor v. 1.8 in R 2.15.1), using a custom CDF (Chip Definition File) for probeset summarization, which updates and maps array oligonucleotides to Entrez gene identifiers^57-59 (http://brainarray.mbni.med.umich.edu/Brainarray/). Each dataset, regardless of platform, was quantile normalized separately. Moreover, each gene was log₂ transformed if not already in log space, and was then unit mean/variance standardized across samples within a given dataset. While alternative microarray normalization methods have been proposed (e.g., RMA⁶⁰, gcRMA⁶¹, fRMA⁶², SCAN-UPC⁶³), for survival analysis we did not observe any significant benefit in comparing Affymetrix data normalized as described above to alternate normalization strategies (data not shown). TCGA RNA-seq and clinical data were downloaded from the TCGA Data Coordinating Center using TCGA-assembler⁶⁴. The gene-level RNA-seq data were pre-processed using TCGA-assembler's ProcessRNASeqData function. RNA-seq and clinical data were matched via the patient barcode provided by TCGA.

For each study, the association of each probe on an array platform with survival outcomes was assessed via Cox proportional hazards regression using the coxph function of the R survival package (v. 2.37). Cox coefficients, hazard ratios with 95% confidence intervals, P values, and z-scores were obtained for each array probe. For datasets that had not been processed with Custom CDF, which yields a unique per-gene expression value, survival z-scores for probes were collapsed to the gene level by averaging z-scores of probes that matched to the same HUGO gene symbol. Z-scores for each gene were summarized across all datasets in each malignancy using Lipták's weighted meta-z test⁶⁵,⁶⁶, with weights set to the square roots of sample sizes⁶⁷. To identify genes with cancer-wide prognostic significance, and avoid bias due to cancers with different sample sizes, we further combined weighted meta-z-scores into a single global meta-z-score for each gene using Stouffer's method (unweighted)⁶⁶.

Validation of z statistics in PRECOG

Using lung adenocarcinoma as a test case, we assessed the relationship between the weighted meta-z-score metric and standard z-scores, the latter of which were derived from a merged expression matrix consisting of GEPs from lung adenocarcinoma studies in PRECOG (Supplementary Table 1). For this purpose, we selected datasets that had at least 40 stage I samples (indicated in Supplementary Table 1). To mitigate batch effects, we standardized each gene in each dataset such that it had unit mean and variance across stage I samples. Sample annotations were manually reviewed to ensure that staging corresponded to American Joint Committee on Cancer (AJCC) version 6 (2002), based on TNM (Tumor-Nodes-Metastasis) information. Many datasets pre-dated version 7 of AJCC, and did not contain the required detail for annotating to that standard. These refinements and standardizations permitted merging of samples from different datasets comprising different array platforms and different distributions of tumor stage across the cohort. In all, we compared lung adenocarcinoma GEPs from n = 1,106 patients, and found that weighted meta-z scores are significantly correlated with merged z-scores (Spearman's R = 0.9, P < 2.2 × 10⁻¹⁶). We observed similar results when comparing the meta- and merged-z statistics for a compendium of 5 AML studies, thus validating our use of the meta-z statistic. Of note, while we applied batch-correction procedures to merge expression datasets prior to calculating cross-study z-scores, these steps were not necessary with the meta-z metric, as z-scores from individual studies were directly integrated. This suggests that the meta-z approach effectively overcomes batch differences across datasets.

We further evaluated the influence of batch effects²¹ within individual datasets using Combat⁶⁸ (Supplementary Fig. 2c–e). Applied to microarray processing dates in four AML studies, we observed only a modest effect on prognostic z-scores, as pre- and post-batch-corrected data were all highly correlated (R ≥ 0.92, P < 2.2 × 10⁻¹⁶; Supplementary Fig. 2c). To test whether batch correction of samples profiled by different study sites would improve data quality, we compared pre- and post-batch-corrected expression data from the NCI director's challenge lung adenocarcinoma dataset (ca00182) with a control dataset consisting of prognostic meta-z scores from a pooled set of all remaining 19 lung adenocarcinoma studies in PRECOG. Little difference in performance was observed for the most prognostic genes, with changes primarily affecting genes whose association with survival outcomes was subtle. (Supplementary Fig. 2d,e).

PRECOG false discovery rate

While z-scores and meta-z scores were analyzed in this work, Q values for global unweighted meta-z and weighted cancer-specific meta-z-scores were estimated using the False Discovery Rate (FDR) method of Storey and Tibshirani⁶⁹, and are available for all analyzed z-score matrices online (http://precog.stanford.edu). Notably, of 23,288 HUGO gene symbols in PRECOG, 4,385 (19%) have a global meta-z significant at Q < 0.05 (|meta-z| > 2.6), and 2,986 (13%) are significant at a Q < 0.01 (|meta-z| > 3.22) (Supplementary Table 1).

Blinding and sample selection criteria

No blinding was used in this work. Duplicate and non-diagnostic (relapse) samples were excluded from analysis.

Statistical analysis of shared prognostic genes

To determine an empirical P value for the fraction of overlapping genes in Figure 1c, we randomized gene labels for every cancer, and for each cancer, calculated the fraction of prognostic genes shared with at least 1 other cancer (again, using |meta-z| > 3.09). We then calculated the median fraction of shared genes across all cancers and repeated the analysis 100,000 times. We performed a similar analysis in Supplementary Fig. 3a,b, but used 1,000 Monte Carlo iterations to test a broad range of P value and Q value thresholds, considering prognostic genes shared by at least 2, 3, or 4 tumor types in PRECOG.

Clustering of the PRECOG z-score matrix

To identify groups of genes with similar prognostic patterns across the entire PRECOG database (166 cancer datasets and 22,461 gene symbols; 827 genes represented in only a small subset of cancers were excluded), we performed unsupervised cluster analysis using AutoSOME²³. Prior to clustering, columns (i.e., cancer datasets) were pre-processed to unit variance, rows (i.e., genes) were median-centered, and sum of squares = 1 normalization was applied to rows and then columns²³. AutoSOME was run with a P value threshold of 0.01, 100 ensemble iterations, and otherwise default settings. In all, 665 clusters were identified, and ~50% genes were assigned to the 55 largest clusters (Supplementary Table 3). Functional annotation for each cluster was assessed using a variety of published gene sets and significance of overlap was determined by a hypergeometric test (for PubMed IDs, see Supplementary Table 3). Gene sets with a Bonferroni-corrected P < 0.01 were considered significant. To further evaluate each cluster, we calculated a global meta-z score for each gene (as in Supplementary Table 1), but considered only absolute meta-z scores to avoid biasing clusters with pan-cancer prognostic genes with the same survival orientation. We then calculated a compound score for each cluster by applying the same formula used for the unweighted meta-z score, which weights cluster-wide meta-z scores by cluster size. Importantly, this formula is not intended to yield a prognostic score in this context, since intra-cluster meta-z scores are not independent variables. Moreover, our approach is not simply a proxy for cluster size, as the top 5 ranking clusters were originally ranked 2, 1, 4,13, and 7 with respect to size. The top 5 clusters ranked by decreasing compound scores are shown in Figure 1d.

Cross-validation of top prognostic genes

The utility of PRECOG as a cancer-wide meta-analysis tool was assessed by cross-validation. Cancer histologies (from Supplementary Table 1) were randomly selected and pooled in a training set, with the number of cancers t ranging from 1 to 31 (of 39) total cancers. For each t, the top 10 adverse and favorable prognostic genes were identified using the mean meta-z score across the training set, and the mean meta-z for the same genes was determined in the remaining (39 – t) cancers (validation set). This process was repeated 100 times for each t, and the results are plotted in Fig. 2a. To determine the top prognostic genes across PRECOG, the top 10 adverse and favorable prognostic genes in each validation set were recorded at t = 31, and the most frequently occurring genes are shown in Fig. 2a.

Construction and assessment of a FOXM1-KLRB1 prognostic model

In brief, we used the coxph function in the R survival package to test the prognostic significance of FOXM1 and KLRB1 in bivariate models across PRECOG (each dataset was normalized as described for PRECOG construction). Bivariate models with a P < 0.05 (Wald test) were considered significant. Derivation of the FOXM1-KLRB1 composite score is described in detail below.

To integrate FOXM1 and KLRB1 into a composite score, we analyzed cancer types with at least two independent datasets having a FOXM1, KLRB1 bivariate P < 0.05. Of these, we extracted bivariate coefficients from a randomly selected group of 20% of the corresponding datasets, which served as the training set (n = 8; also see Supplementary Table 5). The median value of each coefficient in the training set was used as an independent weight in the following formula to define a composite score, (FOXM1 × 0.243) + (KLRB1 × −0.169). The median value of the FOXM1-KLRB1 composite score was determined for each left-out dataset (validation set), and used to stratify corresponding patients into high and low risk groups, which were then aggregated by cancer type, and subjected to Kaplan Meier analysis (Fig. 2b). Importantly, each cancer type in Fig. 2b only includes left-out datasets (i.e., validation data) for which FOXM1 and KLRB1 expression values were both available (see Supplementary Table 1 for datasets used). Moreover, the composite score was not re-optimized for RNA-seq when applied to TCGA data.

PRECOG network connectivity analysis

We evaluated the functional coherence of top prognostic genes in PRECOG by interrogation of human protein-protein associations in STRING version 9.0⁷⁰. STRING integrates evidence from multiple sources, including curated databases, experimentally confirmed physical interactions, co-expression data, and associations inferred via text mining⁷⁰. We compared the connectivity of STRING networks (for edges with a combined confidence level of at least 0.4) for the top 100 adverse and favorable prognostic genes identified from all cancers (by mean meta-z score), individual cancers (meta-z score), and individual datasets (z-score). For each set of proteins, the largest connected component was assessed by two metrics, (1) the average number of associations (i.e., edges) per protein (i.e., node), and (2) algebraic connectivity, a graph theoretic measure of overall network connectedness⁷¹. Only networks with at least 10 nodes were considered. Functional enrichment of global networks shown in Fig. 2c was performed with ToppFun⁷² and P values were determined using a Benjamini-Hochberg-corrected hypergeometric test.

Geneset analysis of PRECOG clusters and ranked z-scores

Enrichment of biological processes with respect to survival z-scores was assessed in two ways. First, cluster memberships defined by AutoSOME were compared to pre-defined gene sets by hypergeometric test. Comparison gene sets were obtained from the Molecular Signatures Database (mSigDB)⁷³, and additional sets of genes that are targets of specific transcription factors were extracted from CHEA²⁹ and ENCODE²⁰. Ranked lists of gene survival z-scores were also analyzed with respect to these gene sets using the PreRanked tool of Gene Set Enrichment Analysis⁷⁴.

Inferring TAL levels in bulk tumor GEPs

The samples profiled within PRECOG primarily represent bulk diagnostic pre-therapy tumor specimens, which often contain a variety of cell types, including diverse TALs. Given the enrichment of lymphocyte markers in favorably prognostic genes across PRECOG (Figs. 1d, 2c), a method to systematically “unmix” or deconvolve bulk tumor GEPs in PRECOG may reveal new insights into tumor immunobiology. We recently developed a new approach for Cell Type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT), a machine-learning method that outperformed other approaches in benchmarking experiments¹⁶. CIBERSORT produces an empirical P value for the deconvolution using Monte Carlo sampling. Like other linear deconvolution methods, CIBERSORT only operates on expression values in non-log linear space⁷⁵.

TAL heterogeneity and prognostic associations

CIBERSORT was applied to all normalized PRECOG GEPs from Affymetrix HGU133 platforms (57 studies and 25 cancers; Supplementary Table 6). In all, 5,782 tumor GEPs were successfully deconvolved (CIBERSORT P < 0.005). For each dataset, estimated mRNA fractions of each leukocyte subset were related to survival using univariate Cox regression. Weighted meta-z scores were determined using the same approach described for PRECOG in order to build an immune-centric version of PRECOG (iPRECOG, Supplementary Fig. 6a), and unweighted global meta-z scores were used to summarize pan-cancer leukocyte associations in Figure 3c.

Immune-PRECOG false discovery rate

To differentiate real from stochastic variation in inferred leukocyte prognostic associations, we first compared P values and meta-z scores in immune-PRECOG (Supplementary Fig. 6b), as any deviation from a standard normal distribution must be considered when drawing statistical conclusions. We generated 1000 null meta-z matrices by (1) shuffling the cell type fractions inferred for each dataset in Supplementary Table 6, and (2) computing z-scores and corresponding meta-z scores to capture relationships to overall survival. We found a tight correspondence between the distribution of null meta-z scores and a standard normal distribution (Supplementary Fig. 6b). Having validated the normality of the meta-z score, we then filtered Supplementary Fig. 6a using a range of statistical significance thresholds, and at each cutoff, compared observed versus expected fractions for all leukocyte prognostic associations (Supplementary Fig. 6c). At a two-sided P value threshold of 0.05 (|z| > 1.96), we found nearly three times more prognostic associations than would be expected by random chance; at P < 0.01, there is a five-fold enrichment, which continues to increase with lower P value cutoffs (Supplementary Fig. 6c).

Separately, we performed a similar analysis on the global meta-z scores shown in Fig. 3c. Here, we integrated the null meta-z scores from Supplementary Fig. 6c into null global meta-z scores and recomputed the analysis shown for pan-cancer leukocyte prognostic associations (plotted as the fraction of leukocyte subsets retained at different significance thresholds; Supplementary Fig. 6d). Taken together, these results explicitly quantify significant versus stochastic variation in leukocyte prognostic associations at different statistical cutoffs, and allow others to tune the nominal statistical threshold to achieve a desired false discovery rate.

Relative PMN levels versus necrotic tissue content

Relative RNA fractions of PMNs inferred by CIBERSORT were not correlated with annotated necrotic content in lung squamous cell carcinoma (TCGA; R² = 0.01; P = NS) or melanoma (microarray dataset GSE8401⁷⁶; R² ~ 0; P = NS).

Flow cytometry versus CIBERSORT

Flow cytometry analysis of non-small cell lung cancer tumor (n = 13) specimens was performed as described below, and median fractions of CD4⁺, CD8⁺, CD19⁺, CD56⁺, and CD14⁺ populations were normalized by overall CD45⁺ content (Fig. 3a). For comparison with CIBERSORT, leukocyte signature matrix populations were grouped into the same cluster of differentiation categories: CD14⁺, monocytes, macrophages, and dendritic cells; CD4⁺, all T cell subsets except CD8 and γδ T cells; CD8⁺, CD8 T cells; CD19⁺, naïve and memory B-cells, CD56⁺, resting and activated NK cells. Median CIBERSORT-inferred fractions for lung adenocarcinoma GEPs, shown in Fig. 3a, were determined from two publicly available microarray datasets, GSE7670⁷⁷ and GSE10072⁷⁸.

Patient samples

All aspects of this study were approved by the Stanford Institutional Review Board in accordance with the Declaration of Helsinki guidelines for the ethical conduct of research, and all patients involved provided informed consent. For Fig. 3a, fresh human lung tumor samples were obtained from Stanford Tissue Bank. For tissue microarray analyses (Fig. 4c,e, Supplementary Figs. 7c–h), patient samples were retrieved from the surgical pathology archives at the Stanford Department of Pathology and linked to a clinical database using the Cancer Center Database and STRIDE Database tools from Stanford.

Human lung dissociation and flow cytometry

Fresh human lung tumor samples were cut into small pieces and dissociated into single cell suspensions by 45 minutes of Collagenase I (STEMCELL Technologies) digestion. Dissociated single cells were suspended at 1 × 10⁷ per mL in staining buffer (HBSS with 2% heat-inactivated calf serum). After 10 minutes of blocking with 10 μg μl⁻¹ rat IgG, the cells were stained for at least 10 minutes with the antibodies listed below. After washing, stained cells were re-suspended in staining buffer with 1μg/ml DAPI, analyzed, and sorted with a FACS Aria II cell sorter (BD Biosciences). Antibodies used for experiments related to Fig. 3a: CD45-A700, CD14-PE, CD8-APC, CD4-FITC, CD56-PE-cy7, and CD19-PerCP-cy5.5. Antibodies used for enumeration of plasmacytic cells: CD45-PE-cy7, CD20-PerCP-cy5.5, CD138-PE, CD38-APC, CD19-A700, and CD27-FITC. All antibodies were obtained from BioLegend.

Tissue Microarray (TMA) cohort

We reviewed patients with lung cancer who had surgically treated disease and paraffin embedded samples from 1995 through June, 2010 for inclusion. Patients with recurrent or metastatic disease samples only were excluded. Medical charts were reviewed to clinically annotate the tumor specimens with demographic, operative procedures, imaging data, and follow-up. Pathology reports were reviewed to confirm specimen type, site, pathology, stage, histology, invasion status and operative procedure. Treated samples (neoadjuvant therapy) were excluded, resulting in a final analysis cohort of 187 pre-treated lung adenocarcinoma tumor specimens with follow-up data.

TMA cohort follow-up

Recurrence was defined by imaging or biopsy and patients with advanced disease or who did not have at least 6 months of follow-up were censored for further analyses. The National Death Index (NDI) was used to define vital status through October 30, 2010. Patients not dead were assumed to be alive except for those who had left the country or were from other countries (who were censored) since the NDI relies on a social security number for vital status assessment. Synchronous tumors resected over time were eligible for prognostic assessment in patients with two primaries.

TMA construction

The Stanford Lung Cancer TMA was developed from surgical specimens that contained viable tumor from duplicate slides that were reviewed by a board-certified pathologist (R.B.W.). The pathologist was not blinded to sample identity. The area of highest tumor content was marked for coring blocks corresponding to the slides. We used 2 mm cores to build the tissue microarray. These cores were aligned by histology and stage and negative controls were taken from the West Lab and included a variety of benign and malignant tissues (65 cores) that included normal non-lung tissue (12 cores), abnormal non-lung tissue (13 cores), placental markers (23 cores) and normal lung (17 cores). Normal lung consisted of a specimen adjacent, but distinct, from tumor over the years 1995 through 2010 to assess the variability of staining by year. OligoDT analysis was performed on the finished array to assess the architecture of selected cores and adequacy of tissue content prior to target IHC analysis. A co-registered Hematoxylin & Eosin (H&E) slide was used as well to verify tumor location for cases where this was unclear on initial inspection.

TMA immunohistochemistry

MPO (DAKO) and CD20 (clone L26, DAKO) immunohistochemistry performed on 4 mm sections using the Ventana BenchMark XT automated immunostaining platform (Ventana Medical Systems/Roche, Tucson, AZ).

TMA RNA in situ hybridization

The RNA in situ hybridization probe for IGKC was designed against chr2: 88,937,790–88,938,290 (hg18) using primer 5′ – CTG TTG TGT GCC TGC TGA AT – 3’ and the T7 promoter-tagged primer 5’ – CTA ATA CGA CTC ACT ATA GGG TTA AAG CCA AGG AGG AGG AG – 3’. RNA in situ hybridizations were performed on TA369, as described previously⁷⁹.

TMA microscopy

All slides were scanned at 20x on an Ariol imaging analysis system (originally built by Applied Imaging).

TMA staining quantification and analysis

To facilitate consistency and reproducibility in quantitating TMA staining patterns, we evaluated the performance of GemIdent⁸⁰, a supervised in silico image segmentation system. As an initial exercise, we trained GemIdent on a single lung adenocarcinoma specimen to recognize both IGKC stains and non-tissue background (white space). GemIdent was then applied to 10 TMA specimens to generate separate image masks of both IGKC localization and non-tissue background (i.e., “empty space”). A custom Perl script was used to process each image mask and quantify the staining area of IGKC for each specimen (by first removing non-tissue white space to calculate the surface area of each tissue). To test the utility of this approach, a board-certified pathologist (RBW) scored IGKC for the same 10 specimens. The pathologist had no knowledge of the results from automated staining, but was not blinded to sample identity. Both assessments were highly correlated (R² = 0.98; Supplementary Fig. 8c). In a separate exercise, two independent operators trained GemIdent on distinct CD20-stained specimens. CD20-stained fractions were then quantified across the entire TMA (n = 187 lung adenocarcinomas) and results were processed as described above. The concordance between independent operators was very high (R² ~ 1; Supplementary Fig. 8d). These data support the utility of GemIdent coupled with image post-processing for automated scoring of TMA specimens. We applied this approach to quantitatively score IGKC, CD20, and MPO for all lung adenocarcinoma TMA specimens (e.g., see Supplementary Fig. 8a).

Comparison between TALs and circulating leukocytes

Among patients with available perioperative circulating leukocyte (lymphocyte and PMN) counts, we analyzed the sample closest to the date of procedure (DOP), within −120 to +28 days, where precedence was given to preoperative samples (total n = 48 lung adenocarcinoma patients). As shown in Supplementary Table 7, no relationships were found between circulating leukocyte (CL) levels and TALs quantified on the TMA. Moreover, while the ratio of MPO to IGKC levels remained significantly prognostic within this patient subset (P = 0.02), CL levels had no significant relationship to survival.