Generating drug chemogenomic profiles

The Connectivity Map (CMap) is a data repository of genome-wide transcriptional expression profiles collected from 6,100 experimental conditions of 1,309 unique compounds applied to human cell lines²³. Each perturbation is represented by a list of differentially expressed genes that we ranked based on fold-change. To capture the consensus profile of a compound across conditions, we merged multiple experiments (i.e. different drug concentrations or cell lines) for the same compound into a single Prototype Ranked List (PRL), using a hierarchical majority-voting scheme^28,29 (Fig. 1a). The collection of PRLs created a comprehensive resource for developing a systematic screening tool to look for connections between drug perturbations and immunological states (Supplementary Table 3).

IP map construction

We created a system-wide interaction map between drugs and immune cells by matching the 1,309 drug perturbation profiles in CMap to the 304 immune cell state changes we curated from the ImmGen compendia. Our matching algorithm evaluates the similarity between two transcriptional expression patterns by comparing the top and bottom ranked genes from both profiles^20,21 (Fig. 1b, Supplementary Fig. 4, and methods). Specifically, we tested the similarity between the immunological state change profiles (state B vs. state A) to each of the drug perturbation profiles (treated vs. untreated) by computing an immunemod score based on the overlap of the top- and bottom-ranked genes in each profile. A positive immunemod score indicates the specific drug treatment profile is similar to immune cell state B and suggests the drug biases the immune cell toward state B, whereas a negative score signals the drug treatment shifts the cell toward state A (Fig. 1b, d).

To evaluate the significance of our predicted drug-cell interactions, we generated random drug perturbation profiles for each compound and repeated the analysis 1,000 times for each immune cell state change (Supplementary Fig. 4a and methods). A complete computational integration of the CMap and ImmGen data sets produced 397,936 potential connections between drugs and immune cell state changes. To assess whether a predicted connection was robust, we varied the set size of the top- and bottom-ranked genes used for the matching algorithm and recalculated all ~400,000 immunemod score-p-value pairs. The proportion each drug-cell interaction was significant amongst all gene set sizes provided a relative weight for each predicted interaction (Supplementary Fig. 4b, c). Larger weights indicate a given drug-cell interaction depends less on the set size chosen to calculate the immunemod score and signifies a robust connection. This selection process enabled discovery of previously unknown interactions while minimizing spurious connections (Supplementary Fig. 5).

Using the significant and robust interactions, we made connections between drugs and immune cell state changes to generate a comprehensive IP map. The IP map contains 69,995 connections (Supplementary Table 4) that are significant at an FDR less than 5% and that appear in > 85% of gene set sizes. Although every drug showed a significant association with at least one of the 304 immunological state changes, the most frequent number of state transitions is 26, and 144 drugs influence 100 or more state changes (Supplementary Fig. 6). Drugs predicted to influence the largest number of immune cell state transitions include potent immunomodulators, many of which induce significant immunosuppression (Table 1). Drugs with immunomodulatory activity (for example, anti-inflammatory agents, anti-histamines, and immunosuppresants) show a significant enrichment for immune cell interactions (E = 1.5, P = 0.002, E = 1.4, P = 0.04, E = 2.1, P = 0.02, respectively, and Supplementary Table 5).

Global properties of the IP map

To examine the global landscape of the IP map, we used the immunemod score as a similarity metric and organized the complete set of drug and immune cell interactions through unsupervised hierarchical clustering (Fig. 2a). We found that drugs with similar therapeutic classes cluster together. For example, anti-psychotics (clozapine, loxapine, haloperidol, and fluphenazine) formed a cluster, as did purine analogs (mercaptopurine and tioguanine), and calcium channel blockers (dexverapamil, bepridil, and perhexiline). These three clusters are predicted to interact with the largest number of immune cell subset transitions. Drug clusters also showed enrichment for the same molecular target. For example, the anti-diabetic drugs troglitazone and rosiglitazone both target PPARG and ACSL4 as part of their mechanism of action for reducing blood glucose. Based on their immunemod scores, these drugs are predicted to influence naïve CD4⁺ T cells and NK cells, providing a potential explanation for the therapeutic benefits observed in patients with autoimmune disease^32,33.

Open in a separate window

Figure 2

Overview of the IP map

(a) Clustered heatmap of the drugs and immune cell state changes organized by the immunemod score. Venn diagrams reflect the number of drugs that shift immune cell states a specific or non-specific direction. (b) Two-dimensional heatmap of the stable cell and drug clusters with a significant enrichment for at least one cell type (rows) or a level 1 therapeutic class (columns). Therepeutic classes include: H = Systemic hormonal preparations, excl. sex hormones and insulins, V = Various, B = Blood and blood forming organs, P = Anti-parasitic products, M = Musculo-skeletal system, L = Anti-neoplastic and immunomodulating agents, G = Genito urinary system and sex hormones, R = Respiratory system, A = Alimentary tract and metabolism, D = Dermatologicals, J = Anti-infectives for systemic use, S = Sensory organs, N = Nervous system, and C = Cardiovascular system. Areas of corresponding circles represent the number of cells or drug per cluster (2–14 cells, 2–13 drugs). Gray colored squares indicate predicted interaction between at least one cell-drug pair in the pairs of clusters.

We identified 28 drugs associated with later stages of lineage development across multiple cell types (for example, stem/progenitor and pre-B cells; Fig. 2a, Supplementary Table 6). These drugs include compounds used to treat diseases of metabolism and the nervous, musculo-skeletal and respiratory systems, or are anti-infectives. Moreover, these 28 drugs are enriched for (i) agents with immunosuppressant, anti-psoriatic, and dopaminergic activity, and (ii) compounds that target the chromatin-associated enzyme PARP1, which is a key regulatory molecule for differentiation and proliferation^34,35. By contrast, 17 compounds influence immature hematopoietic cell subsets (Fig. 2a, Supplementary Table 7). These compounds aren’t associated with a single therapeutic class, yet their molecular targets are enriched for processes associated with oxidoreductase activity and alkylation repair, both of which are important for differentiation and maintaining stem cell integrity^36,37.

To further characterize the drug-immune cell interactome, we performed unsupervised hierarchical clustering with multiscale bootstrap resampling³⁸. Out of the 143 cell type changes, 119 fit into one of 25 stable cell clusters (P < 0.05, multiscale bootstrap analysis) (Supplementary Table 8). Almost half (13 / 25) of the stable cell clusters exhibited a significant over enrichment (E > 2, P < 0.05) for one or more cell types (Supplementary Table 9). By comparison, 1,089 drugs out of the 1,309 in total fell into one of 409 stable drug clusters (P < 0.05, multiscale bootstrap analysis) (Supplementary Table 10). Almost 87% (356 / 409) of the stable drug clusters showed a significant over enrichment (E > 4, P < 0.05) for one or more therapeutic class (ATC classification levels 1–3) or a molecular target, Supplementary Table 11). For example, 88 drug clusters showed significant enrichment for at least one Anatomical Therapeutic Chemical (ATC) Classification System level 1 description. The enrichments are driven by an abundance of anti-thrombotic agents or vitamin K and other hemostatics (B), contrast agents or diagnostic radiopharmaceuticals (V), and alkylating agents, cytotoxic ant-biotics, hormone antagonsits, or immunosuppressants (L) (Supplementary Fig. 7).

To examine the features of stable clusters in greater detail, we identified 53 drug clusters enriched for a therapeutic class and molecular target, and intersected these clusters with the 13 cell clusters enriched for one or more cell type (Fig. 2b). The intersecting clusters revealed that specific immune cell subsets (for example, pre-lymphocytes, monocytes, NKT cells, and gamma-delta T cells) overlap with multiple drug categories, whereas other subsets (for example, B cells and macrophages) intersect with a couple categories. A few drug clusters (for example, cl195, cl335) influenced multiple cell types, whereas other clusters (for example, cl2, cl48, cl49, cl97) influenced a single cell type. The drug cluster with the greatest overlap across immune cells (cl195) was enriched for anti-neoplastic drugs that are cytotoxic antibiotics (for example, doxorubicin and mitoxantrone).

We discovered a strong positive association between the number of molecular targets for a given drug and the number of interactions predicted to influence immunological state transitions (P = 1.7 × 10⁻⁵, linear regression). When we examined adverse drug interactions using the Side Effect Resource (version 2)¹, we found no relationship between the reported number of side effects for a drug and the number of immune cell interactions (P = 0.8, linear regression). However, since side effect data have a broad frequency distribution and are difficult to measure accurately², this lack of correlation may reflect the variation inherent in the bias of capturing and reporting side effects.

Connections and sub-structures in the IP map

To examine possible immunological outcomes that might result from connections in the IP map, we focused on immune cell state changes between cell subsets or tissues. Drugs and cell state transitions were organized by therapeutic class and category of immune cell state transition (e.g., within a subset or between tissue) to provide an overview of all the predicted connections that are statistically significant (Fig. 3a, Supplementary Fig. 8). We identified promiscuous drugs by their interactions with a large number of subset transitions (Table 2). Given the uneven cell type distribution within the subset state changes (Supplementary Table 1), we defined drug hubs based on interactions with the greatest number of cell types, which we hypothesize could have the largest influence on the immune system. By contrast, drug islands were defined on the basis of interacting with the fewest number of cells (Fig. 3b).

Open in a separate window

Figure 3

Predicting immune cell influence from the IP map

(a) Circular layout of the predicted connections between drugs and cell types (all connections statistically significant P < 0.05 with FDR < 5%, see methods). Line widths correspond to the number of interactions. The diagram has been organized by sorting the cell types counter-clockwise (drugs clockwise) in order of decreasing number of connections. Single letter codes for each drug follows the anatomical therapeutic classification system as described in Fig. 2. (b) Subnetwork showing the drug hubs and islands (square nodes in the center and periphery respectively), and their predicted interactions with immune cell subsets (circles). Each square node represents a single drug (for example, the four center nodes labeled L represent, counter clockwise from top, the drugs daunorubicin, azacitidine, vorinostat, and methotrexate). (c) Drugs identified in the IP map based on their immunemod score and predicted to increase neutrophil and monocyte frequencies in the blood. (d) Validation of drug prediction in immune cells from humans. Drug influence (+ = drug, − = no drug) on measured cell counts from lab values of patients in the Mount Sinai Electronic Medical Records (cell counts and frequencies). Bar heights represent median count (neutrophil) or frequency (monocyte) values and bottom/top of the error bars reflect the inter-quartile ranges (25% to 75% of data respectively). Significance assessed via Wilcoxon rank sum test.

Hubs were enriched for anti-neoplastic drugs (E = 13.8, P < 1 × 10⁻⁵) (Supplementary Fig. 9), which could be expected given the influence these compounds have on immune cells^1,2. Hubs were also enriched for nervous system compounds such as the selective serotonin reuptake inhibitor zimeldine (E = 1.3, P < 0.04), which was pulled off the market due to a rare, but severe adverse reaction leading to the autoimmune condition known as Guillain-Barre syndrome³⁹, and the anti-seizure drug topiramate, which was shown to be an effective treatment for inflammatory bowel disease in a pre-clinical model²¹ (Supplementary Fig. 10a). By contrast, drug islands were enriched for metabolic drugs that included the anti-diabetic compounds gliquidone and metformin (Supplementary Fig. 10b). This metabolic island we identified in the IP map mirrors the low connectivity found in metabolic diseases in the human disease network⁴⁰.

To explore drugs predicted to influence a portion of the adaptive immune system, we identified a subnetwork based on the largest magnitude immunemod scores for T cell subsets and tissues (Supplementary Fig. 11). This subnetwork includes more than one thousand compounds predicted to influence CD4⁺ or CD8⁺ subsets, with 113 and 202 compounds unique to each group respectively. The top immunemod score for this subnetwork is between CD4⁺FoxP3⁺ T cells and guanfacine (Supplementary Figs. 12, 13a). Guanfacine is an α2A receptor (ADRA2A) agonist used for lowering blood pressure and treating ADHD⁴¹.

To verify that the immunemod score identifies a drug’s influence on a specific cell subset, we administered the anti-hypertensive drug guanfacine to mice and measured the percentages of regulatory T cell subsets isolated from spleens. Based on the immunemod score direction, we reasoned the frequency of regulatory T cells should increase following treatment with guanfacine. In comparison to untreated mice, the treated mice showed 5% increase in the average frequency of CD62L⁺ cells within the CD4⁺FoxP3⁺ T cell compartment (42.0% vs. 37.0%: n = 15 treated vs. n = 14 untreated, P = 0.01, Supplementary Fig. 13b).

When we examined the CD8⁺ subsets, the top immunemod score was to the anti-parkinsonian drug trihexyphenidyl. One molecular target for trihexyphenidyl is the muscarinic acetylcholine receptor M1 (CHRM1). When the gene that encodes for this molecular target is knocked out in a mouse model, CD8⁺ cells from these mice exhibit defective cytotoxic capability⁴².

Validation of drug-immune cell interactions in humans

Immune cell data were collected from mice and drug perturbation data were gathered from human cancer cell lines. One concern with integrating these data is whether the findings from mice translate to humans. We tested whether interactions predicted by the IP map could be observed in humans by examining immune cell counts of patients administered drug versus untreated patients. To compare patient populations, we examined complete blood counts for more than 2.3 million electronic medical records in the Mount Sinai Hospital System and selected individuals who were treated and had blood cell counts collected within one month of receiving drug.

Given the constraints of routine clinical lab tests found in electronic medical records, we restricted IP map predictions to two common drugs predicted to influence monocytes and neutrophils. The IP map predicted the general anesthetic propofol and the anti-hypertensive spironolactone, would increase neutrophils and monocytes respectively (Fig. 3c). Propofol increased neutrophil counts by 2,500 cells / mm³, and spironolactone increased monocyte frequencies 1.6% (Fig. 3d). Although the cell population changes were small, both shifts were significant (P < 1 × 10⁻¹⁰⁰, Wilcoxon rank sum). Furthermore, the neutrophil increase shifted most patients beyond the upper normal range. To validate these observations, we examined the same drug-immune cell pairings in the electronic medical records of Columbia University Medical Center. This independent data source showed the same direction and significance for both drugs and their predicted influence on immune cells (Supplementary Table 12).

IP map construction

We constructed a matrix of predicted interactions between each of the 1,309 drugs and 304 immunological state changes using a rank-based, pattern-matching strategy described previously²⁰. Briefly, for each trio of drug, cell, and gene set size (d, c, s), we calculated an Immunemod Score (ImS) based on the degree of overlap between drug and immune cell gene sets at the extremes of the two ranked signatures. To obtain a measure of significance for the immunemod score, we shuffled the genes in the drug rank signature and calculated a permuted Immunemod Score (ImS*) for each drug, cell, and gene set triplet [ImS*(d_i, c_j, s_k)]. We calculated the p-value for each ImS by counting the number of randomized scores ImS*(d_i, c_j, s_k) that were greater than or equal to the absolute value of the actual scores ImS(d_i, c_j, s_k) and dividing by the number of permutations (n_perms = 1,000). This permutation strategy sets the lower bound for p-values at 0.001, which yields a biased estimate for the number of false positives given the number of hypotheses under consideration. To provide accurate p-values at the lower range while containing the computational cost, we used the generalized Pareto distribution to model the p-value distribution and calculated improved estimates for low p-values (counts < 1 / 100) based on the distribution of permuted immunemod scores³⁰. We adjusted the p-values⁶⁵ and selected an FDR of 5% as the cutoff for significance. To control for spurious interactions based on the size of the gene set used for matching, we varied the size of the matching set between 100 and 250 genes for each of the top and bottom extremes and recalculated all immunemod score, p-value pairs for every drug-cell interaction. The proportion of times each drug-cell interaction was significant amongst all sizes of gene sets provides a relative strength for any given interaction. A predicted interaction was considered to be strong and stable if it was significant for 85% or more of the set sizes.

Data analysis

To assess the similarity between expression profiles of immune cell subsets, we calculated the Jaccard distance amongst all pairs of extreme fold-change genes, and used an ANOVA to evaluate the differences between immune cell subsets. We investigated a series of diagnostic plots and did not find significant deviations that would violate the assumptions of normality or homoscedasticity.

To organize the drugs and immune cells in an unbiased manner, we applied hierarchical clustering to the full interaction matrix using the computed Pearson correlation coefficient as a distance metric between immunemod scores and complete linkage clustering to agglomerate drugs or cells. We used the pvclust R package³⁸ to compute a bootstrap analysis of the clusters and identified a significant cluster if the approximately unbiased probability was > 95%.

To determine the enrichment of an anatomical therapeutic class category, we calculated the fold-change and p-value. Fold-change enrichment (E) was calculated as a ratio of ratios E = (a / b) / (c / d), where a is the number of drugs with a particular category (for example, “L”) in the cluster of interest, b is the number of drugs with that category in the overall data set, c is the total number of drugs in cluster, and d is the total number of drugs overall. We used the hypergeometric distribution to calculate the p-value and assess the significance of each enrichment calculation.

To examine the association between chemical features (for example, molecular targets and drug side effects) and number of immune cell interactions, we implemented a simple linear regression model. Chemical features followed a skewed distribution so we log-transformed the data, which adjusted the values so they followed a normal distribution. Based on diagnostic plots of the transformed data, we did not find deviations that violated the assumptions of normality and homoscedasticity that are central to regression models.

To test whether drug treatment with clioquinol or amantadine produced any difference in neutrophil cell frequencies in various tissue compartments, we used an ANOVA model to compare the treatment groups. Multiple group testing and p-values were evaluated using Tukey’s honest significant difference. For all ANOVA tests, we generated a series of diagnostic plots to examine: (i) the residual errors for outliers, (ii) the QQ plots for normality, and (iii) the square root of the standardized residuals for heteroscedasticity. In all cases, we did not find significant deviations that would violate the assumptions used in the ANOVA model. For comparison, we evaluated the treatment and control groups directly via the Wilcoxon rank sum test and found the median differences between treatment (Clioquinol, Amantadine) and controls (PEG400, PBS) to follow the exact same pattern obtained using the ANOVA, with a similar maximum p-value for significant differences (P ≤ 0.01). To test whether guanfacine influence regulatory T cell frequencies, we used a meta-analysis strategy⁶³ to combine experimental conditions and groups, which allowed us to ascertain whether the overall differences from each independent experiment were robust and significant.

Electronic medical records

We pulled all patient entries from the Mount Sinai Electronic Medical Records that contained complete blood count information on neutrophils and monocytes (more than 2.3 million entries in total). To determine if either propofol or spironolactone were associated with a change in cell counts, we identified patient entries that had lab values measured within 30 days of drug administration versus patient entries that never received drug. We tested for group-level differences using the non-parametric Wilcoxon rank sum test. The findings were validated using the Electronic Medical Records of Columbia University Medical Center, where we employed the same criteria for patient selection and performed a Wilcoxon rank sum test for group differences.

Visualization

Circos plots created using the circlize R package (version 0.0.7 https://github.com/jokergoo/circlize). Network diagrams produced using Cytoscape⁶⁶ and SPADE trees generated with CytoSPADE⁶⁷. All other plots created using the R statistical package.

Mice and drug treatment

6–12 week old female C57Bl/c mice were obtained from Jackson Laboratories. Mice received intraperitoneal injections 3 times every 12 h with clioquinol, amantadine (both at 30 mg/kg per dose, Sigma Aldrich) or appropriate controls. Dosing level and frequency were chosen based on previous experiments using clioquinol in mice^47,68 and the drug half-life (11–14 h). Clioquinol was dissolved in 8% PEG400/PBS heated to 37 °C; amantadine was dissolved in PBS. Before injection, the solutions were shaken several times. Mice were sacrificed for tissue collection between 2.5 and 3 h after the last treatment. Blood collection was obtained from the tail. For guanfacine treatment, mice received initial injection of 5 mg/kg of drug or vehicle control (PBS) on day 1, followed by two (experiment 1) or six (experiments 2–5) intraperitoneal injections at 2 mg/kg every 12 h starting on day 2. Mice were sacrificed for tissue collection 12 h after the last treatment. All animal procedures were done according to protocols approved by the Mount Sinai School of Medicine Institutional Animal Care and Use Committee.

Flow cytometry

Peritoneal cavity cells were collected by washing with cold PBS containing 4% FBS. Single-cell suspensions of bone marrow were obtained by flushing femurs, followed by filtration through a 100-µm cell strainer (BD Biosciences). Red blood cells were lysed for 2 min at room temperature with RBC lysis buffer (eBioscience). Samples were stained with the following antibodies (all from eBioscience): allophycocyanin-eFluor780-conjugated CD45 (30-F11), peridinin chlorophyll protein_cyanine 5.5-conjugated CD11b (M1/70), phycoerythrin-conjugated Ly6G (RE6-8C5), phycoerythrin-conjugated CD3 (145-2C11), peridinin chlorophyll protein–cyanine 5.5-conjugated CD25 (PC61.5), fluorescein isothiocyanate-conjugated CD62L (MEL-14), allophycocyanin-conjugated FoxP3 (FJK-16s), eFluor450-conjugated CD4 (GK1.5), allophycocyanin-eFluor780-conjugated CD8a (53–6.7), and allophycocyanin-conjugated CD44 (clone IM7, BD Pharmingen). DAPI was used to label dead cells. LSR Fortessa was used for sample acquisition and FlowJo software for data analysis.

RNA isolation and quantitative PCR

Total RNA was extracted from pieces of lung, liver, spleen and bone marrow cells using QIAzol Lysis Reagent (Qiagen) and glycogen blue (Ambion, Life Technologies) according to the manufacturer’s instruction. For cDNA synthesis, 2 µg total RNA was reverse-transcribed for 1 h at 37 °C with an RNA-to-cDNA kit (Applied Biosystems). For quantitative PCR, SYBR green qPCR master mix 2° (Fermentas, Thermo Scientific) and the following primers were used: mouse Actb forward, 5′-CTAAGGCCAACCGTGAAAAG-3′, and reverse, 5′-ACCAGAGGCATACAGGGACA-3′; mouse Cxcl1 forward, 5′-GTGTTGCCCTCAGGGCC-3′, and reverse, 5′-GCCTCGCGACCATTCTTG-3′; mouse Cxcl2 forward, 5′-ACGCCCCCAGGACCC-3′, and reverse, 5′-CTTTTTGACCGCCCTTGAGA-3′; mouse Cxcl5 forward, 5′-CTCGCCATTCATGCGGAT-3′, and reverse, 5′-CTTCAGCTAGATGCTGCGGC-3′; mouse Cxcr2 forward, 5′-CTTTGCCCTGACCTTGCCT-3′, and reverse, 5′-GCACAGGGTTGAGCCAAAA-3′; mouse Cxcr4 forward, 5′-TGGCCTTCATCAGCCTGG-3′, and reverse, 5′-TTGGCCTTTGACTGTTGGTG-3′.

We thank G. Hoffman and B. Readhead for useful conversations about the computational methods and suggestions on the manuscript, L. Li for assistance with the electronic medical records from Mount Sinai, K. Oguntuyo for assistance with the qPCR reactions, and A. Rahman for sample preparation and acquisition on the CyTOF. This work was supported by the US National Institutes of Health 1R01AI104848, 1R33CA182377, and 1R01AI113221 (to B.D.B.), U54CA189201 and R01DK098242, and Jonathan E. Gray TCI Young Scientist Cancer Research Award (to J.T.D.).