NITUMID: Nonnegative matrix factorization-based Immune-TUmor MIcroenvironment Deconvolution

Abstract

Motivation: A number of computational methods have been proposed recently to profile tumor microenvironment (TME) from bulk RNA data, and they have proved useful for understanding microenvironment differences among therapeutic response groups. However, these methods are not able to account for tumor proportion nor variable mRNA levels across cell types.

Results: In this article, we propose a Nonnegative Matrix Factorization-based Immune-TUmor MIcroenvironment Deconvolution (NITUMID) framework for TME profiling that addresses these limitations. It is designed to provide robust estimates of tumor and immune cells proportions simultaneously, while accommodating mRNA level differences across cell types. Through comprehensive simulations and real data analyses, we demonstrate that NITUMID not only can accurately estimate tumor fractions and cell types' mRNA levels, which are currently unavailable in other methods; it also outperforms most existing deconvolution methods in regular cell type profiling accuracy. Moreover, we show that NITUMID can more effectively detect clinical and prognostic signals from gene expression profiles in tumor than other methods.

Availability and implementation: The algorithm is implemented in R. The source code can be downloaded at https://github.com/tdw1221/NITUMID.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

Diagram of NITUMID. (a) We applied CIBERSOT (Newman et al., 2015) and immunoStates (Vallania et al., 2018) on the TCGA SKCM and LUAD data and selected the intersection of the most abundant cell types from these two methods as our component cell types. (b) We curated a list of 53 signature genes for the 11 cell types and obtained their mean expression profiles in each cell type. (c) We trichotomized the signature genes expression profile matrix into matrix A. (d) Illustration of the NMF framework in NITUMID: for the 53 by N gene expression matrix Y, we input the guide matrix A, and factorize matrix Y into a 53 by 11 matrix W and an 11 by N matrix H. (e) We designed and implemented a consistency-based criterion to choose the model’s tuning parameters for different datasets. See Supplementary Section SA.5 for details

Fig. 2.

Comparison of NITUMID, immunoStates, CIBERSORT, xCell and EPIC’s performances on in silico immune cell mixture, measured by Pearson correlation between estimated cell fractions and true cell fractions. (a) In silico immune cell mixture is generated from our microarray training datasets. (b) In silico immune cell mixture is generated from melanoma tumor microenvironment (TME) scRNA-Seq data (Tirosh et al., 2016). (c) In silico immune cell mixture is generated from 10× scRNA-Seq immune cell profiles from Zheng et al. (2017)

Fig. 3.

Estimated mean mRNA levels. (a) Estimated mRNA levels by cell type from the W matrix of the four bulk melanoma datasets. (b) CD8+ T cells mRNA levels for the CTLA4 dataset, by response group. (c) Macrophages mRNA levels for the HUGO dataset by response group. (d) Macrophages mRNA levels for the NIVO dataset by the response group and treatment status, partial responding (PR), stable disease (SD) and progressive disease (PD) groups all showed higher levels compared with the complete responding (CR) group

Fig. 4.

Informative immune cell fractions changes identified by NITUMID. (a) Estimated CD8+ T cell fraction by response group and treatment status (NIVO). (b) Estimated NK/NKT cell fraction by response group and treatment status (NIVO). (c) Survival between high CD8+ T cell and low CD8+ T cell groups in the TCGA SKCM cohort, the high group includes 92 samples from the first CD8+ fraction quantile (top 25%), the low group includes the fourth quantile (bottom 25%). (d) Survival difference between samples with high CD8+ T and high NK/NKT cell proportions and those with low proportions in both (TCGA SKCM). High-high group consisting of patients whose CD8+ T and NK/NKT cell fractions are both in top 25%, the low-low group consisting of patients whose CD8+ T and NK/NKT cell fractions are both in the bottom 25%

Fig. 5.

Additional immune cell fractions-related results. (a) Survival curves for the high and low CD8+ T cell groups for TCGA SKCM stratified by the median CD8+ T cell fraction (184 samples each group). (b) Survival curves for the high and low NK/NKT cell groups for TCGA SKCM stratified by the median NK/NKT fraction (184 samples each group). (c) TCGA SKCM survival stratified by estimated CD8+ T cell and NK/NKT fractions. The high-high group are 19 samples whose both fractions are in top 50%, while the low-low group are 19 samples whose both fractions are in bottom 50%. (d–e) Estimated Macrophage and CD4+ T cell fractions by response group and treatment status for the NIVO. (f) Survival curves for the high and low macrophages groups for TCGA SKCM stratified by the median macrophages fraction. (g) Estimated CD8+ T cell fraction by response group for the CTLA-4 data (Allen et al., 2015)