Systematic investigation of cytokine signaling activity at the tissue and single-cell levels.

Cytokines are critical for intercellular communication in human health and disease, but the investigation of cytokine signaling activity has remained challenging due to the short half-lives of cytokines and the complexity/redundancy of cytokine functions. To address these challenges, we developed the Cytokine Signaling Analyzer (CytoSig; https://cytosig.ccr.cancer.gov/ ), providing both a database of target genes modulated by cytokines and a predictive model of cytokine signaling cascades from transcriptomic profiles. We collected 20,591 transcriptome profiles for human cytokine, chemokine and growth factor responses. This atlas of transcriptional patterns induced by cytokines enabled the reliable prediction of signaling activities in distinct cell populations in infectious diseases, chronic inflammation and cancer using bulk and single-cell transcriptomic data. CytoSig revealed previously unidentified roles of many cytokines, such as BMP6 as an anti-inflammatory factor, and identified candidate therapeutic targets in human inflammatory diseases, such as CXCL8 for severe coronavirus disease 2019.

Introduction

Cytokines are a broad category of intercellular signaling proteins that act in almost every aspect of human immunology, from anti-pathogen immune responses to tissue-damaging inflammation[1,2]. However, the precise characterization of cytokine signaling activities has proven difficult due to two vexing properties of cytokine activity: redundancy and pleiotropy. Many cytokines, especially those with similar cell surface receptors and downstream pathways, have cellular effects that appear redundant within a specific cellular context[3]. At the same time, cytokines often have pleiotropic functions within an organism that depend heavily on cell-type-specific receptor usage and the presence of other signaling components[3].

This apparent redundancy and pleiotropy in cytokine activities are poorly captured by most immunological assays such as the Enzyme-Linked Immuno-Sorbent Assay (ELISA) and Luminex xMAP, which directly measure cytokine release. Cytokine release can be transient, unlike the longer-lasting and more functionally relevant measurement of target activities[4]. Recognizing this limitation, researchers have attempted to create databases of cytokine signaling targets. For example, the “Interferome” identifies interferon target genes in humans and mice through the collection and analysis of microarray data[5]. Gene Set Enrichment Analysis (GSEA) also annotates response genes for selected cytokines based on prior knowledge[6]. However, these databases and approaches cover a small fraction of cytokines, leaving most cytokine-induced target changes unexplored.

The need for systematic profiling approaches that allow modeling of cytokine target activity is urgent because cytokines can trigger life-threatening symptoms in many diseases. For example, COVID-19 mortality has been attributed mainly to a virus-induced cytokine storm, defined by excessive production of pro-inflammatory cytokines that lead to acute respiratory distress and widespread tissue damage[7]. Although pro-inflammatory cytokines help activate the immune response, there does not appear to be a strong relationship between cytokine storm severity and pathogen clearance. For example, successfully recovering COVID-19 patients may not have any inflammatory symptoms[8]. Cytokine release syndrome also causes severe side effects in many cancer treatments, such as immunotherapies[9] and chimeric antigen receptor (CAR) T therapies[10]. Similar to the disconnect between the severity of immune-related symptoms and disease outcomes in COVID-19, complete tumor remission can occur in patients without cytokine release syndromes[11]. While the immunological mechanisms of these observations remain unclear, they imply that if properly modulated, the benefits of cytokine signaling can be realized without substantial pathological effects.

With this goal in mind, and to model cytokine activity generally, we developed Cytokine Signaling Analyzer (CytoSig, https://cytosig.ccr.cancer.gov), a data-driven infrastructure hosted by the National Cancer Institute (NCI). CytoSig includes both a database of target genes modulated by cytokines and a predictive model of cytokine signaling activity and regulatory cascade from transcriptomic profiles. To build the CytoSig platform, we first created the Framework for Data Curation (FDC) to assist expert annotations on metadata deposited on databases through natural language processing functions (https://curate.ccr.cancer.gov). Using FDC, we analyzed 9,271 published studies and curated 20,591 transcriptomic profiles for human cytokine, chemokine, and growth factor response to create the CytoSig database and predictive model. We validated CytoSig by showing that it can reliably predict cytokine target activities in both human clinical studies and our in-vivo experiments. Further, CytoSig identified CXCL8 signaling as potential COVID-19 therapeutic targets that may alleviate adverse inflammation without undermining protective immunity.

Results

The Framework for Data Curation on Public Repositories

We hypothesized that the large number of cytokine-treatment datasets available publicly could serve as a knowledge base to model signaling activities in diverse biological contexts. However, two hurdles must be overcome to transform this body of data into a useful model. First, the experimental design behind each published dataset is unique, requiring labor-intensive expert interpretation of the metadata and standardization of the data into a format suitable for automated analysis. Second, one must identify and exclude experiments that involve cell models, stimuli, doses, or time intervals that are not physiologically relevant. More broadly, such challenges exist for many other biological topics that could be addressed by data aggregation. To overcome these hurdles, we established the Framework for Data Curation (FDC), which couples large-scale automatic data processing with natural language processing functions to assist expert annotation of experimental design (Online Methods, Fig. 1a).

Fig. 1:

Curation of Human Cytokine Response Data

a, The Framework for Data Curation (FDC). The FDC can automatically process RNA-Seq and MicroArray transcriptomic data from public repositories. Then, with FDC’s natural language processing functions, curators read the metadata of each sample to annotate experimental conditions, including cytokine treatment, cell model, dose, and duration. The output is differential log2-fold change (logFC) upon treatment.

b, Two Uses of the CytoSig framework: 1, Query a gene name to view upstream cytokine regulators or downstream target genes (if the query is a cytokine); 2, Predict cytokine signaling activities through transcriptomic profiles using a linear regression model. (Input: the input transcriptomic profile of the sample as the response variable in regression; Signature: The average response signature of cytokines as explanatory covariates; Activity: the regression coefficients reflecting signaling activities.)

c, Count of treatment response profiles with biological replicates for different molecule types.

d, Example correlation between the expression of IL10 target genes and its ligand or receptor. Each dot represents a TCGA lung adenocarcinoma sample (n=513). The X-axis shows the expression of IL10 or receptor (IL10RA + IL10RB). The Y-axis presents the expression scores of IL10 targets from a monocyte treatment experiment. Pearson correlation (r) indicates the human physiological relevance of the current data.

e, Distribution of target score correlations. Correlations were computed using all cytokine response profiles and TCGA or GTEx expression matrices. Distributions of correlations are shown by violin plots, smoothed by a kernel density estimator, for both real and randomized data through gene label permutations. The p-value, calculated with the one-sided Wilcoxon signed-rank test, represents the statistical significance of correlations being higher than zero. (n = 112 ligands and n = 111 receptors).

f, Similarity of signaling response profiles. We created a composite signature for each cytokine that consisted of the median logFC across all experiments and then calculated pairwise correlations between composite signatures for the hierarchical clustering. Red branches highlight the clusters of similar cytokines discussed in the manuscript.

The FDC automatically extracts RNASeq data from the Sequence Read Archive[12] and European Nucleotide Archive[13], along with automatically extracting MicroArray data from the Gene Expression Omnibus[14] and Array Express[15]. For metadata annotation, the FDC interacts with curators in iterative cycles. If the metadata structure and experimental designs differ drastically across studies, as was the case for cytokine-response data, the initial cycle of curation relies heavily on human expertise. However, based on the initial curations, the curators may specify automatic annotation rules, including highlighting text patterns that drive annotation decisions, translating aliases to standard names, and implementing controlled vocabularies. These natural language processing functions will dramatically reduce the human effort required after iterative cycles. The FDC is suitable for a wide range of data collection projects and is available at https://curate.ccr.cancer.gov.

Generating the CytoSig Database of Cytokine-modulated Genes

The Cytokine Signaling analyzer (CytoSig) aims to provide both a database of target genes modulated by cytokines and a predictive model of cytokine signaling activities from a sample’s transcriptomic profile (Fig. 1b). Both goals depend on an extensive data collection of cytokine-induced target genes. We first queried Array Express (AE) and Gene Expression Omnibus (GEO) databases with names and aliases of human cytokines, chemokines, and growth factors. Note that for brevity, we use the term “cytokine” at times in this manuscript to refer to these three types of signaling molecules generally. The cytokine name search yielded 9,271 candidate studies. Out of 9,271 candidates, 5,186 studies had genome-wide expression matrices and could be automatically processed by FDC.

After automatic data extraction, PhD scientists with immunology training conducted a curation of the 5,186 selected experiments (Fig. 1a). Each dataset was assigned two curators, such that the secondary curator would proofread annotations of the primary curator and correct any errors. This initial manual curation was time intensive because the metadata structure and experimental designs differed drastically across studies. However, based on the rules learned from the initial curation, the natural language processing functions from the FDC accelerated the annotation process such that minimal human effort was subsequently required. This semi-automated extraction system ensures that CytoSig will remain updated and relevant as new datasets are released.

Of the 5,186 experiments examined, 962 experiments were designated as cytokine response studies, which comprised 20,591 non-redundant individual samples (Supplementary Table 1). Curators then labeled each sample with the treatment cytokine, the cell model, treatment dose, and duration. We combined these human annotations with automatically parsed matrices of gene expression values and merged biological replicates, which generated 2,056 differential expression signatures between cytokine treatments and controls (Fig. 1c). Certain cell models and experimental conditions tend to be more frequently used than others (Extended Data Fig. 1a–f).

Extended Data Fig. 1

Curation of Human Cytokine Response Data

a, The histogram of experimental cell-type models among collected treatment profiles. For each cell-type model, we counted the number of differential expression profiles upon cytokine treatment in that model as the Profile Count.

b, The count of treatment response profiles for top-20 most frequently used cell models.

c, The histogram of treatment durations among collected treatment profiles. For each duration, we counted the number of differential expression profiles upon cytokine treatment with that duration as the Profile Count.

d, The count of treatment response profiles for top-20 most frequently used durations.

e, The histogram of experimental doses among collected treatment profiles. For each dose, we counted the number of differential expression profiles upon cytokine treatment with that dose as the Profile Count.

f, The count of treatment response profiles for top-20 most frequently used doses.

g, The association between cytokine response profiles’ quality and treatment duration. We evaluated the quality for each treatment response profile as the correlation of expression levels between target gene scores and its ligand or receptor. Each dot represents one cytokine response profile with its treatment duration on the X-axis and median correlation across all TCGA or GTEx cohorts on the Y-axis. Activin A has high quality data with long treatment durations, while EGF has high quality data with transient treatment duration.

h, The composite profile quality depends on the number of independent experiments merged. Four cytokines have more than 100 response profiles passing quality filters. We evaluated the quality of composite profiles after down-sampling the number of independent experiments merged. The quality metric is the correlation of expression between target genes and its ligand (blue) or receptor (yellow). The dots and error bars represent the median and standard deviation from 100 randomizations. Black triangle dots represent the down-sampling point that achieves 90% of the highest correlation.

i, Similarities among IFNG response signatures from diverse cell models. An average response signature was computed for each cell model. Then, average response signatures were hierarchically clustered based on Pearson correlations.

For target genes, a differential signature presents the direction of the expression change (up or down) and the magnitude of that change, expressed as log2 fold change (logFC), under each experimental condition. These differential signatures have continuous magnitude values. So, rather than using cutoffs to define cytokine targets, the differential magnitudes were used in our further analysis.

CytoSig Data Reflect Signaling Activity in Human Physiology

Since our datasets are generated through treatment experiments in cell cultures, we evaluated whether our collected cytokine targets are target genes under human physiological conditions. We measured the Pearson correlation between expression levels of the cytokine and its candidate targets in independent human tissue data. For example, we defined IL10 targets based on an IL10 treatment profile conducted in monocytes[16] and then measured the correlation between the IL10 expression and average expression scores of its candidate targets across tumors in a lung adenocarcinoma cohort[17], which we found to be 0.68 (Fig. 1d). We also found that the expression correlation between IL10 targets and IL10 receptors (IL10RA + IL10RB) is 0.62 (Fig. 1d). As a control to evaluate correlations expected by random, permuting gene identities of the IL10 treatment profile ten times resulted in a low average correlation of 0.04 between IL10 and its target genes and a low correlation of 0.05 between the IL10 receptors and its targets.

We computed correlations in this way between the expression levels of each cytokine and its candidate targets for all 2,056 cytokine treatment profiles across the Cancer Genome Atlas (TCGA)[18] and the Genotype-Tissue Expression (GTEx)[19] cohorts. The distribution of correlations between the candidate target gene expression with the respective ligands and receptors is significantly higher than expected by chance (Fig. 1e). The correspondence between the expression of target genes and a cytokine ligand or receptor in independent human tissue samples suggests that our data collection is useful for modeling cytokine signaling events in human physiology.

Although most of the cytokine response profiles derived from cell culture models are relevant in human in-vivo settings, experimental conditions of some cytokines may not reflect physiological kinetics (Extended Data Fig. 1g and Supplementary Table 2). For analyses presented from this point forward, we only use differential expression signatures with significantly positive correlations between expression levels of target genes and ligands or receptors in both TCGA and GTEx cohorts (FDR < 0.05, Online Methods). This criterion was met by 1307 out of 2056 signatures.

For each cytokine, merging independent signatures can create a composite profile with superior performance than individual signatures as measured by the correlation metrics described above (Extended Data Fig. 1h). Each cytokine’s composite signature is composed of the median logFC across all experiments for each gene, reflecting target genes induced or repressed in most conditions. We compared the overall similarities of response by performing hierarchical clustering of the composite signatures of the 43 cytokines that have at least five high-quality independent expression profiles (Fig. 1f). A few sub-clusters contain cytokines with very high correlations. For example, the composite response signature of IL27 is very similar to interferon gamma (IFNG), and to a lesser extent, to interferon type-one (IFN1) and type-three (IFNL). This observation is consistent with the downstream transcriptional similarity between IL27 and IFNG, because they both act through STAT1 signaling[20,21]. Another cluster with high similarity contains TNFA, IL1A/B, and CD40L (CD40L is both a soluble ligand and a cell surface molecule[22]), all of which activate NFKB signaling[23,24].

Although many cytokines have highly similar target responses, the same cytokine may also present context-specific differences in target response patterns. For example, the IFNG response signatures formed distinct clusters based on their cell origins. Macrophages and monocytes are clustered together and have different responses than other clusters, such as fibroblast (Extended Data Fig. 1i).

Two Regulatory Cascades from Primary to Secondary Cytokines

The hierarchical cascade of cytokine regulation is a paradigm in cellular signaling. For example, CXCL9, 10, and 11 (CXCR3 ligands) are immune-activating chemokines induced by IFNG[25], which itself is regulated upstream by IL1226. Within a signaling cascade, a cytokine can also inhibit downstream signals. For example, IL4 can block IL1 and TNFA signaling in human monocytes[27]. These hierarchical activations and inhibitions are essential to ensure rapid clearance of different pathogen classes while at the same time preventing an overzealous immune response[25-27]. The activation and repression relationships among CXCL9/10/11, IFNG, IL12, IL4, and IL1 examples discussed above are all statistically significant in our dataset (Fig. 2a and Extended Data Fig. 2a).

Fig. 2:

Cross-Regulation Hierarchies within Cytokines

a, Differential expression of target cytokines (x-axis) in response to primary cytokines (y-axis). Each dot represents a treatment profile. Blue solid points are profiles that pass quality controls, with numbers labeled on the top of box-plots. The thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. The p-value, testing whether group values are different from zero, was calculated through the two-sided Wilcoxon signed-rank test.

b, Inter-cytokine regulation hierarchy. Between each pair of regulators (X-axis) and targets (Y-axis), the log2-fold change (logFC) values are the medians among profiles that passed our quality controls. The plot only includes regulators and targets with at least one logFC larger than two. The upper heatmap includes target ligands, and the lower part includes target receptors with ligands in brackets. The NFKB and interferon transcriptional groups are labeled with boxes.

c, Differential logFC values of GZMA, GZMB, and PRF1 upon TGFB1 treatment, shown by violin plots, smoothed by a kernel density estimator. Each dot represents an independent profile. The p-value, testing whether the T-cell group values are different from zero, was calculated through the two-sided Wilcoxon signed-rank test (n = 11 independent treatment profiles for each gene).

d, Perforin intracellular levels measured by flow cytometry upon TGFB1 treatments in human primary CD8 T cells. The X-axis shows the signal intensity. The Y-axis shows the T-cell fraction with modal normalization, scaling the maximum Y-axis to 100%. Vertical line indicates percentage of T cells with signal intensity above the gate threshold (defined in Extended Data Fig. 2d).

e, Granzyme and Perforin protein levels upon TGF-beta treatments in human and mouse CD8 T cells. The mean percentage of T cells with an intensity above the gate threshold is shown, with standard deviations as error bars (n = 3 cell-culture replicates per condition). The two-sided Wilcoxon rank-sum test was used to compare groups.

Extended Data Fig. 2

Target Genes in Response to Cytokine Treatments

a, Differential expression of target cytokines (x-axis) regulated by IL4 or BMP6 (y-axis). The thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range. The difference between group values and zero among treatment profiles that passed quality controls was tested using the two-sided Wilcoxon signed-rank test, with p-values and sample counts labeled.

b, Gene set enrichment among cytokine response targets. The normalized enrichment score from the gene set enrichment analysis (GSEA)[6] represents the overall enrichment of pathway members among target genes from each cytokine’s composite signature.

c, Normalized enrichment scores from panel b among pairs of hallmark categories and cytokines.

d, Gating strategy of flow analysis. The forward scatter height (FSC-H) and side scatter height (SSC-H) are used for dead cell and cell debris removal. FSC-H width and SSC-H width are used to gate the single cells. We include unstained cells to define the gate threshold that separates positive populations and negative control cells.

We systematically examined the changes induced by primary cytokines on secondary signals and identified two distinct pro-inflammatory clusters (Fig. 2b). Each group activates a distinct set of secondary targets. The first group, including TNFA, IL1, IL17A, IL36, and CD40L, triggers IL1, 6, CXCL1, 2, 5, 6, 8, and CCL20 (Fig. 2a and Fig. 2b, left box). These primary cytokines have target genes enriched in the NFKB signaling pathway (Extended Data Fig. 2b). The target chemokines of this group may attract or activate pro-inflammatory immune cells, such as neutrophils, fibroblasts, and T cells[28]. We hereafter refer to this first group the “NFKB transcriptional group”.

The second group, including IFNG, IFN1 (IFNA and IFNB), IFNL, and IL27, trigger CXCL9, 10, 11, and TRAIL (Fig. 2b, right box). These primary cytokines are known to be interferon-related because IL27 has similar downstream transcriptional profiles with IFNG[20] through STAT1 signaling[21]. Their secondary targets are chemokines for activated T cells (CXCL9, 10, 11)[25], or pro-apoptotic signals released by effector T cells (TRAIL)[28]. We hereafter refer to this second group as the “interferon transcriptional group”.

Besides regulating a ligand, cytokines can also modulate receptor activity as an alternative means of cascade regulation. For example, CytoSig finds that Activin A activates CXCR4 while GMCSF represses CXCR4 (Fig. 2a, b). This result is consistent with previous studies of the CXCR4 regulation[29,30]. Regulation of receptors in this way appears less frequent than ligand regulation: we observed just 4 out of 183 (2.2%) annotated receptors that had a logFC greater than two, whereas 33 out of 253 (13%) annotated ligands have a logFC larger than two (Fig. 2b).

The CytoSig Data Reveal Anti-Inflammatory Cytokines

In contrast to cytokines that induce secondary cytokines and chemokines, IL4 and BMP6 repress other pro-inflammatory molecules, such as IL1B, CXCL1, 8, and CCL2 (Fig. 2b and Extended Data Fig. 2a). Gene Set Enrichment Analysis (GSEA) on target genes of IL4 and BMP6 revealed a depletion of inflammatory response pathways (Extended Data Fig. 2b, c). IL4 is a well-known anti-inflammatory cytokine that inhibits certain immune processes, although it may also cause allergic inflammation in a context-dependent manner[31].

Besides IL4 and BMP6, which directly suppress the transcription of downstream cytokines and chemokines, other anti-inflammatory molecules may counteract inflammation by alternative mechanisms. For example, a previous study in mouse models demonstrated that TGF-beta signaling directly targets cytotoxic T-cell functions in mice[32]. Indeed, our collected data shows that TGFB1 treatment in human T cells significantly down-regulates granzyme A (GZMA), B (GZMB), and Perforin (PRF1), which induce cell death in target cells attacked by T cells (Fig. 2c). Flow cytometry analysis in human and mouse primary T cells validated the inhibitory effect of TGFB1 on GZMA, GZMB, and PRF1 (Fig. 2d, e and Extended Data Fig. 2d (gating strategy)). Therefore, our data can reveal broad categories of anti-inflammatory cytokines.

We next examined how cytokines cooperate with and antagonize each other with respect to target genes across the human genome. To test whether pairs of cytokines co-regulate target genes, we enumerated genes with significant logFC values from both cytokines under analysis. Then, we compared the gene counts against values when gene labels are shuffled to compute false discovery rates (FDR) (Online Methods). We defined significant results with an FDR threshold of 0.05 (Fig. 3a). In 86% of statistically significant cases, cytokine pairs either enhance or repress target genes in concert (Fig. 3b). For example, TNFA and IL1B induce a similar set of genes and also repress a similar set of genes when investigating average targets across all models in our data collection (Fig. 3a, left panel).

Fig. 3:

Antagonizing Interactions between Cytokines on Target Genes

a, Co-regulation between IL1B and other cytokines on target genes. Each dot represents a target gene. The X and Y axes show the median log2-fold change (logFC). Blue (co-enhance) represents targets with significant positive values on both axes (false discovery rate < 0.05). Green (co-repress) represents targets with significant negative values on both axes. Orange (antagonize) represents targets with a significant positive score from one signal and a negative one from the other. Selected antagonized targets are labeled.

b, Fractions of cytokine co-regulation type. For each cytokine pair and their targets, we counted the fractions of significant co-regulation in the three categories illustrated in panel a.

c, Antagonistic relationships among cytokines on target genes. Each flat-end edge indicates that the first cytokine represses targets of the second one. The edge width is proportional to the number of antagonized targets. The node pie chart represents the in and out degrees.

d, Antagonizing regulatory relationships between cytokines and anti-inflammatory signals. The heatmap presents the number of target genes induced by a cytokine (row) but repressed by an anti-inflammatory signal (column). The similarities are shown with the hierarchical clustering trees with Pearson correlation as the distance metric.

In 14% of the statistically significant cases, cytokines exhibit an antagonistic relationship, meaning that they have opposite signaling effects on downstream targets (Fig. 3b). For example, IL4 and BMP6 downregulate many targets induced by IL1B and TNFA (Fig. 3a and Extended Data Fig. 3a). We also observed a similar relationship among four other cytokines: the IFNG target genes are antagonized by IL10 and GCSF but enhanced by IL27 (Extended Data Fig. 3b). Thus, our target co-regulation analysis identified four major anti-inflammatory regulators (IL4, BMP6, IL10, and GCSF), which antagonize chiefly pro-inflammatory molecules in two groups (Fig. 3c, d), which we referred to as the NFKB and interferon transcriptional group (Fig. 2b).

Extended Data Fig. 3

Co-regulation between Cytokines on Target Genes

a, Target co-regulation between TNFA and other anti-inflammatory cytokines.

b, Target co-regulation between IFNG and other cytokines.

BMP6 Antagonizes Pro-Inflammatory Cytokine Target Genes

Previous work demonstrated that IL4 could inhibit NFKB transcription programs[33], explaining the antagonistic relationship between IL4 and cytokines in the NFKB group. A previous study demonstrated that BMP6 could inhibit the CCL2 mRNA level induced by TNF[34]. However, to the best of our knowledge, no previous studies have reported BMP6 as an anti-inflammatory molecule that antagonizes many proinflammatory targets.

Our analysis indicated that BMP6 may antagonize the effect of IL1B through down-regulating IL1B-induced pro-inflammatory chemokines, with CXCL8 and CCL2 as the most significant targets (Fig. 3a). To validate our prediction, we first evaluated the intracellular protein levels of CXCL8 and CCL2 upon IL1B and BMP6 treatments by flow cytometry in two human lung epithelial cell lines A549 and NCI-H1299. Consistent with our data analysis, BMP6 treatment significantly inhibits the IL1B induction of CXCL8 and CCL2 (Extended Data Fig. 4a–c). ELISA assays also indicated that levels of soluble CXCL8 and CCL2 are consistently lower in cells treated with BMP6+IL1B compared to IL1B alone (Extended Data Fig. 4d, one-sided Wilcoxon signed-rank p-value = 0.016).

Extended Data Fig. 4

BMP6 represses CXCL8 and CCL2 induced by IL1B

a, Representative plots of CXCL8 and CCL2 protein levels upon BMP6 and IL1B cotreatments in A549 cells. We utilized flow cytometry to measure CXCL8 and CCL2 intracellular protein levels after 12-hour treatments with IL1B, IL1B+BMP6, and media control. The X-axis shows the signal intensity measured by flow cytometry. The Y-axis shows the A549 fraction distribution with modal normalization, scaling the maximum Y-axis value to 100%. The percentage of cells with signal intensity above the gate threshold (vertical line, panel c) is indicated.

b, Summary plots of CXCL8 and CCL2 protein levels upon BMP6 and IL1B treatments. In three cell-culture replicates, A549 and H1299 cells were treated either simultaneously or sequentially with combinations of BMP6 and IL1B. The mean fraction of cells with an intensity above the gate thresholds (defined in panel c) is plotted with standard deviations as error bars (n = 3 cell-culture replicates per condition). The two-sided Wilcoxon rank-sum test p-values were computed to compare groups.

c, Gating strategy of flow analysis. The forward scatter height (FSC-H) and side scatter height (SSC-H) are used for dead cell and cell debris removal. FSC-H width and SSC-H width are used to gate the single cells. We include unstained cells to define the gate threshold that separates positive populations and negative control cells.

d, CXCL8 and CCL2 soluble protein levels by ELISA in A549 cells. A549 cells were treated with BMP6 and IL1B in combinations for 24 hours. The mean soluble cytokine levels were measured by ELISA with standard error of the mean as error bars (n = 3 independent experiments).

CytoSig Predicts Signaling Activities from Expression Data

Since the cytokine response data in our collection reflect signaling relationships in human tissues (Fig. 1e), we created the CytoSig model to predict signaling activities using the transcriptome profile of an input sample. The output of CytoSig is different from standard cytokine assays such as ELISA, which measure cytokine levels instead of cytokine target activities. CytoSig utilized the ridge regression to search for features in an input transcriptome profile that can be explained by a cytokine’s influence on its target gene expression (Extended Data Fig. 5 and Online Methods).

Extended Data Fig. 5

CytoSig model selection

a, Prediction performance measured as the cross-validation (CV) R2 in TCGA cohorts. For each algorithm, median 5-fold CV R2 metrics across all TCGA datasets were shown at different penalty values with standard error of the mean (SEM) as error bands. The horizontal line indicates 70% of the optimal CV R2 of ridge regression. The vertical line marks the lambda value of 10,000, the penalty used in the CytoSig model.

b, Prediction performance measured as the cross-validation (CV) R2 metrics in GTEx cohorts, shown as panel a.

c, Inference performance measured as the correlation between model coefficients and cytokine expression in TCGA cohorts. For each algorithm, we computed the median correlation values between model coefficients and cytokine ligand or receptor expression at different penalty values, with SEM as error bands. The vertical line marks the Lambda value of 10000, which is the penalty reaching 80% of the optimal correlation. XG Boost with tree learner cannot be evaluated for inference performance because its tree structure cannot provide coefficients as cytokine response.

d, Inference performance measured as the correlation between model coefficients and cytokine expression in GTEx cohorts, shown as panel c.

As described in the introduction, redundancy and pleiotropy are major obstacles to modeling cytokine activity. To account for complications from signaling pleiotropy, our model only aims to predict each cytokine’s overall activity, instead of its effects on individual genes or pathways. We analyzed each cytokine’s composite signature averaged across at least five independent experiments. Significant enrichment of the composite signature of a cytokine in the input sample’s transcriptome should indicate the presence of signaling events. To address signaling redundancy, we utilized a penalized linear model that avoids reporting a cytokine as active if other cytokines with similar composite signatures have influenced target gene expression to a greater extent. For any input profile, our model reports a significant score for a cytokine only if predicted activities are significantly higher than expected by chance (Online Methods).

Accuracy Validation by Cytokine-Blocking Clinical Response

To test the model accuracy, we reasoned that the patient’s clinical response upon cytokine blocking therapies should reflect authentic cytokine activities in human tissues. Therefore, we compared CytoSig predictions of cytokine activities with transcriptomic data before and after cytokine-blocking therapies in inflammatory diseases (Supplementary Table 3). For example, a microarray study measured the whole blood transcriptome of arthritis patients at baseline and day 3 after anti-IL1B Canakinumab treatmentand evaluated the therapy response at day 15 post-therapy[35]. Upon IL1B neutralizing therapy, the IL1B activity reduction at day 3 predicted by CytoSig correlates significantly with the patient clinical response evaluated at day 15 (Fig. 4a). For another example, an IFNA vaccine trial among systemic lupus patients profiled both whole blood transcriptomes and clinical response as the titer of IFNA neutralizing antibodies in blood after immunization[36]. The IFN1 activity reduction predicted by CytoSig correlates significantly with the clinical response across patients (Fig. 4b).

Fig. 4:

CytoSig Predicts Cytokine Activities in Human Diseases

a, IL1B activities in blood predict anti-IL1B therapy response in arthritis. Each dot represents a blood sample[35]. The X-axis shows patients with similar therapy responses. The Y-axis shows the IL1B activities predicted by CytoSig at an early time point, shown by violin plots smoothed by a kernel density estimator. Spearman correlation between clinical response in patient groups and the median IL1B activity with a two-sided t-test p-value is indicated.

b, IFN1 activity in blood correlates with antibody titer upon IFNA vaccine in systemic lupus. Each dot represents a blood sample. The X-axis presents the titer of IFNA antibody post immunization[36]. The Y-axis presents the IFN1 differential activity predicted by CytoSig. Spearman correlation between X and Y axes with a two-sided t-test p-value is indicated (n=36).

c, Predicted activity change upon anti-cytokine therapies. Each dot represents an anti-cytokine therapy study in inflammatory diseases, with targets on the X-axis. Grey labels represent mouse model studies for cytokines without clinical data. The Y-axis presents the average differential activity between post- and pre-treatments. The accuracy is computed as the fraction of cytokines with median activity reduction smaller than one, with p-value from the two-sided Wilcoxon signed-rank test.

d, TGFB activity changes upon neutralizing antibody treatments. The first antibody inhibits all TGFB isoforms (123, 7 mice) and the second antibody inhibits only TGFB1 and 2 (12, 6 mice). The anti-TGFB3 profile (123/12) was the differential profile between the pan-TGFB and the TGFB1,2 specific groups. The TGFB1 and TGFB3 activities predicted by CytoSig were shown by bar plots, with two-sided p-values from the permutation test with 10,000 randomizations.

e, VEGF activities in pre-treatment tumors predict anti-VEGF therapy response, in Sunitinib[38] and Bevacizumab clinical studies[39]. The Kaplan-Meier plot presents patient fractions (Y-axis) with survival length (X-axis, progression-free or overall) among pre-treatment tumors with high and low VEGFA activities predicted by CytoSig. The activity cutoff is selected through maximizing the difference between high and low groups. The p-value was from the one-sided Wald test using continuous values without cutoffs.

Among all cytokine-blocking studies collected from GEO and ArrayExpress databases, CytoSig predicted the activity reduction score to be at least negative one (one standard deviation below zero) for 85% cytokines (Fig. 4c). The accuracy dropped to 0% when gene labels are permuted in the model. These results supported the reliability of CytoSig on cytokine activity prediction in human tissues and demonstrated the clinical utility to guide therapy decisions.

Accuracy Validation on TGF-beta Isoform-Specific Activities

CytoSig predicts different activities for cytokines from the same family sharing receptors, such as TGFB1 and TGFB3. The validation in the previous section established that CytoSig can perform with high accuracy on a broad set of cytokines. To validate the accuracy of CytoSig’s predictions of signaling activities among cytokine isoforms sharing the same receptors and similar downstream pathways, we performed in-vivo experiments with the 4T1 breast cancer mouse model using neutralizing antibodies against TGF-beta isoforms.

Specifically, we profiled the transcriptomes of mouse 4T1 tumors treated with neutralizing antibodies targeting all TGF beta isoforms and antibodies targeting only TGFB1, 2 (but not TGFB3). The differential profile between pan-TGFB and TGFB1,2 antibodies can reflect the anti-TGFB3 effects because the TGFB3 isoform is the differential target between two antibodies. CytoSig predicts a significant reduction of TGFB1 activity based on the differential transcriptomic profiles upon treatments for both anti-TGFB antibodies, and a significant reduction of TGFB3 activity only for the anti-TGFB3 profile (Fig. 4d).

Accuracy Validation in Tumors and Cancer Therapy Response

To further evaluate CytoSig model accuracy, we utilized the International Cancer Genome Consortium (ICGC) tumor cohort[37], which has no overlap with the previous TCGA and GTEx data in model training. We assumed that tumors with ligand or receptor expression levels higher than one standard deviation above the average level in the entire dataset reflect positive activity for that cytokine. Under this assumption, we evaluated the accuracy of CytoSig on predicting samples with significant cytokine signaling activities. Based on the Receiver Operating Characteristic (ROC) curve and area under the ROC curve (AUC), 35 out of 43 cytokines have a performance significantly better than chance (Extended Data Fig. 6a, b). Therefore, the CytoSig model can predict target activities of most cytokines.

Extended Data Fig. 6

CytoSig Predicts Cytokine Activities in Tumors and Cancer Therapy Response

a, Receiver Operating Characteristic (ROC) curve of IFNG activity prediction. For each dataset from the ICGC, the ROC curve presents false-positive rates against true-positive rates at different IFNG activity thresholds.

b, Area under the ROC Curve (AUC) as the prediction accuracy. We computed the ROC curves for all cytokines following the procedure in panel e. Each bar represents the median AUC among all ICGC cohorts, with standard errors of the mean as error bars (n=9 independent datasets). The AUC baseline is 0.5, representing a random prediction. We applied the one-sided Wilcoxon signed-rank test to evaluate whether AUC values are higher than 0.5 for each signal, and converted p-values to false discovery rates (FDR) through the Benjamini-Hochberg correction.

c, IFNG activity predicts overall survival upon Atezolizumab treatment in urothelial carcinoma[43]. The Kaplan-Meier plot presents patient fractions (Y-axis) with different overall survival (X-axis) among pre-treatment tumors with high and low IFNG activities predicted by CytoSig. The activity cutoff is selected through maximizing the difference between high and low groups. The p-value was from the two-sided Wald test using continuous values without cutoffs.

d, Signaling activity computed by CytoSig better predicts clinical outcome than other metrics. We compared three approaches to compute cytokine activities, including expression levels of ligand, receptors, and CytoSig predictions. For IFNG activity, we also utilized a geneset signature developed by Merck to predict checkpoint blockade response[41], as well as the PDL1 expression. The association between activity and survival outcome was computed as the Wald test z-score in Cox-PH regression.

We also evaluated the CytoSig model in predicting the clinical outcome of anticancer therapies that inhibit cytokine signaling. Vascular endothelial growth factor (VEGF)blocking is a category of treatments inhibiting either VEGF ligands or receptors from promoting abnormal angiogenesis in tumors[38,39]. As the cancer driver, the pre-treatment target pathway activity may predict targeted therapy efficacy and patient survival after treatment[40]. We found that high VEGF signaling activities predicted by CytoSig in pre-treatment tumors, using data from two clinical studies[38,39], are highly predictive of longer survival outcomes upon blocking the VEGF pathway through either ligand (Bevacizumab) or receptors (Sunitinib, inhibitor of multiple receptor tyrosine kinases, including VEGF receptors) (Fig. 4e).

Immune checkpoint blockade is another treatment category whose responses depend on cytokine signaling by IFNG[41]. PDL1 is a target gene induced by IFNG signaling[42]; thus, we evaluated the association between IFNG activity in pre-treatment tumors and the anti-PDL1 therapy response, using data from an anti-PDL1 clinical trial in urothelial cancer[43]. IFNG activity predicted by CytoSig is highly predictive of overall survival outcome upon anti-PDL1 (Extended Data Fig. 6c). Moreover, for both anti-VEGF and anti-PDL1 clinical studies, the CytoSig predictions had better associations with the clinical outcome than other approaches, such as ligand or receptor expression and gene set signatures (Extended Data Fig. 6d).

Accuracy Validation in Single-cell Transcriptomic Data

Encouraged by the reliable performance on bulk data, we further evaluated the capability of the CytoSig to predict signaling activities in single cells. The ideal evaluation standard for CytoSig predictions in single cells would be a method providing systematic measurements of both transcriptome and cytokine activities in each single cell. However, no such method currently exists to our knowledge. To validate the accuracy in single-cell data, we utilized transcription factor (TF) activities as indicators of active cytokine signaling (Supplementary Table 4).

We computed TF activities for a single-cell transcriptomic profile using the RABIT framework, which leverages an extensive collection of ChIP-Seq profiles to predict TF activities through transcriptional patterns of TF target genes[44]. For example, using data from a COVID-19 single-cell study[45], RABIT predicted that most CD8 T cells have a positive STAT1 TF activity, reflected as a higher expression level of STAT1 ChIP-Seq target genes compared to other genes (Fig. 5a). A minor CD8 T cell population has negative STAT1 TF activity. Consistent with the dependence of interferons and IL27 on STAT1 signaling (Supplementary Table 4), cells with positive TF activities have significantly higher signaling activities from the CytoSig model than cells with negative TF activities (Fig. 5a).

Fig. 5:

CytoSig Predicts Cytokine Activity in Single-Cell RNASeq Data

a, Single-cell signaling activities for STAT1-group cytokines among CD8 T cells using data from a COVID-19 study[45]. For each single-cell (dot), the X-axis shows the STAT1 transcription factor (TF) activities. The Y-axis presents the average signaling activities among interferons and IL27. The color brightness indicates the local dot density. The p-value was calculated through the two-sided Wilcoxon rank-sum test, comparing cytokine activities of single cells with positive (n=1323) and negative (n=117) STAT1 activities.

b, Receiver operating characteristic (ROC) curves for CD8 T cell effector cytokines. The ROC curve presents false positive rates against true positive rates of TF activity prediction at different cytokine activity thresholds using the data in panel a. The diagonal line represents random expectation.

c, Accuracy of single-cell cytokine activities prediction. The area under the ROC curve (AUC) for each cell type was computed for each TF and its cytokine families, using the COVID-19 single-cell dataset analyzed in panel a and b. Each dot represents a cell type (n=12 cell types per boxplot). The thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range.

d, Accuracy of single-cell cytokine activities prediction among many single-cell datasets. Each dot presents the median AUC among all cell types (Y-axis) in a single cell dataset. The results were shown for all cytokine and TF combinations using box plots as panel c (n=10 cytokine-TF combinations per boxplot).

We utilized a receiver operating characteristic (ROC) curve to measure the CytoSig’s ability to predict TF activities based on the predicted cytokine activity. We found that the activity of CD8 T cell effector cytokines, including interferons and TNF, all predict downstream TF activities better than random (Fig. 5b). Using the area under the ROC curve (AUC), we next evaluated cytokine activities’ predictive performance on their downstream TF activity for all cell types in the COVID-19 single-cell study. The AUC metrics are consistently higher than expected by chance for ten out of eleven pairs of cytokines and downstream TFs (Fig. 5c and Supplementary Table 4). We observed similar high performance in another cancer study (Extended Data Fig. 7a). We performed such evaluation on 18 single-cell datasets and found that AUC metrics are consistently higher than random (Fig. 5d).

Extended Data Fig. 7

CytoSig Reliably Predicts Cytokine Activity in Single Cells

a, Result for a liver cancer cohort[72] (n=11 cell types per boxplot). The thick line represents the median value. The bottom and top of the boxes are the 25th and 75th percentiles (interquartile range). The whiskers encompass 1.5 times the interquartile range.

b, Result for a COVID-19 peripheral blood cohort[46] (n=15 cell types per boxplot). The boxplot is defined as panel a.

CytoSig Identifies Signaling Markers of Severe COVID-19

The global spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an urgent health crisis. The symptoms of coronavirus disease 2019 (COVID-19) range from mild fever, cough, and difficulty breathing to respiratory failure and death[7]. While more severe outcomes have been associated with an exaggerated immune response referred to as the cytokine storm, the immune-response mechanisms underlying the dramatic difference of disease severity remain unclear.

We applied CytoSig to analyze single-cell RNASeq data from bronchoalveolar lavage fluid[45] and peripheral blood[46] of COVID-19 patients (Fig. 6a). These datasets are used earlier in this manuscript to establish reliable prediction accuracy in single-cell RNASeq data (Fig. 5c, d and Extended Data Fig. 7b). Many cytokine signaling activities are significantly associated with the severity differences in COVID19 symptoms. For example, among macrophages from lavage, IL10 activity is significantly higher in severe patients than mild and healthy controls (Fig. 6b, c, and Supplementary Table 5). In contrast, among CD8 T-cells from blood, the type I interferon (IFN1) activity is highest in mild patients relative to severe cases or healthy controls (Fig. 6d, e). This is consistent with previous studies reporting a lack of type I interferon response among both patients[47] and cultured cells upon SARS-CoV-2 infection[48].

Fig. 6:

Signaling Features underlying the COVID-19 Symptom Difference

a, Sample source sites. We analyzed single-cell RNASeq datasets from the broncho-alveolar lavage and peripheral blood samples from COVID-19 patients.

b, IL10 activities in Macrophage from lavage. Each dot presents a single cell from the COVID-19 study with t-SNE (t-distributed stochastic neighbor embedding) coordinates computed as the original publication[45]. The color represents the predicted IL10 activity. The circles highlight the cluster of macrophages.

c, Enrichment of IL10 activities among severe patients in lavage macrophage. The violin plots present IL10 activity distributions in different patient groups, smoothed by a kernel density estimator. The color legend is available at the bottom of panel b. The two-sided Wilcoxon rank-sum p-value is computed to compare severe (n=6 patients) and mild (n=3 patients) group activities.

d, IFN1 activities in CD8 T cells from blood. The IFN1 activities in CD8+ T cells (circles) are shown as panel b, with UMAP (Uniform Manifold Approximation and Projection) coordinates from the original publication[46].

e, Depletion of IFN1 activities among severe patients in peripheral blood CD8 T cells. The two-sided Wilcoxon rank-sum p-value is computed to compare severe (n=14 patients) and mild (n=13 patients) group activities within each cell type.

f, Differential signaling activities in macrophage from lavage. The heatmap presents cytokines whose predicted signaling activities are significantly different between severe, mild, healthy individuals.

g, Summary of signaling activities in COVID-19. The heatmap includes cell types with significant differences in signaling activities among severe, mild, and healthy individuals. Each cell shows the median value across all individuals in each group. Only cytokines with at least three significant values in at least three cell types are included.

Analysis of differential activity among severe, mild, and healthy individuals revealed a few cytokines as signaling markers of COVID-19 symptom severity. Among severe COVID-19 patients, we found elevated IL1B activities in macrophages, the most populous cell type from bronchoalveolar lavage fluid[45] (Fig. 6f). This result of IL1B activation in macrophages is consistent with previous reports[7]. Meanwhile, IL10 has high activities among severe patients in lavage and blood myeloid cells, including macrophage, monocytes, and neutrophils (Fig. 6g and Extended Data Fig. 8), a result also consistent with published studies[49,50]. IL10 is known to inhibit antigen presentation by dendritic cells to cytotoxic T cells, thus impairing T cell-mediated antiviral immunity[51]. IL10 can also suppress macrophage activation for fighting against intracellular pathogens[52]. Thus, in severe COVID-19 cases, the cytokine environment may compromise the antiviral immune response while triggering pathological inflammation. We observed that monocytes and macrophages have high expression levels of IL10 and IL1B, thus potentially serve as cytokine producing cells (Extended Data Fig. 9).

Extended Data Fig. 8

Differential signaling activities from COVID-19 samples

The heatmap presents cytokines whose predicted signaling activities are significantly different among severe, mild, healthy individuals in a COVID-19 peripheral blood cohort[46] (false discovery rate < 0.05).

Extended Data Fig. 9

Gene expression of cytokines with differential signaling activities from COVID-19 samples

IL1B and IL10 gene expression in diverse cell types. The violin plots present the expression of IL1B and IL10 in patient groups, with distributions smoothed by a kernel density estimator.

CXCL8 as Candidate Therapeutic Target for Severe COVID-19

Direct IL1B or IL10 blockade may compromise the antiviral response[53] or fail to alleviate inflammation[54]. Therefore, inhibition of downstream targets could serve as alternative approaches. With this in mind, we analyzed the downstream targets induced by IL1B or IL10. We found that IL1B target CXCL8 have higher expression levels in macrophages from lavage samples of patients with severe COVID-19 (Extended Data Fig. 10). Meanwhile, in blood neutrophils, IL10 target CXCL1 has higher expression levels among severe patients than mild cases and healthy controls (Extended Data Fig. 10).

Extended Data Fig. 10

Potential Therapeutic Targets to Overcome COVID-19 Induced Tissue Damage

a, Target genes of cytokines with elevated activities among severe patients. Arrow-headed edges indicate up-regulation, and flat-headed edges indicate down-regulation. For each cell type, red square nodes represent cytokines with high activity scores among severe patients, and blue diamonds represent anti-inflammatory signals with low activities. The network only includes targets whose expression values are significantly higher compared between severe and mild patients, and between disease and healthy controls.

b, Example of gene expression in different patient groups. The violin plots present gene expression distributions in patient groups, smoothed by a kernel density estimator. Examples from lavage macrophages[45] are in the left, and peripheral blood neutrophils[46] are in the right. Y axis indicated values from individual cells.

CXCL1 and CXCL8 all bind to the CXCR2 receptor and serve as primary chemokines in neutrophil recruitment. A high neutrophil-to-lymphocyte ratio in peripheral blood indicates severe disease and organ failure[55]. Aberrant formation of neutrophil extracellular traps may contribute to severe damage to the lung parenchyma in COVID-19[56]. A phase I clinical trial evaluated the CXCL8 blocking antibody in treating solid tumors and has not observed any dose-limiting toxicities[57], which indicates the potential of therapy repurposing.

Discussion

We have introduced CytoSig, a data-driven platform to model cytokine activities. CytoSig complements existing cytokine release assays because it can predict cytokine target activities from bulk transcriptomic data available from many large-scale cohorts and single-cell RNA-seq data that provides resolution down to individual cells. The acquisition of both types of data is now routine, making CytoSig useful to a broad spectrum of research questions.

CytoSig offers particular advantages in analyzing single-cell data because it is not affected by the absence of cytokine-producing cells or zero read counts for ligand or receptor genes. This advantage is especially important because current single-cell technologies have difficulty capturing some cell types, such as neutrophils[58]. Many studies also sort cells using markers such as CD45, which may exclude cytokine-producing cell populations. Moreover, the dropout events, reflected as zero read counts, on transiently expressed cytokine genes further complicate analysis. The CytoSig model uses a reliable alternative strategy, analyzing receiver cells’ transcriptional patterns across many cytokine target genes.

A limitation of CytoSig is the ascertainment bias of public datasets, which leads to many experiments on a few cytokines, cell models, or experimental conditions, and a lack of data on others. There are currently 67 human cytokines, 42 chemokines, and 133 growth factors annotated in the literature[28]. However, our collection from public databases captures high-quality profiles for 43 of the 242 documented molecules (17.8%) due to the lack of data available for most signaling molecules. Also, most datasets were generated through a few models, such as monocytes or fibroblasts, without sufficient coverage on diverse cell lineages. Such a gap indicates a need for attention on a broad range of cytokines and cell models beyond a few deeply studied molecules and systems.

Despite these limitations, the CytoSig platform provides biologists and clinicians with a powerful resource to study signaling activities in laboratory or clinical samples. Furthermore, independent of CytoSig, the Framework for Data Curation (FDC) is a general resource for data scientists to accelerate data curation projects. Using the FDC, our plans to continuously integrate new datasets will provide the community with an ever-growing repository for generating new biological insights.

Online Methods

The Framework for Data Curation (FDC) from Public Databases

The FDC aims to automate the data curation process as much as possible with two components: 1, semi-automated metadata annotation; 2, automatic gene expression matrix extraction. Databases processed by FDC include ArrayExpress (AE)[15], Gene Expression Omnibus (GEO)[14], European Nucleotide Archive (ENA)[13], and Sequence Read Archive (SRA)[12]. The FDC server is built using Python 3 and Dango 3 frameworks with MySQL 8 as the database backend. The natural language processing functions are created in the web browser frontend using the jQuery 3 Javascript library.

The first FDC component for metadata annotation utilizes a three-stage approach. In the first stage, the users should query the GEO and AE databases with keywords related to their biological topic. The SRA and ENA metadata are available through the GEO and AE, respectively. The database query will generate a list of candidate datasets. After uploading the candidate list to FDC, the users can define pattern-matching rules, implemented as regular expressions[59], to narrow down query results. In the second stage, users should browse study summaries and determine which datasets are relevant to their study topic. To accelerate the process, users could define a set of highlighting rules, implemented as regular expressions[59], so that curators only need to focus on the most relevant texts.

In the third stage, users will extract metadata fields for each experimental profile. The FDC aims to reduce human manual edits as much as possible, with automatic rules and text transformation functions defined by users. FDC will automatically parse each sample’s study design information and summarize all potentially relevant fields in a candidate table. The users can define a set of automatic mapping rules to convert aliases, such as biological molecule and cell model names, to their standard names. FDC also provides automated functions to extract and transform text information from candidate metadata columns. Based on these functions, curators will standardize metadata columns into controlled vocabularies.

The second FDC component can automatically extract MicroArray and RNA-Seq public databases. For Affymetrix MicroArray data on the AE and GEO, we generated expression matrices from CEL files through the R Oligo package[60]. For other MicroArray platforms, we downloaded the processed data through the R GEOquery[61] and Python Orange[62] packages for GEO and ArrayExpress, respectively. For RNASeq data from the ENA and SRA databases, we downloaded the fragments per kilobase of transcript per million mapped (FPKM) data through the RNASeq-ER Application Programming Interface[63]. In total, our framework extracted 27,181 independent human transcriptomics datasets deposited before February 2nd, 2020. Many datasets from public repositories are not gene expression studies, thus cannot be automatically analyzed by the second FDC component. However, the first FDC component can still assist the metadata annotation for non-transcriptomic data.

Besides the two primary components introduced above, FDC also provides other assistant modules, such as curator management for a project and result proofreading panels if project managers want to review the annotations from curators. With the standardized metadata matrix from the first component and gene expression matrix from the second component, users can perform algorithmic data analysis for their biological topics.

Collection of Human Cytokine Response Data Based on FDC

In the first step, we queried names and aliases of human cytokine, chemokine, growth factors[28], and a few immune-suppressive signals in the tumor microenvironment[64] through the query interface of GEO and ArrayExpress. The SRA and ENA metadata are available through the GEO and ArrayExpress, respectively. Our query returned 9271 candidate series, and 5186 of them have processed data matrices from FDC. The other 4085 datasets do not have FDC-processed data due to several reasons. Some studies using NanoString platforms only focus on hundreds of genes instead of genome-wide. Some MicroArray or RNA-Seq studies may have corrupted raw data, leading to FDC extraction failures. We also excluded all micro-RNA and non-coding RNA studies.

In the second step, we recruited Ph.D. scientists with immunology training for data annotations based on the FDC. Curators have focused on data curation for several months. A second curator proofreads all annotations of the first curator and corrects errors. Most studies among the 5186 candidates only mentioned cytokines in their description but did not study the signaling response. The curators read descriptions of 5186 experiments and identified 962 of them as cytokine response studies, including 20591 non-redundant samples (Supplementary Table 1). Then, curators read descriptions of 20591 samples and label cytokine treatment, cell model, dose, and duration, using the semi-automated functions on the FDC. We established a set of control vocabularies about signal names, cell models, concentrations, and duration units.

Together with data matrices extracted and expert annotations, we generated differential gene expression profiles, defined as the log2 fold change (logFC) between treatment and control conditions. We only kept experiments with biological replicates and acquired 2056 logFC vectors after merging biological replicates. Meanwhile, we merged IFNA and IFNB as IFN1, representing type 1 interferons, due to the high Pearson correlation of 0.698 between their composite profiles. We also combined IL36A and IL36B as IL36 due to the high correlation of 0.938 between their composite profiles.

Data Quality Control

To test the human physiological relevance of data collection, we defined a quality control metric as the Pearson correlation between expression levels of cytokine target genes and the ligands or receptors in independent human tissue data. We utilized the Cancer Genome Atlas (TCGA)[18] and Genotype-Tissue Expression (GTEx)[19] datasets. Each TCGA and GTEx sample measures a bulk tissue’s average expression that contains both producer and receiver cells for cytokines.

To measure the overall expression of target genes, we performed a linear regression for each pair of cytokine response and tissue expression profile as “tissue expression = A * cytokine profile + B,” and computed the cytokine profile’s target score as A / standard error(A) using the ordinary least squares[65]. The target score represents the enrichment of a response signature in the tissue expression profile. We then analyzed the Pearson correlation between the target score and ligand or receptor expression across tissue samples (Fig. 1e).

TCGA has 33 tumor cohorts, and GTEx has 27 tissue cohorts. For each cytokine profile, we utilized the one-sided Wilcoxon test to evaluate whether the correlations with the ligand or receptor are higher than zero across both TCGA and GTEx cohorts (false discovery rate < 0.05 with the Benjamini-Hochberg correction). We only included 1307 profiles that passed the threshold in further analysis.

Target Cooperation and Antagonization Analysis between Cytokine Pairs

We compute false discovery rates (FDR) between each cytokine pair to test the statistical significance of co-regulating target genes. For each target gene C, there are three types of co-regulations from a cytokine pair.

cytokine A and B both induce target gene C.

cytokine A and B both suppress target gene C.

cytokine A (or B) induces target C, but the other cytokine B (or A) represses target C.

First, we defined two log2-fold change (logFC) thresholds of cytokine A and B for the FDR computation. For type 1 coregulation (Co-Enhance), we computed the FDR (thres,thres) as The gene count derives directly from the data. The random count is equal to N as the total number of genes. We compute both probabilities from the logFC rank of each gene. In summary, the FDR computations are as follows:

Similarly, for type 2 coregulation (Co-Repress), we compute the FDR as

For type 3 coregulation (Antagonize), we compute the FDR as

After computing the FDR at each threshold combination (thres_A, thres_B), we adjust FDRs into monotonically decreasing values with respect to increasing threshold values, following the q-value procedure[66]. Finally, for the triplet of each cytokine pair and target gene, its statistical significance is the FDR (logFC_A, logFC_B) under each coregulation category.

Penalized Linear Model to Predict Cytokine Target Activities

The CytoSig linear model is programmed through a combination of Python 3 and GNU Compiler Collection (GCC) 4 C++. We only included 43 cytokines with at least five high-quality experiments (details in the Methods section Data Quality Control). We utilized a linear model to identify each signaling molecule’s signature patterns in an input sample’s expression profile. Composite profiles of cytokine response are the explanatory variables, and an input sample’s transcriptomic profile is the response variable. The regression coefficients represent cytokine target activities. The linear regression with all cytokine composite profiles as explanatory variables will reduce a cytokine’s coefficient if other cytokines with similar response profiles have more extensive impacts on the sample’s transcriptomic pattern[67].

The expression values, from either RNASeq or MicroArray, should be transformed by log2(x+1). We also recommend quantile-normalization across conditions. Some software packages, such as RMA or DESeq, will automatically include all normalizations. We recommend input differential profiles between the two conditions. If data is from a sample collection without pairs, please mean-centralize the value of each gene across all samples.

Many cytokine profiles are highly similar (Fig. 1f); such signature collinearity will create large result variance in a regular linear regression[65]. Therefore, we utilized the penalized ridge regression, which trades off the result bias to reduce the variance. The vector y is the input sample’s expression profile. The matrix X contains composite profiles of 43 cytokines. The parameter λ is the penalty. The ridge regression aims to minimize the objective function (y – Xβ)(y – Xβ) + λ * ββ. The coefficient βrepresents signaling activities.

To optimize parameters, we evaluated two types of model performance:

1, Prediction performance evaluates how the fitted model of cytokine activities predicts a sample’s gene expression profile. We use the five-fold cross-validation (CV) R2 ratio as the prediction performance metrics.

2, Inference performance evaluates whether coefficients on cytokine covariates of the fitted model represent the actual cytokine activities in a sample. We used the correlation between model coefficients and the ligand or receptor expression across samples as the inference performance metrics for each cytokine.

Typically, the training of ridge regression models only evaluates the prediction performance through cross-validation to determine the optimal penalization factor and coefficients[65]. However, we also evaluated the inference performance because the goal of CytoSig model is to infer cytokine signaling activities. The collinearity among cytokine response profiles may not affect the prediction performance but will induce significant variance on model coefficients[65], thus undermining the inference performance. The penalty factor in the ridge model will reduce the model variance at the cost of lowering prediction performance. Thus, we aim to find a penalty factor as a trade-off between two performance aspects.

On average, the CV R2 metric reaches its maximal point at a low penalty factor and deteriorate while the penalty factor is increasing (Extended Data Fig. 5a, b). In contrast, the inference performance, measured as correlation values, monotonically increases with increasing penalties (Extended Data Fig. 5c, d). Therefore, we select a value of 10000, which is the minimal lambda to achieve 80% best inference performance and 70% best prediction performance. Such a penalty will control both result variance and bias in the ridge regression.

We also evaluated XG Boost, a popular machine learning algorithm[68] (Extended Data Fig. 5). The XG boost with tree learners outperformed ridge regression in prediction but does not provide any coefficients on cytokine covariates for the inference purpose due to the tree structure of learners. The prediction performance of XG boost with linear learners quickly deteriorates to zero with increasing penalties although it has a high inference performance. Ridge regression is the only method with reasonable performance in both prediction and inference metrics.

We utilized a permutation test to estimate ridge coefficients’ standard errors after shuffling gene identities 1000 times. The z-scores (coefficient - random_average_coefficient) / standard_deviation on each cytokine represent its target activity.

T-cell activation and TGFB1 treatment assay

Human primary T cells are sourced from Hong Kong Red Cross Transfusion Service. Peripheral blood mononuclear cells (PBMC) were isolated from healthy donors using the Ficoll Paque Plus (GE healthcare # 17-1440-03) via density gradient centrifugation. CD8 T cells were purified from fresh PBMC by magnetic negative selection using the human CD8 T cell isolation kit (Miltenyi Biotec, Cat# 130-096-495). Isolated cells were stimulated with the human T cell TransAct (Miltenyi Biotec, Cat# 130-111-160) in the presence or absence of human recombinant TGFB1 (R&D systems, 240-B-002) at 5 ng/mL for 72 hours. Cells were cultured in MACS GMP medium, which is TexMACS GMP medium (Miltenyi Biotec, Cat# 170-076-309) supplemented with 10% inactivated fetal bovine serum (Gibco, 10082147), 50uM 2-mercaptoethanol (Gibco, 21985023), 10 mM N-Acetyl L-cysteine and 1% Penicillin-Streptomycin (P/S, Gibco, 15140122) at 1E6 cells/mL.

Mice CD8 T cells were isolated from splenocytes of one 8-week male C57BL/6J mouse using the CD8a T cell isolation kit (Miltenyi Biotec, Cat# 130-104-075) by magnetic negative selection. Isolated CD8+ T cells were stimulated with plate-bound anti-mouse CD3 (Biolegend, Cat# 100202, clone 17A2) at 5 ug/mL (1:100 dilution) and soluble anti-mouse CD28 (Biolegend Cat# 102102, clone 37.51) at 2 ug/mL (1:250 dilution) in the presence or absence of human recombinant TGFB1 (R&D systems, 240-B-002) at 5 ng/mL for 72 hours. Cells were cultured in complete RPMI 1640 medium, which is RPMI 1640 Medium (Gibco, 11875119) supplemented with 10% inactivated fetal bovine serum (Gibco, 10082147), 20 mM HEPES (Gibo, 15630080), 1 mM sodium pyruvate (Gibco, 11360070), 50uM 2-mercaptoethanol (Gibco, 21985023), 2 mM L-glutamine (Gibco, 25030024), and 1% Penicillin-Streptomycin (P/S, Gibco, 15140122) at 1E6 cells/mL.

Human inadequate whole blood was collected with informed consent and protocols approved by the ethics committee at the University of Hong Kong and the Hong Kong Red Cross Blood Transfusion Service. Animal experiments were approved by the committee of the Use of Live Animals in Teaching and Research at the University of Hong Kong (HKU) and performed strictly according to the animal protocol 5310–20. C57BL/6J mice were purchased from the Laboratory Animal Unit (LAU) of HKU.

BMP6 and IL1B in-vitro treatment combinations

NCI-H1299 (ATCC #CRL-5803) and A549 (ATCC #CCL-185) cells were purchased from American Type Culture Collection (ATCC). NCI-H1299 cells are cultured in high glucose DMEM medium (Gibco) supplemented with 10% Fetal Bovine Serum (FBS, Gibco BRL) and 100IU/mL penicillin/streptomycin (P/S). Human A549 cells are cultured in F12-K medium (ATCC, 30–2004) supplemented with 10% FBS and 100IU/mL penicillin/streptomycin (P/S).

NCI-H1299 and A549 cells were seeded in 6-well-plate at the density of 2E5 cells/well. On the next day, cells were treated with human recombinant IL1B (R&D systems, 201-LB-005, 10 ng/mL) alone or in combination with human recombinant BMP6 (R&D systems, 507-BP-020, 10 ng/mL) for 12 hours. In an alternative sequential treatment schedule, cells were pre-treated with IL1B first for 12 hours, then BMP6 or Media control for another 12 hours. Reconstitution buffers of the IL1B (PBS containing 0.1% BSA) and BMP6 (4 mM HCl containing 0.1% BSA) were used as negative controls.

Flow cytometry

For the evaluation of intracellular markers on A549 and H1299 cells, the following antibodies were used at indicated dilutions:

PE anti-human PRF1 (Biolegend, 353303, clone B-D48, 1:50 dilution),

PE anti-mouse PRF1 (Biolegend, 154305, clone S16009A, 1:50 dilution),

PE anti-human GZMA (Biolegend, 507206, clone CB9, 1:50 dilution),

PE anti-mouse GZMA (Biolegend, 149703, clone 3G8.5, 1:100 dilution),

FITC anti-human/mouse GZMB (Biolegend, 515403, clone GB11, 1:50 dilution),

APC anti-human/mouse MCP-1 (CCL2, Biolegend, 505909, clone 2H5, 1:200 dilution),

FITC anti-human CXCL8 (Biolegend, 511406, clone E8N1, 1:50 dilution).

Cells were fixed before permeabilization per the manufacturer’s instructions of wash buffer (Biolegend, 421002), and followed by intracellular staining with the antibodies. Flow cytometry was performed on ACEA NovoCyte Quanteon and raw data was analyzed using Flowjo software (Version 10.7).

To determine the gating threshold to detect marker positive cells, we used the forward scatter height (FSC-H) and side scatter height (SSC-H) for dead cell and debris removal. FSC-H/width and SSC-H/width are used to select single cells. We include unstained cells to define the threshold that separates positive populations from negative control cells (Extended Data Fig. 2d and Extended Data Fig. 4c).

CXCL8 and CCL2 detection by Enzyme-linked immunosorbent assay (ELISA)

A549 cells were seeded in 6-well-plate at the density of 2E5 cells/well. On the next day, cells were treated with human recombinant IL1B (R&D systems, 201-LB-005, 10 ng/mL), human recombinant BMP6 (R&D systems, 507-BP-020, 10 ng/mL), and combinations of IL1B and BMP6 for 24h. Reconstitution buffers of the IL1B (PBS containing 0.1% BSA) and BMP6 (4 mM HCl containing 0.1% BSA) were used as negative controls.

The amount of released CCL2 and CXCL8 from tumor cells in the supernatants was measured by ELISA assay using human CCL2 DuoSet ELISA kit (R&D systems, DY279) and human CCL2 DuoSet ELISA kit (R&D systems, DY208). Optical density (OD) value was determined using a microplate reader (TECAN, Infinite 200) at 450 nm wavelength with the correction wavelength set at 570 nm.

Upon treatment combinations of IL1B and BMP6 after 24 hours, supernatants from different conditions were 200X diluted and measured. The experiment was repeated independently in three batches. In each batch, a standard curve was created to measure the relationship between fluorescence values and seven 2X concentration dilutions from 2000pg/ml and 1000pg/ml for CXCL8 and CCL2, respectively. The In(concentration + 1) and fluorescence values follow a linear relationship. We fitted a linear regression model to convert the fluorescence measurements to concentrations.

Anti-TGFB Animal Studies

XOMA068 (pan-TGF-β1,2,3), XOMA089 (TGF-β1,2 selective), and anti-KLH (control) antibodies, supplied by XOMA Corp (Berkeley, CA), were all fully human IgG2(κ) antibodies generated by phage display and affinity maturation in our previous study[69]. Briefly, fully human antibody phage display libraries were used to discover a number of antibodies that bind and neutralize various combinations of TGFβ1, 2, or 3. The primary panning did not yield any uniformly potent pan-isoform neutralizing antibodies; therefore, an antibody that displayed potent TGFβ 1, 2 inhibition but more modest affinity versus TGFβ3, was affinity matured by shuffling with a light chain sub-library and further screening. This process yielded the high-affinity pan-isoform neutralizing clone. Antibodies were diluted in 10mM histidine, 142 mM L-arginine, pH6.0 buffer “vehicle” for in vivo studies.

Animal studies were conducted under protocol LC-070, approved by the Animal Care and Use Committee of the National Cancer Institute (NCI), Bethesda, MD. The animals are on a 12:12 light: dark cycle. The ambient temperature is 72F +/− 2F, and the humidity is kept between 30–70%. 40,000 4T1 mouse mammary tumor cells were surgically implanted into four mammary fat pads of 8 week old female BALB/c mice. From day one after surgery, mice were treated with anti-TGFB antibodies at 5mg/kg intraperitoneally three times per week for two weeks. Tumors were surgically resected on day 13 when they reached 0.8–1cm in diameter and were snap-frozen for molecular analysis. The NCI ethics committee requires that animals must be euthanized at the time of observation if the tumor size is approaching 20 mm, in any dimension. None of the tumors in our experiment exceeded this limitation.

RNA was isolated from tumor samples using the RNeasy method (Qiagen) according to the manufacturer’s instructions following tissue lysis with a Precellys 24 Homogenizer (Bertin Instruments). Tumor RNA that passed quality control (RIN>7) was sequenced on HighSeq2500 using Illumina TruSeq v4 chemistry, generating 50–100 million pass-filtered reads per sample. There are 6 mice in the XOMA089 group, 7 mice in the XOMA068 group and anti-KLH group. No data points were excluded from the analysis.

Identification of Signaling Signatures in COVID-19 Severe Symptoms

For each single-cell dataset, we computed the cytokine activities for individual cells using the CytoSig model and got the mean value for each cell type in each patient. Then, for each cell type, we compared activities between different patient groups using the two-sided Wilcoxon rank-sum test and converted the p-values to false discovery rates (FDR) by the Benjamini-Hochberg correction. FDR 0.05 is the threshold for the result significance.

For the COVID-19 study on bronchoalveolar lavage[45] and peripheral blood[46] samples, we performed comparisons between severe and mild patients, and between disease (severe+mild) and healthy individuals. Our analysis only reported results identified in both comparisons with FDR less than 0.05. We made an exception for the analysis of neutrophils. Among neutrophils from peripheral blood, a few cytokines’ signaling differences between severe and mild patients achieved a statistical significance of FDR 0.051. We believe these results are still significant, thus reported them in our analysis.

We used the original coordinates of two-dimensional embedding from each publication (plotted in Fig. 6b, d). The bronchoalveolar lavage study[45] utilized t-SNE (t-distributed stochastic neighbor embedding) that projects the single-cell RNASeq profiles in two dimensions with distances between dots representing the profile similarities. The peripheral blood study[46] utilized the UMAP (Uniform Manifold Approximation and Projection), a dimension reduction approach.

Statistics & Reproducibility

All comparisons between two groups utilized the two-sided Wilcoxon rank-sum test, a non-parametric test without any assumptions on the data distribution. Similarly, all comparisons between group values and zero used the non-parametric Wilcoxon signed-rank test. No data were excluded from any analyses.

No statistical method was used to predetermine sample size. Instead, we selected a fixed sample size in the following experiments. In the in-vitro validation of TGFB1’s inhibitory role (Fig. 2e) and BMP6’s anti-inflammatory role (Extended Data Fig. 4), we utilized a sample size of three, the minimum number to achieve statistical significance of p-value <= 0.05 in the two-sided Wilcoxon rank-sum test. All cell-culture replicates lead to reproducible successful results (Fig. 2e and Extended Data Fig. 4b, d). In the TGF-beta blocking in-vivo experiment (Fig. 4d), we used a minimal mouse number of six, suggested by a previous study to detect differential expression events through RNA-Seq[70]. Our recent study demonstrated that four tumors (smaller than our sample size 6) per condition would be sufficient to detect differentially expressed genes between conditions[71].

Mouse identities were randomized before in-vivo experiments. Randomizations were not performed for in-vitro cell cultures because all conditions were derived from a homogeneous cell-line population. Blinding is not performed in our experiments because the robust phenotype of our results is based on strictly objective measurements by equipment instead of any human estimations. The outcome assessments include flow cytometry in Fig. 2e and Extended Data Fig. 4a–c, ELISA assay plate reader in Extended Data Fig. 4d, and RNA-sequencing in Fig. 4d. None of these measurements involve human subjective perception.