TIMER2.0 for analysis of tumor-infiltrating immune cells.

Tumor progression and the efficacy of immunotherapy are strongly influenced by the composition and abundance of immune cells in the tumor microenvironment. Due to the limitations of direct measurement methods, computational algorithms are often used to infer immune cell composition from bulk tumor transcriptome profiles. These estimated tumor immune infiltrate populations have been associated with genomic and transcriptomic changes in the tumors, providing insight into tumor-immune interactions. However, such investigations on large-scale public data remain challenging. To lower the barriers for the analysis of complex tumor-immune interactions, we significantly improved our previous web platform TIMER. Instead of just using one algorithm, TIMER2.0 (http://timer.cistrome.org/) provides more robust estimation of immune infiltration levels for The Cancer Genome Atlas (TCGA) or user-provided tumor profiles using six state-of-the-art algorithms. TIMER2.0 provides four modules for investigating the associations between immune infiltrates and genetic or clinical features, and four modules for exploring cancer-related associations in the TCGA cohorts. Each module can generate a functional heatmap table, enabling the user to easily identify significant associations in multiple cancer types simultaneously. Overall, the TIMER2.0 web server provides comprehensive analysis and visualization functions of tumor infiltrating immune cells.

INTRODUCTION

Tumor-infiltrating immune cells are important in cancer treatment efficacy and patient prognosis (1–5). The composition of immune cells in the tumor microenvironment contributes to tumor heterogeneity, and creates interesting yet challenging complexities when investigating dynamic interactions between immune and cancer cells (6). Gene expression profiling using microarrays or RNA sequencing (RNA-seq) are mature methods for tumor characterization, which have been widely used to generate a wealth of transcriptomics profiles in many cancer types. While informative, tumor transcriptomics data do not immediately indicate immune cell compositions, which instead require computational inference.

The available computational algorithms (7–12) for immune infiltration estimation fall into two main categories: gene signature- and deconvolution-based approaches (13). Gene signature-based approaches utilize a list of cell-type-specific gene sets. By using the expression values of these signature gene sets in tissue samples, these models infer every cell type independently, either by performing enrichment analysis of the gene sets (7) or by aggregating them into an abundance score (8). Deconvolution methods define the problem as mathematical equations that model the gene expression of a tissue sample as the weighted sum of the expression profiles from the cells in the population mix (9–12). These two complementary categories of algorithms have demonstrated variable performance advantages in estimating specific immune cell types in different tumors. We hope that evaluating estimations from multiple algorithms might help the user gain more comprehensive and robust results. This motivated us to improve TIMER, a web server that we developed previously (14).

The original TIMER web server included only one immune infiltration estimation algorithm (14). Guided by user feedback, we have developed an enhanced version of TIMER (TIMER2.0) that integrates multiple state-of-the-art algorithms for immune infiltration estimation. These algorithms were applied to the expression profiles of the Cancer Genome Atlas (TCGA) tumors, allowing users to explore various associations between immune infiltrates and genetic features in the TCGA cohorts (15). We also implemented several modules to facilitate investigation of features of interest in different cancer types. TIMER2.0 also provides an ‘estimation component’ to infer immune infiltrates on user-provided expression profiles using multiple algorithms.

New developments

Cellular composition estimation

To make reliable immune infiltration estimations, TIMER2.0 utilizes the immunedeconv (16), an R package which integrates six state-of-the-art algorithms, including TIMER (9), xCell (7), MCP-counter (8), CIBERSORT (10), EPIC (11) and quanTIseq (12). These algorithms have been systematically benchmarked (16), and each was found to have unique properties and strengths. For example, different tissue types may induce distinct cancer-cell intrinsic expression and create different immune contexts (17). TIMER is the only method that takes tissue specificity into account when estimating immune cell populations, although it only makes estimations on six immune cell types. xCell can make estimations on the higher number of different immune cell types but may fail to detect signals from homogeneous samples. CIBERSORT deconvolves more detailed subsets of T-cell signatures. EPIC and quanTIseq have the advantage of directly generating scores interpreted as cell fractions. Uniform pre-processing of all the TCGA samples using these algorithms enables users to visualize all the estimations together to reach more confident conclusions.

Functional heatmap table

Users commonly wish to explore associations in all cancer types simultaneously to contrast the cancer type of interest. Using the old version of TIMER, the user had to manually query cancer type one-by-one. We improved this feature in TIMER2.0 by presenting a color-coded heatmap table to combine all calculation information for each of the input features and cancer types. In this table, red blocks indicate statistically significant positive values and the blue blocks indicate statistically significant negative values. Gray denotes non-significant results. The user can order rows by each column through clicking the column header. Detailed information about a relationship can be found intuitively by clicking on the corresponding entry. Thus, the functional heatmap table enables the user to quickly identify significant associations and find common trends among different cancer types or a unique relationship in a specific cancer type.

Estimation for mouse data

Mouse models are widely used in research of human disease, although there are few computational algorithms available to estimate immune infiltration for mouse tissues. Recently Petitprez et al. developed a new version of MCP-counter for mouse, called mMCP-counter (18), which we have incorporated into the mouse model of the ‘estimation component’. To perform analyses with the other five algorithms that were developed to quantify the immune infiltrates for human samples, TIMER2.0 converts mouse gene IDs to the orthologous human gene IDs. We note that estimations on mouse samples should be interpreted with caution since human and the mouse tumors may have different immune properties.

RESULTS

Overview of TIMER2.0 web server

TIMER2.0 is a freely available web server to the research community, and it is implemented by the Shiny web framework for R (version 3.6.1). TIMER2.0 consists of three major components: immune, exploration and estimation. The immune component contains four modules that allow users to investigate the association between estimated immune infiltrates and gene expression, somatic mutations, somatic copy number alterations (sCNAs) and clinical outcomes in the TCGA cohorts (Figure 1). The exploration component has four modules that allow users to find cancer-related associations in TCGA. Users can compare the expression of a gene between tumor and normal and find the associations between the expression of one gene with patient survival, or the mutation status and expression level of other genes. Given the inputs for each module in the immune component and the exploration component, a functional heatmap table presents the association between each input feature and cancer type. All TCGA tumor data, including transcriptome profiles, somatic mutation calls, somatic copy number variations and patient clinical outcomes, are collected from the GDAC firehose website (http://firebrowse.org/). TIMER2.0 extracts raw counts and Transcripts Per Million (TPM) from RSEM results (19). The estimation component can infer immune cell infiltration of user-provided expression profiles using the six estimation algorithms (Figure 2). Instruction for each module is provided on the module page, and an overview of the entire website with a tutorial video is available on the home page.

Figure 1.

An illustration of immune component outputs. Upon a query to the immune component, TIMER2.0 evaluates associations between immune infiltrates and genetic or clinical features and displays the results as a functional heatmap (A). Upon a user clicking on the heatmap, each module subsequently generates figures (B–E) showing the detailed information about the corresponding relationship. (A) An example of functional heatmap table generated by the ‘Gene Module’ shows the association between PDCD1 expression and immune infiltration level of multiple CD8+ T cell types estimated by all six algorithms across TCGA cancer types. The red indicates a statistically significant positive association, and the blue indicates a statistically significant negative association. Gray denotes a non-significant result. (B) An example of scatter plots from the ‘Gene Module’. Correlation of PDCD1 expression with tumor purity (left) and with the infiltration level of CD8 T cell estimated by TIMER (right) in lung squamous carcinoma. (C) An example of the violin plot from the ‘Mutation Module’ displays the difference in TIMER-estimated CD8 T cell infiltration levels between tumors with mutant or wild-type TP53 in bladder cancer. (D) An example of the violin plot from the ‘sCNA Module’ visualizes the difference of CD8 T cell infiltration level estimate among tumors with different sCNA status of PIK3CA gene in head neck cancer. (E) An example of the Kaplan–Meier plot from the ‘Outcome Module’ shows the difference of overall survival among patients stratified by both the estimated infiltration level of CD8+ T cell and PDCD1 expression level in breast cancer.

Figure 2.

An illustration of outputs for the estimation component. TIMER2.0 provides estimations of immune infiltration levels for user-provided tumor profiles using six state-of-the-art algorithms. An illustration of expression profiles (TPM-normalized) of five randomly selected TCGA LUAD samples into the estimation component. (A) The data table presents the different immune cell type infiltration estimated by multiple algorithms. (B) For the eight immune cell types for which all six algorithms can estimate the abundance, TIMER2.0 draws a multi-panel bar plot showing differences of their infiltration level estimated by different algorithms among samples. Of the six algorithms, three algorithms (CIBERSORT in original mode, quanTIseq and EPIC) generate values comparable within the same sample. Based on estimations from these algorithms, TIMER2.0 also presents (C) a multi-panel pie plot showing the proportion of immune cell types in each sample. The eight previously described immune cell types are highlighted in different colors.

Immune component

Gene Module

The relationship between tumor gene expression and immune infiltration is of great interest to many researchers. The ‘Gene Module’ allows a user to identify such relationships in a fast, comprehensive and unbiased way. If a user selects a gene and an immune cell type, TIMER2.0 generates a heatmap table of the Spearman’s correlations between the expression of the input gene and the abundance of the immune cell type as well as its subtypes in the 59-cell hierarchy across all cancer types (Figure 1A). This display allows the user to view whether the expression of their input gene is consistently associated with the immune cell type across different cancer types using different estimation algorithms. When ‘purity adjusted’ option is selected, clicking on any of the numbers in the table will return two scatter plots showing (i) the correlation of the given gene expression with tumor purity (the proportion of cancer cells in a sample) (20) and (ii) the association of the gene expression with the immune cell type (Figure 1B).

Mutation Module

Somatic mutations play crucial roles in tumorigenesis and profoundly influence patient treatment response and survival (21,22). Researchers are often interested in the effect of gene mutations on tumor immune environment, specifically immune cell infiltration changes. The ‘Mutation Module’ allows users to analyze and visualize the effect of non-synonymous somatic mutations on immune cell infiltration across multiple cancer types and immune cell types simultaneously. Given the input gene, TIMER2.0 displays a bar plot showing mutation frequency of the gene for each TCGA cancer type. The heatmap table with embedded violin plots presents differential immune infiltration levels between tumors with the input gene mutated and tumors without the input gene mutated (Figure 1C).

sCNA Module

Genes with somatic copy number alterations (sCNAs) are hallmarks of tumorigenesis and progression (23), and could influence immunotherapy response (24). The ‘sCNA Module’ allows a user to compare immune infiltration distribution by the sCNA status of a gene across TCGA cancer types. Given the input gene, TIMER2.0 displays a stacked bar plot showing the relative proportion of different sCNA states of the gene for all TCGA cancer types. TIMER2.0 requires the user to specify a ‘deep deletion’ or ‘high amplification’ alteration status of the gene, as defined by GISTIC2.0 (25), to compare with the ‘diploid/normal’ status. The heatmap table displays the log-fold changes of immune infiltration levels between the specified alteration group and the normal one. With a click on an entry of the table, the user can view the immune infiltration distribution between different sCNA status of the gene on a violin plot for pairwise comparisons of normal group with each alteration group (Figure 1D).

Outcome Module

Tumor-infiltrating T cells (26,27) and B cells (28,29) have been reported to affect patient clinical outcome. The clinical relevance of the immune infiltrates can be learned by examining associations between immune infiltration level and patient clinical outcome. The ‘Outcome Module’ allows a user to quickly evaluate the association of immune subset abundance and patient survival across TCGA cancer types. Briefly, TIMER2.0 uses a Cox proportional hazard model to evaluate the effects of immune infiltration level, specific gene expression on patient clinical outcome. User-selected variables, such as patient age, tumor stage, gene expression, with the exception of immune cell types, are considered as confounding factors for adjustment in the regression model. The heatmap table shows the normalized coefficients of immune cell types. Each entry on the heatmap table stands for an independent Cox model that includes only one immune cell type and the confounding factors for one cancer type. The user can click on an entry to examine how the immune cell (high versus low) and gene expression (high versus low, if the user inputs a gene of interest) are associated with patient survival on Kaplan–Meier (KM) curves (Figure 1E).

Exploration component

Tumor expression of certain genes such as oncogenes (30) or cancer/testis genes (31,32) may also be associated with tumor therapy response and patient clinical outcome. To investigate cancer-related associations between tumor features in the TCGA cohorts, TIMER2.0 provides the ‘Exploration Component’ with four modules.

Gene_DE Module

To evaluate whether a gene has therapeutic target potential, the gene is commonly uniquely expressed or highly expressed in tumors than most normal tissues. The ‘Gene_DE Module’ allows users to rapidly compare a gene’s expression level between tumor and matched normal tissues across all TCGA cancer types. When a user inputs a gene, TIMER2.0 displays box plots of the gene expression distribution in tumors versus normal tissues across all cancer types. Significance is determined by differential gene expression analysis using edgeR (33).

Gene_Outcome Module

Gene expression-based approaches have been developed to predict survival and therapy response (32,34,35). The ‘Gene_Outcome Module’ allows a user to simultaneously evaluate the gene expression association of multiple genes on patient survival across TCGA cancer types. Similar to the ‘Outcome Module’, TIMER2.0 utilizes the Cox proportional hazard model for association evaluation. Variables belonging to clinical or histological features are considered as confounding factors in the regression model. Once a user submits the request, TIMER2.0 draws a heatmap table of the normalized coefficients of the input gene expression in the Cox model. Each entry on the heatmap represents an independent Cox regression coefficient, which associates each gene expression and other co-variant to the hazard ratio for the corresponding cancer type, respectively. The user can click an entry on the heatmap to investigate the association between the gene expression and survival in the KM curves. For example, high expression levels of CD8A, CD8B and the cytotoxic T lymphocytes (CTLs) marker genes GZMA, GZMB and PRF1 are associated with improved survival in skin cutaneous melanoma (SKCM) patients, which is consistent with the anti-tumor function of CTLs.

Gene_Mutation Module

Identifying the functional impact of somatic mutations in the gene expression pattern is critical for precision oncology and drug discovery (36). The ‘Gene_Mutation Module’ allows a user to assess how the mutation status of one gene influences the expression of the mutated gene and other genes across TCGA cancer types. When a user selects one recurrently mutated gene X and enters gene Y on the expression query, TIMER2.0 displays a heatmap of log-fold changes of gene Y expression between tumors with the gene X mutation and tumors without the gene mutation. The user can click an entry on the table to view detailed violin plots of gene expression distribution in the mutant versus wild-type tumors.

Gene_Corr Module

This module allows a user to discover the co-expression pattern of genes across TCGA cancer types. Given one initial gene X of interest and up to 20 other genes, TIMER2.0 generates a heatmap table of Spearman’s correlation of gene expression between gene X and the other input genes. The user can click an entry on the heatmap to view scatter plots of the expression correlation between the two genes.

Estimation component

Given the prevalence of transcriptome studies across immuno-oncology, we improved the capacity for immune infiltration estimation and interpretation of user-provided expression profiles. The ‘estimation component’ accepts a gene by sample expression matrix in a ‘csv’ or ‘txt’ file formatted with standard delimiters. The values in this matrix should be TPM-normalized without log-transformation to fit the requirements of all estimation methods. Once the input file is uploaded, TIMER2.0 automatically shows the number of samples so that the user is able to confirm that the file was uploaded correctly. TIMER2.0 adopts the immunedeconv package (16) to infer immune cell component abundance of the input gene expression data using the six state-of-the-art algorithms described previously. The TIMER algorithm takes tissue-specific effects into account and requires cancer-type information to improve the accuracy of the estimation. If a user selects ‘AUTO’ as the cancer type, TIMER2.0 will automatically assign the sample to the TCGA cancer type with the highest expression similarity.

Once the user uploads the gene expression file and clicks ‘RUN’, TIMER2.0 will run immune infiltration estimation and show a progress bar. The results are summarized in a downloadable data table with the estimated values of all immune infiltrates for each sample from each algorithm (Figure 2A). It is worth noting that the estimation from different algorithms have different interpretations. While all six methods allow comparisons of the same cell type between samples, CIBERSORT, quanTIseq and EPIC also allow comparisons of the different cell types within the same sample. To illustrate this difference, we examine the two images that TIMER2.0 draw when input sample size is <10. The first image is a multi-panel bar plot comparing the infiltration level between samples for eight common immune infiltrates estimated by each method (Figure 2B). The second image is a pie plot showing the proportion of the eight immune infiltrates estimated by CIBERSORT, quanTIseq and EPIC (Figure 2C).

DISCUSSION

TIMER2.0 is an updated web server with unique features that enable analyses and visualization of tumor immunity and its association with other tumor molecular and clinical features. First, TIMER2.0 estimates immune cell infiltration in TCGA tumors or user-provided transcriptome data using six different computational methods so users might gain more confidence on the results that are consistently observed using different methods. Second, TIMER2.0 helps to find associations between immune infiltration, gene expression, mutations and survival features in the TCGA cohorts. Finally, TIMER2.0 provides user-friendly and interactive visualization to facilitate data exploration. With tumor profiling data from single-cell technologies increasingly available, there will undoubtedly be improvements in defining signatures of tumor-infiltrating immune cells and in methods for computational estimation. We envision continued development of the TIMER web resource to benefit the cancer immunology research community.

DATA AVAILABILITY

TIMER2.0 is freely available at http://timer.cistrome.org.