Evaluation of immune responses in individual immune cell types is important for the development of new medicines. Here, we propose a computational method designated ICEPOP (Immune CEll POPulation) to estimate individual immune cell type responses from bulk tissue and organ samples. The relative gene responses are scored for each cell type by using the data from differentially expressed genes derived from control- vs drug-treated sample pairs, and the data from public databases including ImmGen and IRIS, which contain gene expression profiles of a variety of immune cells. By ICEPOP, we analysed cell responses induced by vaccine-adjuvants in the mouse spleen, and extended the analyses to human peripheral blood mononuclear cells and gut biopsy samples focusing on human papilloma virus vaccination and inflammatory bowel disease treatment with Infliximab. In both mouse and human datasets, our method reliably quantified the responding immune cell types and provided insightful information, demonstrating that our method is useful to evaluate immune responses from bulk sample-derived gene expression data. ICEPOP is available as an interactive web site (://vdynamics.shinyapps.io/icepop/) and Python package (://github.com/ewijaya/icepop).
The immune system consists of many different types of cells, which function in concert to mediate immune responses. However, the interaction complexity hinders the evaluation of the individual immune cell type responses to drug treatment. For example, especially in a bulk sample like a whole organ or biopsy sample, the evaluation requires in vitro analysis of each isolated immune cell type, which is usually effortful.
Recent large-scale multi-omics technologies have accelerated the development of new in silico methods to understand complex immune responses1,2 and immune-drug interactions3. The CTen (Cell Type Enrichment)4 and GSEA (Gene Set Enrichment Analysis)5 programs have been used to analyse immune cell responses using large-scale gene expression data. These methods calculate enrichment scores according to the presence of the query genes over the gene set references. Various defined gene sets are available and have been widely used for a variety of biological analyses including immune responses.
These methods use the annotation enrichment of various biological processes for their primary output. They usually do not account for gene expression information in their calculations. The methods are qualitative and hence are not suitable for quantitative estimation analyses, especially of immune cell responses from whole-organ or tissue-derived samples, because whole samples usually contain many different cell types. In cases where a single cell population responds specifically to drug treatment, these qualitative methods are also practical for sample characterisation because most of the gene responses derive from a single cell population (see Supplementary Fig. S14a). In other cases, different cell populations respond simultaneously to drug treatment. In such cases, CTen and GSEA only provide a mixture of annotations from the different cell populations together; these mixed qualitative results are usually difficult to interpret, especially to identify the key responding cell population(s) (see Supplementary Fig. S14b–d and Supplementary Fig. S15).
Whole tissue/organ sample microarray analysis also poses a technical hindrance. During the sample preparation, the whole organ samples can be contaminated with neighbouring organs or tissues. In the ongoing Adjuvant Database Project (://adjuvantdb.nibiohn.go.jp), we have sporadically experienced those unexpected contaminations despite very careful handling. In such cases, the contaminants are unavoidably co-purified with the primary target RNA and introduce detrimental false results in the gene expression analyses.
To solve these problems, we have developed a computational program designated ICEPOP (Immune CEll POPulation) to quantitatively estimate the relative immune cell response for each cell type. ICEPOP uses gene names and the associated expression values from microarray data of mouse and human organ or tissue samples (Fig. 1a). Based on the given gene expression fold-changes and names, ICEPOP calculates the relative response scores for each immune cell type defined by the reference matrix made from public gene expression databases including the ImmGen6 (for mouse) and the IRIS7 (for human), which encompass a wide-range of immunological cell types and states (Fig. 1b). ICEPOP is intended to analyse non-purified bulky samples such as whole organs or biopsy tissue. We integrated a nearby organ-derived gene contamination removing filter. ICEPOP analysis provides two major outputs: (i) bar graph presentation of the relative response score for ImmGen or IRIS defined immune cell types; and (ii) comparison of differently treated samples with circular map presentation, allowing more detailed analysis and interpretation.
In this study, we tested ICEPOP on adjuvant administration studies in mice and two human studies including vaccination of virus-like particles (VLPs) based on human papilloma virus (HPV), and infliximab treatment of inflammatory bowel disease. Insightful characterizations of the immune cell responses are described.
ICEPOP was first tested on two publicly available microarray datasets—GSE63332 and GSE7768—from the Gene Expression Omnibus (GEO) repository (://www.ncbi.nlm.nih.gov/geo/) (Table 1). These data were derived from mouse spleen, normalized as previously described8,9. Contaminating genes from nearby organs, such as the pancreas, were removed by a coefficient of variation (CV) filter (see Methods and Fig. 1b). Then, activation of each cell type was quantified by ICEPOP with several analytical parameters (see Methods).
Examination of the GSE63332 dataset (Table 1) focused on hydroxypropyl-β-cyclodextrin (bCD). bCD is a potent vaccine adjuvant in mice8. Local high concentrations of bCD temporally damage host cells at the administration site. The damaged cells release DNA, which is key in bCD’s adjuvant character8. In this dataset, mice received phosphate buffered saline (PBS) or bCD via intradermal (ID) or intraperitoneal (IP) injection. Six hours later the spleen was removed for microarray analysis of gene expressions. Differentially expressed genes (DEGs) for paired PBS and bCD injected samples were compared (Supplementary Fig. S1a).
In bCD.ID samples, neutrophils showed the highest ICEPOP score among 10 different cell types (Fig. 2a), implicating neutrophils as the most activated cell type in the spleen in this condition. A bar plot of ICEPOP analysis quantitatively presents three analytical scores (see Methods for detailed descriptions). The ICEPOP score represents the relative gene response in each cell type (y-axis of the bar graph). The cell type response threshold (CRT) indicated by the red horizontal line represents the threshold value of cell type response. For example, in bCD.ID, the CRT was 0.049. The ICEPOP scores of six cell types including dendritic cells (ICEPOP score = 0.062), macrophages (0.096), monocytes (0.088), natural killer (NK) cells (0.050), neutrophils (0.398), stem cells (0.055), and stromal cells (0.108) were higher than this CRT, the latter two barely, indicating that the cells responded to bCD ID injection. With the same criteria, B cells, alpha beta T cells, and gamma delta T cells did not respond (Fig. 2a). The third analytical score, the sample response (SR) score, represents the overall sample response (in this case, whole spleen). Based on our null data simulation (Supplementary Fig. S13b, discussed later), an SR score exceeding about 0.3 in the mouse spleen indicates that the sample is substantially responding to the drug. In bCD.ID, SR was 0.723, suggesting bCD ID injection induced considerable biological responses in the spleen (Fig. 2a).
In contrast, in the bCD.IP condition, ICEPOP scores of all cell types and CRT were about 0.1 and 0.098, respectively, indicating virtually no responses (Fig. 2d). The SR of 0.025 also indicated no response to the IP route of bCD administration. Of note, the CRT value almost equaled the theoretical maximum (0.100) of the 10 cell type analysis (see Methods for more detailed description), bolstering the lack of detectable immunological cell responses.
The difference between administration routes was confirmed by histology examination. In the spleen cryosections stained with markers of neutrophils (Ly-6G) and B cells (B220), Ly-6G staining was prominent in bCD.ID (Fig. 2b) but not in bCD.IP (Fig. 2e), indicating neutrophils were activated in the spleen 6 hours after bCD ID administration, but not after IP administration. Correspondingly, there was no difference in B cell staining after these two treatments (Fig. 2c and f). The results were consistent with ICEPOP analysis (Fig. 2a and d). The collective results demonstrated that ICEPOP analysis can provide reliable evaluation of the sample responses at both cell type and whole organ levels.
The examination of the GSE7768 compared lipopolysaccharide (LPS) and monophosphoryl Lipid A (MPL) 6 hours after injection (Table 1). LPS is a prototypic component of the Gram-negative bacterial cell wall. It can be toxic. MPL is a derivative of LPS with low toxicity9,10. In this dataset9, administrations were all performed via intravenous (IV) injection. The PBS, LPS, and MPL groups each comprised three mice; data from the total of nine spleens were used for ICEPOP analysis. DEGs were calculated as PBS vs LPS or PBS vs MPL (Supplementary Fig. S1b). The CV filter (described in later section) was applied and cell type responses evaluated. With the LPS treatment (Fig. 3a), six of ten cell types (dendritic cells, macrophages, monocytes, NK cells, neutrophils, and stromal cells) surpassed the CRT (0.060), indicating LPS responsiveness. The same cell types responded to MPL with CRT of 0.060 (Fig. 3b). Stromal cells showed a stronger response to LPS than to MPL (Fig. 3a vs b).
We further evaluated the gene differences between LPS and MPL samples using the gene level response score of ICEPOP (Supplementary Fig. S2e,m; see Methods for detail). This score was obtained for all genes, and the top 100 scoring genes were identified (Supplementary Fig. S2n). These genes and cell type paired information were plotted in a circular plot (Fig. 3c). In the plot, some genes appeared more than once, since many genes are expressed by multiple cell populations; in the plot, these co-occurring genes are linked with lines. Coloured lines indicate that they appeared in both LPS and MPL samples, whereas those appearing in LPS or MPL only are denoted by grey lines. This circular plot strongly indicated that many of the feature genes with a relatively high gene level response score were in both sample types, suggesting that LPS and MPL induce very similar gene responses in the mouse spleen. On the other hand, LPS induced stronger gene responses in stromal cells for Ubd, Cxcl5, Cxcl1, Serpine1, Atf3, Mmp3, and Timp1, while dendritic cells, macrophages, monocytes, neutrophils, and NK cells showed similar gene responses to LPS and MPL (Fig. 3c). We are unaware of any prior report of the different sensitivity of stromal cells toward LPS and MPL; this may explain the lower biological toxicity of MPL.
Microarray analysis is highly susceptible to contamination. When harvesting the spleen, fragments of neighbouring tissues including the pancreas and duodenum might be included. We observed heavy contamination of at least two MPL data of GSE7768 with pancreas- and duodenum-derived genes, which detrimentally affected the analysis (Fig. 3b vs e).
To lessen the effect of contaminated tissue-derived genes in silico, we applied our originally developed coefficient of variation (CV) reference table (see Fig. 1b). CV is widely used to compare the amount of variation of one variable11. Higher CV in the microarray data of multiple replicates indicates greater dispersion of gene expression among replicates. Therefore, genes with high CV values in control (PBS-treated) samples were more likely to originate from other, sporadically contaminating tissues. To calculate the CV for each gene, we utilized samples of the ongoing Adjuvant Database Project, which contains microarray expression data from PBS-treated control mouse spleen (33 samples), liver (33 samples), and lymph nodes (18 samples) (unpublished data). By analysing these samples, we empirically determined that the genes with a CV value exceeding 1 were more likely to be derived from tissue contamination. This control sample derived CV values strongly correlated with entropy scores of pancreas or intestine in BARCODE 3.012, suggesting that high CV genes are likely from pancreatic or intestinal (duodenum and small intestine) tissues (Supplementary Fig. S3). The same BARCODE analysis of liver and lymph node control samples was not as evocative as spleen (Supplementary Fig. S3), but removing those high-CV genes from the following analysis would expect to stabilize the analysis result, especially for expression value utilizing analysis like ICEPOP. Subsequently, we integrated this ‘CV more than 1 cut-filter’ in ICEPOP as default when the CV reference table is available (see Fig. 1b).
For the GSE7768 contamination analysis, we plotted the gene expression data (x-axis) together with our reference CV values (y-axis) for all nine samples (Supplementary Fig. S4). MPL samples in the third row, especially the first and second replicates, displayed dots at the upper right areas of the plot, and the other samples did not (Supplementary Fig. S4). These genes with high expression and high CV were considered a signal from the contaminating tissues in this dataset. To validate the reliability, we generated a gene expression heat map of these presumed signal genes. We also placed a heat map of tissue specificity entropy score obtained from the BARCODE 3.0 database12, and placed them side-by-side (Supplementary Fig. S5). Most of the genes with high expression and high CV were strongly correlated with the genes derived from pancreatic tissues (pancreas and Islets of Langerhans) and intestinal tissues (duodenum, small intestine, jejunum, ileum, colon, intestinal epithelial cells and caecal tissue). Although not all of the high-CV genes were derived from contaminated tissues (Klk1 and Ighg1 were expressed in all 9 samples), most of the high-CV genes were correlated with pancreas- and duodenum-derived genes. In addition, the same plot of bCD samples (Supplementary Fig. S6a) indicated that they were free of detectable contamination, and their high-CV genes (upper left area in the graph) correlated mostly with spleen, as with the bCD samples (Supplementary Fig. S6b). These results strongly suggest that our CV filtering method is a reliable approach to remove contaminated genes in silico.
The contaminated genes unusually caused high fold-changes and profoundly affected the ICEPOP analysis, especially GSE7768 analysis. With or without of the contamination removal CV filter, the SR of the LPS did not change (before CV filter = 0.504 and after = 0.519) indicating that LPS was not contaminated (Fig. 3a vs d). However, in the MPL, without contamination removal, ICEPOP analysis resulted in almost no response with SR of 0.045 (Fig. 3e). Interestingly, with CV filtering, the same MPL-treated spleen data showed clear responses with SR of 0.503 (Fig. 3b), very similar to that of LPS (Fig. 3a). Taken together, these results confirmed the effectiveness of our CV filtering method, and showed that removing the contaminated tissue-derived genes provides more reliable results.
We also applied our method to human samples, in this case a peripheral blood mononuclear cell (PBMC) microarray dataset from the HPV-16L1 VLP (HPV-VLP) vaccine study13. This dataset (GSE13587) was also obtained from the GEO repository (Table 1). It contains gene expression microarray data of 4 groups: either placebo or HPV-VLP vaccine injection, at two-different times: pre– and post– (see Supplementary Fig. S1c). For the microarray analysis, PBMCs were stimulated in vitro with either medium or HPV-VLP in this dataset. We calculated the DEGs for medium vs HPV-VLP pairs, producing 4 input DEG data including pre-placebo, post-placebo, pre-vaccine, post-vaccine. This dataset focused on the induced genes difference after in vitro stimulation with HPV-VLP in these groups.
ICEPOP analysis indicated dendritic cell and macrophage responses in all 4 groups (Fig. 4). The vaccine status and time point-independent responses were likely to be induced by in vitro HPV-VLP stimulation. HPV-VLP induces dendritic cell activation in vitro14,15,16,17, indicating that HPV-VLP itself is immune-stimulatory, consistent with our results.
In the placebo group, there were no obvious differences between pre- and post- groups (Fig. 4a vs b). This response similarity was also demonstrated by the coloured lines in the circular map (Fig. 4c). In contrast, the same analysis of the vaccine group showed that T cells, NK cells, macrophages, and dendritic cells were more strongly activated in the post-vaccine group compared to that in the pre-vaccine group (Fig. 4d vs e, and f). Importantly, T cell responses were only observed in the post-vaccine group (Fig. 4e and f). This result suggested that HPV-VLP vaccination can successfully induced T cell responses. The listed genes of T cells and NK cells in the circular map (Ifng, Csf2, Gzmb, Socs1, Lta, Lif, Tyms) of the post-vaccine group (Fig. 4f) also suggested that HPV-VLP vaccination induced HPV-specific cellular immune responses18,19,20. Interestingly, gene responses in dendritic cells were also expanded in the post-vaccine (Fig. 4f). This may imply that the vaccination-dependent activation of T cells and NK cells after in vitro VLP stimulation also influenced the activation status of dendritic cells and macrophages in the same culture, possibly through a cytokine and/or a cell-cell contact mechanism(s) among them at the secondary level.
We applied ICEPOP to another dataset obtained from gut mucosal biopsy (GSE16879) derived from patients with inflammatory bowel disease (IBD)21. IBD is a chronic remitting and relapsing disease characterized by mucosal inflammation of the gastrointestinal tract. Tumour necrosis factor-alpha (TNF-α) is pivotal in the pathogenesis of IBD, therefore TNF-α blocking antibody therapy is utilized for IBD treatment22. The GSE16879 study involved three IBD diseases: ulcerative colitis (UC), Crohn’s ileitis (CDi), and Crohn’s colitis (CDc). All patients received TNF-α blocking monoclonal antibody infliximab (IFX). IFX responders (R) and non-responders (NR) were determined and reported previously21. Each patient had a biopsy before and after IFX treatment. Healthy control colon and ileum biopsy samples were included in this dataset. For calculation of DEGs, healthy colon control data was used for UC and CDc data (since they were colon biopsies), and healthy ileum control data was used for CDi data (since CDi samples were ileum biopsies) (Supplementary Fig. S1d).
In the ICEPOP analysis of UC patients, NRs showed almost identical results before and after IFX treatment (Fig. 5a vs b). This is also supported by the circular plot as many responsive genes were connected between “before” and “after”, indicating that their gene profiles were not changed with IFX treatment (Fig. 5c).
In contrast, the responses were noticeably different in UCR patients (Fig. 5d–f), where a notable decrease of macrophage activation was observed after IFX treatment (Fig. 5d vs e). A similar finding has been described23. The circular plot analysis further indicated that the number of listed genes was greatly reduced after treatment (Fig. 5f). Characteristically, TNFAIP624 was present in the pre-treatment but not post-treatment (Fig. 5f). This decrease of macrophage activation and the absence of TNFAIP6 after IFX treatment strongly suggest the effectiveness of IFX treatment in UC responders.
Concerning CDi, ICEPOP score of neutrophil considerably decreased after IFX treatment in CDi responders (Supplementary Fig. S7d vs S7e) but not in NR (Supplementary Fig. S7a vs S7b). Previous studies suggested the involvement of neutrophil activation in CD25,26. The circular map analysis further indicated that the gene responses in neutrophils, monocytes, dendritic cells, and macrophages were substantially decreased after IFX treatment in the CDiR group (Supplementary Fig. S7f) but not as much in the CDiNR group (Supplementary Fig. S7c). These tendencies were further supported by an alternative DEG calculation (see Supplementary Fig. S1d), in which data obtained prior to treatment was divided by data after treatment in the same group, namely healthy controls were not used for this alternative DEG calculation. Neutrophil and monocyte responses were preferentially decreased in CDiR (Supplementary Fig. S9b) but not in CDiNR (Supplementary Fig. S9a). Of note, in this analysis, increased scores in the bar graph (Supplementary Fig. S9b) and more listed genes in circular plot (Supplementary Fig. S9c) indicated that the actual responses decreased with IFX treatment.
ICDc data revealed a similar tendency to CDi in the alternative before-after treatment analysis. In the CDcR group, the neutrophil score increased (indicating decreased neutrophil response with IFX treatment; Supplementary Fig. S9e). This was not evident in the CDcNR group (Supplementary Fig. S9d). This was similarly demonstrated by the circular plot (Supplementary Fig. S9f; increased gene list means that the actual gene responses were decreased). Interestingly, the non-alternative bar graph analysis shown in Supplementary Figure S8 (this is the same analysis for CDi group as shown in Supplementary Figure S7) provided no obvious clue of neutrophil decrease in CDcR with IFX treatment (Supplementary Fig. S8d vs S8e). Instead, circular analysis of CDc (Supplementary Fig. S8f) already suggested the decrease of neutrophil, dendritic cell, and macrophage gene responses after IFX treatment.
Taken together, these results of CD biopsy indicate that the decrease in neutrophil responses in CD patients after IFX treatment may be a general feature of the IFX responsiveness and is consistent with a recent report that neutrophil surface marker status is associated with IFX treatment effectiveness27. To a lesser degree, there was a small decrease in monocytes in CDiR and CDcR (Supplementary Fig. S7 and Supplementary Fig. S8), which may implicate monocyte involvement in CD pathogenesis, consistent with that suggested by recent studies28,29.
The collective data indicate that ICEPOP analysis can provide practical and insightful information for a variety of mouse and human derived data.
To examine the robustness of ICEPOP analysis, we first examined whether the selection of cell type affects the ICEPOP score stability. In ImmGen, 214 subtypes are categorized in 10 cell types (Supplementary Fig. S10a). Similarly, in IRIS, 22 subtypes are categorized in 7 cell types (Supplementary Fig. S10b). For this examination, we explored three different cases to make an immune cell scoring matrix. Case 1 was a matrix of all individual subtypes. Case 2 was a matrix of every cell type containing with three randomly selected subtypes of each (Supplementary Fig. S11a,b). Case 3 was a matrix of 100 randomly selected subtypes irrespective of cell type, with the number of cell type potentially being fewer than 10 depending on the selected subtypes (Supplementary Fig. S11c,d). In case 2 and 3, 1000 sampling procedures were performed for ImmGen dataset, yielding 1000 scoring matrices for each case. This permutation analysis was not performed in IRIS dataset due to the small number of available subtypes.
ICEPOP analysis using the case 1 matrix provided finer sub-cell-type responses for the same mouse and human datasets used in this study (Supplementary Fig. S10). This fine matrix may be useful when focus on very specific immune subtypes is needed. ICEPOP analysis using case 2 and 3 matrices proved nearly consistent with the default matrix as a whole (Supplementary Fig. S11), suggesting that the results of ICEPOP analysis are stable to the random subtype selections, at least for these two different permutation approaches. In the GSE6332 dataset, neutrophils in bCD.ID remained the most responsive, whereas in bCD.IP, the immune response was almost flat (Supplementary Fig….