Defining tissue-and disease-associated macrophages using a transcriptome-based classification model

Macrophages are heterogeneous multifunctional leukocytes which are regulated in a tissue-and disease-specific context. Many different studies have been published using in vitro macrophage models to study disease. Here, we aggregated public expression data to define consensus expression profiles for eight commonly-used in vitro macrophage models. Altogether, we observed well-known but also novel markers for different macrophage subtypes. Using these data we subsequently built the classifier macIDR, capable of distinguishing macrophage subsets with high accuracy (>0.95). This classifier was subsequently applied to transcriptional profiles of tissue-isolated and disease-associated macrophages to specifically define macrophage characteristics in vivo. Classification of these in vivo macrophages showed that alveolar macrophages displayed high resemblance to interleukin-10 activated macrophages, whereas macrophages from patients with chronic obstructive pulmonary disease patients displayed a drop in interferon-γ signature. Adipose tissue-derived macrophages were classified as unstimulated macrophages, but resembled LPS-activated macrophages more in diabetic-obese patients. Finally, rheumatoid arthritic synovial macrophages showed characteristics of both interleukin-10 or interferon-γ signatures. Altogether, our results suggest that macIDR is capable of identifying macrophage-specific changes as a result of tissue-and disease-specific stimuli and thereby can be used to better define and model populations of macrophages that contribute to disease.


Introduction 50
Macrophages are multifunctional innate immune cells that play a central role in the spatiotemporal 51 regulation of tissue homeostasis between pro-inflammatory defense and anti-inflammatory tissue repair. 52 Dysregulation of macrophages has been implicated in a variety of disorders. As in vivo macrophages are 53 often difficult to obtain and study, in vitro peripheral blood monocyte derived macrophages (MDMs) have 54 been used extensively as model systems for assessing the transcriptional and functional regulation in 55 response to various stimuli. 56 To mimic in vivo macrophages encountering different microenvironmental signals (1, 2), MDMs, 57 differentiated with for example macrophage-(M-CSF), or granulocyte-macrophage stimulating factor (GM-58 CSF), are activated in vitro with bacterial lipopolysaccharides (LPS) or Th1 cytokine interferon gamma 59 (IFNγ) to generate pro-inflammatory macrophages (M1). Anti-inflammatory macrophages (M2) are often 60 generated by activating the cells with Th2 cytokines, interleukin-4 (IL4), or other anti-inflammatory stimuli, 61 such as interleukin-10 (IL10), tumor growth factor (TGF) and glucocorticoids (1-6). By applying these pro-62 or anti-inflammatory stimuli to the MDMs, researchers sought to obtain proper models for studying the 63 A multinomial elastic net classifier distinguishes macrophage activation states with high accuracy 126 Subsequently, we sought to perform feature selection to identify genes capable of distinguishing 127 macrophage activation states. We merged the data from different microarray and RNA-seq platforms and 128 performed multinomial elastic net regression on the 5986 overlapping genes present in all datasets. The 129 expression data was randomly split into a training (2/3) and a test (1/3) set whereupon the training set was 130 used to build a model through ten-fold cross-validation. We repeated the cross-validation procedure 500 131 times and took the median thereof to stabilize the log odds ratios (35). Genes were considered stable 132 predictors if their log odds ratio was non-zero in more than 50% of the iterations (Fig. S1). Subsequent 133 classification was performed using the median log odds ratio per gene across all iterations ( Fig. 2A). 134 Altogether, our classifier was composed of 97 median-stabilized predictor genes, and was compiled as an 135 R package called macIDR (https://github.com/ND91/macIDR). 136 To validate our model, we tested macIDR against the previously withheld test set, which included 137 a newly-generated RNA-seq experiment containing all included activation states. Classification of the test 138 set revealed an accuracy above 0.95 with both high sensitivity (>0.98) and specificity (> 0.83; Table 2). In 139 total, 75 out of 79 test samples were correctly classified (Fig. 2B). Notably, for three of the four 140 misclassified samples, the second-best prediction was the subset as reported by the authors. Investigation  Table S3). A follow-up 148 transcription factor motif analysis on the promoters of the predictor genes showed significant enrichment 149 for macrophage transcription factors, the E26 transformation-specific PU.1 (Spi1) and SpiB (Fig. 2E). 150 MacIDR generalizes to granulocyte-macrophage colony-stimulating-factor differentiated 151 macrophages 152 Besides M-CSF, granulocyte-macrophage colony-stimulating-factor (GM-CSF) is often used to differentiate 153 monocytes to macrophages, where it is thought to evoke a more pro-inflammatory phenotype (36-38). We 154 investigated whether macIDR was capable of discerning GM-CSF-differentiated MDMs (GM-MDMs) 155 exposed to various stimuli. In total, we obtained 31 datasets from three microarray studies (5, 16, 39) 156 where GM-MDMs were activated with LPS (GM-LPS early , GM-LPS late ), IFNγ (GM-IFNγ), and IL4 (GM-IL4), 157 or remained unstimulated (GM0). These studies were selected based on their similarity in activation 158 duration with the M-CSF differentiated macrophages used in the training set. We observed that all  MDMs were classified as their M-CSF counterparts (Fig. 3B) despite the fact that some predictive genes 160 were absent from two studies (Table S4) (Table S4). By contrast, we found that the FLS, T, B, NP and NK 174 cells were all classified as M0 (Fig. 3A), suggesting that the M0 class is used as a label for expression 175 profiles of non-monocyte-derived cells. By looking at the distribution of the log odds we observed that true 176 M0 classifications generally displayed higher log odds than the M0 classifications of most non-macrophage 177 cells, but found this not to be definite (Fig. 3D). The increased variance for the M0 signal observed for GM-178 MDMs and MoDCs was due to the different GM-MDMs activation and DC maturation states, as GM0 and 179 up to some extent DC0, depicted M0 log odds similar to true M0 samples. 180 Alveolar macrophages from COPD and smoking individuals show a reduced M-IFNγ signal

181
We next attempted to classify macrophages derived from patient tissues to study their semblance to in 182 vitro generated MDMs. We investigated alveolar macrophages (AMs) obtained through bronchioalveolar 183 lavage from smoking individuals, COPD patients, asthma patients (31), and healthy control (HCs). Overall, 184 we found that the AMs were classified primarily as M-IL10 (Fig. 4A and Table S10), which appears to be 185 driven by the MARCO signal. This corroborates the observation that lung tissue, specifically AMs, display 186 significantly higher gene expression of MARCO relative to surrounding cells and other tissues (41, 42). 187 Moreover, MARCO was found to be necessary in AMs for mounting a proper defense response (41-43). 188 Further comparisons of the AMs derived from different patient groups, we found that the macrophage 189 signal could be stratified according to health status where COPD-and smoker-derived AMs displayed a Furthermore, while no clear differences were observed when comparing AMs from asthma 196 patients to HCs, AMs from steroid-sensitive asthma patients displayed a decreased M-dex signal and 197 increased M-LPS early , M-IFNγ and M-LPS+IFNγ signal relative to steroid resistant asthma AMs (Fig. 4A). 198 The apparent difference in steroid sensitivity appeared to be caused by a difference in TNF and CXCL9 199 signals contributing to the M-IFNγ and M-LPS+IFNγ classification respectively (Fig. S2). This observation 200 corroborates previous studies where IFNγ signaling was found to suppress glucocorticoid-triggered 201 transcriptional remodeling in macrophages leading to the macrophage-dependent steroid-resistance, 202 thereby reflecting a higher level of IFNγ in steroid-resistant relative to steroid-sensitive asthma patients. 203 Adipose macrophages show a M0 and M-IL4 classification

204
We next investigated visceral adipose tissue macrophages (ATM) derived from diabetic and non-diabetic 205 obese patients (34). Classification analysis suggested visceral ATMs showed most similarity with M0 206 followed by M-IL4 (Fig. 4B). While we were not capable of defining a set of genes responsible for the M0 207 classification, we observed that the concordant expression of MAOA and C-C chemokine ligand 18 208 (CCL18) contributed the most to the M-IL4 signal (Fig. S3). MAOA encodes a norepinephrine degradation 209 enzyme and is expressed more in sympathetic neuron-associated macrophages isolated from the 210 subcutaneous adipose tissue. In these cells, MAOA's norepinephrine clearance activity has been linked to 211 obesity (44). Interestingly, when comparing ATMs from diabetic obese with non-diabetic obese patients, 212 we observed a stronger signal of M-IL10 and M-LPS early driven by CCL18 and TNF respectively (Fig. S3). 213 Notably, CCL18 expression in both visceral and subcutaneous adipose tissue has been associated with 214

237
In this study, we performed a macrophage characterization study by integrating public datasets of eight in 238 vitro macrophage activation states. Our meta-analysis returned both well-known and novel markers for 239 activated macrophages. At a genome-wide level, we observed separation according to the conventional 240 pro-and anti-inflammatory macrophages. We subsequently built a classification model capable of 241 discriminating macrophage activation states based on their transcriptomic profile and made this available 242 as an R package called macIDR. By applying macIDR to in vivo macrophages, we projected the latter onto Similarly, immature MoDCs (DC0) displayed a slightly higher response for M0. All other non-macrophage 267 cells were classified as M0. While the log odds of proper M0 classifications were slightly higher than 268 improper M0 classifications, we acknowledge that no clear threshold could be defined. We are unsure why 269 M0 was predicted as a class for non-MDM and non-MoDC datasets, but we speculate that it might be 270 related to the M0 class being represented by most predictor genes relative to the other classes. We 271 therefore recommend potential users of macIDR to determine a priori that their dataset of interest 272 represents macrophages either experimentally or using in silico cell-composition estimation methods, such 273 as CIBERSORT (54) or xCell (55). SMs suggested that they represent multiple subpopulations, such as the CD163 + anti-inflammatory tissue-293 resident macrophages and the S100A8/9 + pro-inflammatory macrophages recruited from peripheral 294 monocytes (56). 295 As in vivo macrophages likely represent a more heterogeneous population compared to the in vitro 296 MDMs used in building the classifier, prediction of macrophage mixtures using bulk RNA-seq or microarray 297 returns signals from multiple different subsets. The classification therefore mainly reflects altered subset 298 proportions due to disease progression or association. It is likely that some in vivo macrophage 299 subpopulations do not share transcriptomic signatures with any of the eight in vitro MDMs preventing the 300 exact characterization thereof. Future studies using single cell RNA-sequencing should aim at defining 301 novel additional in vivo macrophage subsets. Nonetheless, we were capable of extracting signals from the 302 in vivo macrophage classifications that not only corroborated previous findings, but also provided novel 303 features for future research. It is essential to analyze macrophage subsets within the context of their in 304 vivo environment and therefore we provide this quantitative method to aid researchers in better defining 305 and modelling macrophages in tissue and disease. The initial screen was limited to studies that investigated primary macrophages and excluded the 315 stem cell derived macrophages or immortalized cell lines. Then we categorized macrophage subsets 316 based on the stimuli and the treatment time. For each subset, we sought to obtain at least 4 studies 317 including at least 2 biological replicates. Further, as a background control, only studies including 318 unstimulated control macrophages were selected. After this screening, we investigated macrophages 319 stimulated with control medium, LPS (either for 2 to 4 hours or 18 to 24 hours), LPS with IFNγ, IFNγ, IL4, 320 IL10 or dexamethasone for 18 to 24 hours. Microarray datasets generated on platforms other than 321 Illumina, Affymetrix or Agilent were excluded to ensure comparability. As we only investigated genes that 322 were measured in every single study, datasets that displayed limited overlap in the measured genes with 323 the other studies were removed. The in vivo macrophage datasets were samples obtained from clinical 324 specimen. In total, we obtained 206 datasets belonging to 19 studies for the meta-analysis and 325 classification. 326

327
A random effects meta-analysis was performed on the normalized data using the GeneMeta (v1.52.0) (57) 328 package, which implements the statistical framework outlined by Choi et al. (25). In short, the standardized 329 effect size (Cohen d adjusted using Hedges and Olkin's bias factor) and the associated variance were 330 calculated for each study by comparing each activated macrophage with the unstimulated macrophage 331 within each study. The standardized effect sizes were then compared across studies by means of random 332 effects model to correct for the inter-study variation, thereby yielding a weighted least squares estimator of 333 the effect size and its associated variance. The estimator of the effect size was then used to calculate the 334 Z-statistic and the p-value, which was corrected for multiple testing using the Benjamini-Hochberg 335 procedure. We modified the GeneMeta functions to incorporate the shrunken sample variances obtained 336 from limma (58) for calculating the standardized effect sizes. 337

338
Raw microarray and RNA-seq data were log 2 transformed where necessary after which the data was 339 (inner) merged and randomly divided into a training (2/3) and test set (1/3). As the test set should remain 340 hidden from the training set, raw microarray and RNA-seq data was used instead of normalized data to 341 prevent data leakage (59). Elastic net regression was subsequently performed using the R glmnet (v2.0) 342 (60) package. As penalized regression approaches are sensitive to the magnitude of each feature, we 343 investigated the use of standardization. As the lowest deviance appeared to be higher in the inter-fold 344 standardized training set relative to raw data, we did not include any standardization (Fig. S1). We also 345 investigated the optimal alpha for minimizing the deviance by means of an initial grid search approach, 346 which was found to be 0.8. 347 The training set was subjected to ten-fold cross-validation for tuning the penalty regularization 348 penalty parameter lambda. This process was subsequently repeated 500 times to stabilize the 349 Page 11 of 23 randomness introduced during the splitting step for cross-validation (35). We considered genes to be 350 stable classifiers if they displayed a non-zero log odds ratio in at least 50% of the 500 iterations. The final 351 log odds ratio was selected by taking the median of each stable predictor gene across the 500 iterations. 352 We subsequently validated our classifier on the withheld test set. 353 In addition to the training and test data, we downloaded and imported additional datasets from 354 GM-MDMs, non-macrophage cells, and in vivo macrophages. Subsequent classification was performed 355 using the macIDR package. Unlike the studies included for training and testing, some of the included 356 studies were performed on platforms that did not measure the expression of some of the predictor genes. 357 To that end, the relative log odds ratios were calculated, which represent the log odds ratio present relative 358 to the total log odds ratio had all predictor genes been present.       genes ranked by p-value. Columns represent the consensus sequence, the motif name, the p-value, the 572 median-stabilized log(OR) EBI3 CCL18 MAOA  IDO1  TNFAIP6  IL7R  TNFRSF21  CD9  CLEC10A  CCL22  CTNNAL1  G0S2  SLC2A5  NAMPT  ENPP2  RPL37  S100A9  PLTP  SOCS3  RPL23A  THBD  PTX3  TNF  CYTH4  CXCL10  EXOG  MX2  SERPINA1  STAB1  FCMR  TIMP3  CD101  PLEK  H1F0  IL1B  OLFML2B  SEMA4D  HLA-DOA  CLIC1  ITGB7  RPL6  AVPI1  SGK1  IFI44L  RNASE1  JUND  SERPING1  BCL2A1  ASAP2  EMP1  MMP19  RGS1  CD52  BMP6  A2M  MS4A6A  CXCL1  SORL1  FCGBP  IER2  CCL5  MT1G  ALOX5AP  SLC39A8  AMIGO2  SLAMF1  TFPI  PALLD  GPR153  NCAPH  APBB1IP  FOS  RPL9  FKBP5  CD47  GCLM  DDIT4  RHOBTB1  ALOX15B  TPST1  FSCN1  RAB3IL1  USP49  JAKMIP2  CXCL9  VSIG4  ADORA3  HMGB2  HTR2B  TAF10  FUCA1  CAMK1  TNFSF15  GM2A  TREM2      Pairwise gene correlations were generated by calculating the Spearman correlation to minimize 640 the effect of outliers. Pathway and regulator analyses were performed using the Canonical Pathways and 641 Ingenuity Upstream Regulator found in the Ingenuity Pathway Analysis software package (QIAGEN Inc., 642 https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis) and Metascape 643 (http://metascape.org/gp/index.html). Transcription factor motif analysis was performed by using HOMER 644 (v4.10) using the following parameters: -start -200 -end 100 -len 8, 10, 12. Table S1 meta-analysis results. A table containing the unbiased estimator of the effect size and the 671 variance as calculated by the meta-analysis ranked by the p-value. Separate tabs contain the results for 672 the different comparisons of the meta-analysis. Columns "Entrez" and "Gene" represent the Entrez gene 673 ID and the HGNC symbol respectively. The column "mu" represents the unbiased estimator of the effect 674 size and the "mu_var" represents the unbiased estimator of the variance. The columns "Z", "Z_pval", and 675 "Z_pval"BH" represent the Z-statistic, the associated p-value and the Benjamini-Hochberg adjusted p-676 value. 677 678  INHBA  AQP9  TNPO1  CXCL9  HSD17B7  PPIB  MAN1C1  VKORC1  SHCBP1  RNASET2  S1PR4  USP13  DHRS3  CXCL10  A2M  MT1G  IFI44L  RNASE1  SORL1  CD101  TIMP3  STAB1  FCMR  TAF10  FUCA1  CAMK1  GM2A  TREM2  TNFSF15