Correspondence | Open | Published:
Data normalization considerations for digital tumor dissection
Genome Biologyvolume 18, Article number: 128 (2017)
In a recently published article in Genome Biology, Li and colleagues introduced TIMER, a gene expression deconvolution approach for studying tumor-infiltrating leukocytes (TILs) in 23 cancer types profiled by The Cancer Genome Atlas. Methods to characterize TIL biology are increasingly important, and the authors offer several arguments in favor of their strategy. Several of these claims warrant further discussion and highlight the critical importance of data normalization in gene expression deconvolution applications.
Computational approaches for enumerating cell subsets from bulk tissue expression profiles have significant potential for studying tumor cellular ecosystems, including tumor-infiltrating leukocytes (TILs) [1,2,3]. We therefore read with interest the recent Genome Biology article by Li and colleagues in which they introduce TIMER, an in silico method for TIL deconvolution . TIMER relies on prior knowledge of immune signature genes as input and consists of three major steps: (1) gene expression normalization across platforms and sample types; (2) selection of immune signature genes that are negatively correlated with tumor purity; and (3) deconvolution of RNA admixtures using a previously described technique for iterative linear least squares regression (LLSR) . They apply TIMER to the inference of six distinct immune subsets (B cells, CD4 T cells, CD8 T cells, neutrophils, macrophages, and dendritic cells) in The Cancer Genome Atlas (TCGA) bulk tumor expression profiles and investigate links between TIL heterogeneity, tumor genomic features, and survival in 23 cancer types.
Several groups, including ours, have also proposed methods for gene expression deconvolution [1, 3, 5, 6]. We recently described CIBERSORT, an in silico tissue dissection approach that is robust to noise, unknown mixture content, and closely related cell types (collinearity) . Notably, in benchmarking experiments CIBERSORT outperformed other deconvolution methods, including LLSR, and revealed complex associations between 22 distinct immune subsets and outcomes in a pan-cancer meta-analysis [6, 7]. We were therefore surprised by several claims in relation to CIBERSORT.
First, the authors assert that CIBERSORT succumbs to statistical collinearity (i.e., cell subsets with highly correlated expression profiles), leading to biased estimations. Evidence for this argument is primarily based on a simple experiment in which inferred levels of each immune subset were compared by Pearson correlation. After aggregating CIBERSORT results from 22 phenotypes into the same six subsets, the authors compared cross-correlation matrices between TIMER and CIBERSORT on four cancer types. Leukocyte levels estimated by TIMER were almost always positively correlated. According to the authors, positive correlations make intuitive sense because “immune cells work in synergy.” In contrast, the correlations among the six phenotypes estimated by CIBERSORT were largely negative. The authors also observed negative correlations when analyzing more than six cell types with TIMER (i.e., LLSR), stating that negative correlations indicate a technical artifact due to collinearity.
In fact, CIBERSORT mitigates such bias through regularization, as was rigorously demonstrated through a battery of validation experiments [6, 7]. These analyses included an assessment of “deep deconvolution” in which in silico predictions of closely related leukocyte subsets were directly compared against flow cytometry . Moreover, an independent study confirmed that regularization improves the performance of gene expression deconvolution .
We were also surprised by the authors’ claim that hematopoietic cell types should generally track together, especially in light of their diverse functions (innate or adaptive, stimulatory or suppressive, etc.) and migration patterns (circulating, tissue-infiltrating, or tissue-resident) [9, 10]. For instance, while specific hematopoietic subsets infiltrating tumors can be positively correlated in a given tumor type , the expectation of universally positive correlations does not extend to all tumor types, or to all infiltrating immune cells . Separately, age-related lymphomyeloid lineage skewing of hematopoiesis  would be expected to further confound this assumption. Finally, while acute and chronic inflammation can be substrates for tumor initiation, a number of distinct tumor-infiltrating immune cells are known to have either tumor-promoting or anti-tumor properties, and are associated with inverse prognostic correlations with cancer outcomes [9, 14].
We therefore reanalyzed previously published flow cytometry data of leukocyte subsets directly enumerated in peripheral blood mononuclear cells (PBMCs) from healthy donors and in tumor biopsies obtained from patients with lung squamous cell carcinoma (LUSC) [6, 7]. When we quantified each immune subset as a fraction of total leukocyte content, many of the pairwise correlations were negative, as were the mean correlation coefficients (Fig. 1a and b), consistent with CIBERSORT. However, when we instead considered absolute TIL levels in the same lung tumors, most of the correlations were positive (Fig. 1c), likely reflecting differences in tumor purity. Thus, data normalization in solid tumors significantly impacts the assessment of TIL heterogeneity and composition.
Given these results, we suspected that tumor purity would explain the discrepancy between TIMER and CIBERSORT. Indeed, after examining the TIMER source code, we found that, unlike most previous deconvolution methods including CIBERSORT, TIMER solves the regression problem without normalizing inferred cell subset frequencies to 1. TIMER results are therefore directly influenced by total leukocyte content, which is inversely correlated with tumor purity across TCGA (Fig. 2a). As a result, all six cell types strongly correlate with total leukocyte abundance in nearly every analyzed tumor type (Fig. 2b), making it difficult to discern the intercellular heterogeneity among the leukocyte subsets that variably infiltrate these tumors (Fig. 2c). When TIMER results were instead normalized in relative space (i.e., summing to 1) for each sample, all mean cross-correlation coefficients were negative (Fig. 2d). The inverse held true for CIBERSORT: mean cross-correlation coefficients became positive when we either (1) omitted the sum-to-1 normalization step, or (2) multiplied the normalized results by a separate estimate of overall immune content (Fig. 2e). While we acknowledge that TIMER estimates were not intended to be analyzed in relative space (as described below), the same reasoning should have been applied by Li et al. to CIBERSORT; that is, CIBERSORT relative abundance estimates should not have been directly compared with absolute leukocyte abundance (as in Table S6 from ). Collectively, these data highlight the importance of data normalization in comparing gene expression deconvolution methods.
We and others have previously shown that regression-based gene expression deconvolution can robustly quantify cell-type proportions [5,6,7, 15, 16]. Therefore, we were surprised by the claim that “levels of different cell types are not comparable” in the output of TIMER . Upon further examination of this output, we found disproportionately high levels of rare dendritic cells (DCs) across all 23 cancer types (approximately 50% by inferred fractional abundance; Fig. 2f), suggesting problems with marker gene selection and/or data normalization. We hypothesized that this result might be due to the authors’ use of ComBat  to purge batch effects between two highly distinct sample types: bulk tumors profiled by TCGA and a knowledgebase of purified leukocytes used for signature genes. In support of this hypothesis, we found that the number of DC reference profiles in the knowledge base was strongly correlated with the expression of DC marker genes in normalized tumors (Fig. 2g). Further analysis revealed a strong association between predicted abundance in TCGA (Fig. 2f) and representation in the knowledge base for all six leukocyte subsets (Fig. 2h). Thus, misapplication of ComBat distorted important biological signals that correlate with experimental batches , preventing TIMER from estimating cell type proportions.
Separately, we wish to address the claim that CIBERSORT is only applicable to microarray data [2, 4]. While microarray datasets were indeed the focus of our previous studies, this is not an inherent restriction of the deconvolution algorithm itself, which is platform agnostic. In fact, the analytical assumptions made by CIBERSORT are likely to hold for any mixture that can be modeled as a linear sum of its parts and for which an appropriate signature matrix exists. Such mixtures include RNA-Seq datasets, as others have already shown for bulk tumor profiling [19,20,21], and for single-cell RNA-Seq profiling , as well as other genomic features associated with cell lineage . For example, CIBERSORT was recently used to enumerate hematopoietic subsets in bone marrow biopsies from healthy and diseased patients based on genomic patterns of nucleosome accessibility profiled by ATAC-Seq , demonstrating its broad applicability.
Finally, in response to this correspondence, Li et al.  have made a number of new claims that warrant clarification. In order to comprehensively address these claims, we have included a detailed point-by-point response, including new analyses, in Additional file 1: Figures S1 and S2). We summarize three key points:
The authors continue to ignore the significant impact of data normalization on deconvolution results, stating that CIBERSORT produces nonbiological negative correlations mainly due to collinearity. They dismiss the notion that regularization can help combat collinearity (despite significant literature on the topic , e.g. ridge regression, and ), and offer a flawed analysis to support their claim consisting of synthetic mixture datasets that are improperly defined since the mixed populations do not sum to 100% and are therefore unsuitable for addressing this topic (Additional file 1: Figure S1a and b).
Furthermore, Li et al. use a single flow cytometry experiment (Fig. 1a) to argue that closely related immune cell types should be positively correlated in abundance, whether in blood or in tumors. Since the default version of CIBERSORT produced negative correlations for the same cell types when enumerated in solid tumor biopsies, Li et al. claim these results contradict our own experimental data in Fig. 1a and are likely due to collinearity. In making these arguments, Li et al. disregard some fundamental immunological principles governing leukocyte migration patterns (see above and Additional file 1), relevant prior literature (e.g., Fig. 3a in ), and the main point of this correspondence (e.g., Fig. 2e). To further illustrate the impact of data normalization on deconvolution results, we extended our CIBERSORT analysis in Fig. 2e to 22 immune subsets (i.e., LM22 ) in TCGA LUSC tumors. As expected, the majority of pairwise correlations were positive when relative abundance estimates were scaled by total immune content, including correlations between closely related cell types (e.g., naïve versus memory B cells; Additional file 1: Figure S2). Moreover, we observed no significant association between pairwise correlations of leukocyte estimates in tumors and pairwise correlations of corresponding expression profiles in LM22 (Additional file 1: Figure S2). Therefore, leukocyte behavior is highly complex and unlikely to be distilled into simplistic comigration patterns without significant further investigation, especially without consideration of data normalization.
Finally, Li et al. claim that up to 25% of LM22 genes are positively correlated with tumor purity and, as a result, they contend that CIBERSORT’s model is “frequently violated” when applied to tumors. Unfortunately, the authors ignore critical details of the algorithm and the LM22 signature matrix design. They also fail to consider many important factors in the interpretation of their own analyses, including the statistical significance, magnitude, and distribution of correlation coefficients, and the impact of positively correlated LM22 genes on CIBERSORT results. When considering these variables, most of the significant positive correlations are of modest magnitude (e.g., 70% with r < 0.2) and only a small minority of LM22 genes are significantly positively correlated with tumor purity (3% with r > 0.2, approximately 0% with r > 0.4; Additional file 1: Figure S1c). Furthermore, since exclusion of all significantly positively correlated genes from LM22 had virtually no impact on tumor deconvolution performance (Additional file 1), we observed no empirical evidence consistent with the above claim.
In summary, our results address key conclusions in Li et al. [4, 24] and emphasize the importance of data normalization in deconvolution analyses. In particular, deconvolution methods cannot be meaningfully compared without taking normalization differences into account. By focusing on relative measures of TIL content in previous work , we avoided the confounding impact of tumor purity . This approach has precedence in prior literature, particularly since many prognostic associations are more robust when defined as ratios of functionally distinct TILs (e.g., CD8 T cells versus Tregs, lymphocytes versus neutrophils, etc.) [7, 27, 28]. Whether absolute or relative measures of TIL abundance better capture tumor immunology in clinical settings remains an important consideration for future studies.
Linear least squares regression
Lung squamous cell carcinoma
Peripheral blood mononuclear cell
The Cancer Genome Atlas
Newman AM, Alizadeh AA. High-throughput genomic profiling of tumor-infiltrating leukocytes. Curr Opin Immunol. 2016;41:77–84.
Aran D, Butte AJ. Digitally deconvolving the tumor microenvironment. Genome Biol. 2016;17:175.
Shen-Orr SS, Gaujoux R. Computational deconvolution: extracting cell type-specific information from heterogeneous samples. Curr Opin Immunol. 2013;25:571–8.
Li B, Severson E, Pignon JC, Zhao H, Li T, Novak J, et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol. 2016;17:174.
Abbas AR, Wolslegel K, Seshasayee D, Modrusan Z, Clark HF. Deconvolution of blood microarray data identifies cellular activation patterns in systemic lupus erythematosus. PLoS One. 2009;4, e6098.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7.
Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med. 2015;21:938–45.
Mohammadi S, Zuckerman N, Goldsmith A, Grama A. A critical survey of deconvolution methods for separating cell-types in complex tissues. arXiv. 2015: 1510.04583. https://arxiv.org/abs/1510.04583.
Fridman WH, Pages F, Sautes-Fridman C, Galon J. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12:298–306.
Shiao SL, Ganesan AP, Rugo HS, Coussens LM. Immune microenvironments in solid tumors: new targets for therapy. Genes Dev. 2011;25:2559–72.
Gao Q, Qiu SJ, Fan J, Zhou J, Wang XY, Xiao YS, et al. Intratumoral balance of regulatory and cytotoxic T cells is associated with prognosis of hepatocellular carcinoma after resection. J Clin Oncol. 2007;25:2586–93.
Stoll G, Bindea G, Mlecnik B, Galon J, Zitvogel L, Kroemer G. Meta-analysis of organ-specific differences in the structure of the immune infiltrate in major malignancies. Oncotarget. 2015;6:11894–909.
Rossi DJ, Bryder D, Zahn JM, Ahlenius H, Sonu R, Wagers AJ, Weissman IL. Cell intrinsic alterations underlie hematopoietic stem cell aging. Proc Natl Acad Sci U S A. 2005;102:9194–9.
Terzić J, Grivennikov S, Karin E, Karin M. Inflammation and colon cancer. Gastroenterology. 2010;138:2101–14.e2105.
Gong T, Hartmann N, Kohane IS, Brinkmann V, Staedtler F, Letzkus M, et al. Optimal deconvolution of transcriptional profiling data using quadratic programming with application to complex clinical blood samples. PLoS One. 2011;6, e27156.
Zhong Y, Wan YW, Pang K, Chow LM, Liu Z. Digital sorting of complex tissues for cell type-specific gene expression profiles. BMC Bioinformatics. 2013;14:89.
Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8:118–27.
Jaffe AE, Hyde T, Kleinman J, Weinbergern DR, Chenoweth JG, McKay RD, et al. Practical impacts of genomic data "cleaning" on biological discovery using surrogate variable analysis. BMC Bioinformatics. 2015;16:372.
Tuong ZK, Fitzsimmons R, Wang SM, Oh TG, Lau P, Steyn F, et al. Transgenic adipose-specific expression of the nuclear receptor RORα drives a striking shift in fat distribution and impairs glycemic control. EBioMedicine. 2016;11:101–17.
Mehnert JM, Panda A, Zhong H, Hirshfield K, Damare S, Lane K, et al. Immune activation and response to pembrolizumab in POLE-mutant endometrial cancer. J Clin Invest. 2016;126:2334–40.
Srinivasan S, Su M, Ravishankar S, Moore J, Head PE, Dixon JB, et al. TLR-exosomes exhibit distinct kinetics and effector function. arXiv. 2016 :1608.08565v1. https://arxiv.org/abs/1608.08565 .
Baron M, Veres A, Wolock SL, Faust AL, Gaujoux R, Vetere A, et al. A single-cell transcriptomic map of the human and mouse pancreas reveals inter- and intra-cell population structure. Cell Syst. 2016;3:346–60.e4.
Corces MR, Buenrostro JD, Wu B, Greenside PG, Chan SM, Koenig JL, et al. Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution. Nat Genet. 2016;48:1193–203.
Li B, Liu JS, Liu XS. Revisit linear regression based deconvolution methods for tumor gene expression data. Genome Biol. 2017. doi:10.1186/s13059-017-1256-5.
Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference and prediction. 2nd ed. New York: Springer; 2009.
Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun. 2015;6:8971.
Sato E, Olson SH, Ahn J, Bundy B, Nishikawa H, Qian F, et al. Intraepithelial CD8+ tumor-infiltrating lymphocytes and a high CD8+/regulatory T cell ratio are associated with favorable prognosis in ovarian cancer. Proc Natl Acad Sci U S A. 2005;102:18538–43.
Templeton AJ, McNamara MG, Seruga B, Vera-Badillo FE, Aneja P, Ocaña A, et al. Prognostic role of neutrophil-to-lymphocyte ratio in solid tumors: a systematic review and meta-analysis. J Natl Cancer Inst. 2014;106:dju124.
Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30:413–21.
Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Abbas AR, Baldwin D, Ma Y, Ouyang W, Gurney A, Martin F, et al. Immune response in silico (IRIS): immune-specific genes identified from a compendium of microarray expression data. Genes Immun. 2005;6:319–3.
AMN and AAA wrote the manuscript with input and edits from all authors. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.