Tang et al.

Cell Discovery (2024)10:39

https://doi.org/10.1038/s41421-024-00650-7
Cell Discovery
www.nature.com/celldisc

ARTICLE Open Access

Integration of 3D bioprinting and multi-algorithm

machine learning identified glioma susceptibilities
and microenvironment characteristics
Min Tang1,2 ✉, Shan Jiang3, Xiaoming Huang4,5,6,7,8, Chunxia Ji4,5,6,7,8, Yexin Gu9, Ying Qi4,5,6,7,8, Yi Xiang2, Emmie Yao2,
Nancy Zhang10, Emma Berman11, Di Yu10, Yunjia Qu12, Longwei Liu 13, David Berry2,14 and Yu Yao4,5,6,7,8 ✉

Abstract
Glioma, with its heterogeneous microenvironments and genetic subtypes, presents substantial challenges for
treatment prediction and development. We integrated 3D bioprinting and multi-algorithm machine learning as a
novel approach to enhance the assessment and understanding of glioma treatment responses and microenvironment
characteristics. The bioprinted patient-derived glioma tissues successfully recapitulated molecular properties and drug
responses of native tumors. We then developed GlioML, a machine learning workflow incorporating nine distinct
algorithms and a weighted ensemble model that generated robust gene expression-based predictors, each reflecting
the diverse action mechanisms of various compounds and drugs. The ensemble model superseded the performance
of all individual algorithms across diverse in vitro systems, including sphere cultures, complex 3D bioprinted
multicellular models, and 3D patient-derived tissues. By integrating bioprinting, the evaluative scope of the treatment
expanded to T cell-related therapy and anti-angiogenesis targeted therapy. We identified promising compounds and
1234567890():,;
1234567890():,;
1234567890():,;
1234567890():,;

drugs for glioma treatment and revealed distinct immunosuppressive or angiogenic myeloid-infiltrated tumor
microenvironments. These insights pave the way for enhanced therapeutic development for glioma and potentially
for other cancers, highlighting the broad application potential of this integrative and translational approach.

Introduction undergoing radiotherapy and temozolomide (TMZ)

Glioma is a complex central nervous system cancer maintenance treatment2. However, mismatch repair
exhibiting significant genetic and phenotypic hetero- defects and hypermutation following the TMZ treatment
geneity among patients1. Investigations focusing solely on contribute to drug resistance, and GBM patients often
molecular alterations driving neoplastic events or pre- experience high relapse rates and develop chemoresis-
clinical studies have translated into limited survival tance3. Immunotherapies that demonstrate promising
advantages in clinical practice for glioblastoma (GBM), effects in other cancers also show varied efficacies in
the most lethal type of glioma with a five-year survival rate GBM4–6. Despite extensive attempts in glioma investiga-
of 6.9%1. Methylation of the O-6-methylguanine-DNA tions and drug development, therapeutic improvement
methyltransferase (MGMT) promoter has been associated appears stagnant, necessitating a further understanding of
with better prognosis in newly diagnosed GBM patients the disease and more accurate treatment evaluations to
optimize patient survival rates7.
Growing evidence highlights that complex tumor
Correspondence: Min Tang ([email protected]) or microenvironments, including myeloid cells, significantly
Yu Yao ([email protected])
1
impact glioma treatment responses8,9. Clinical studies
Shanghai University of Traditional Chinese Medicine, Shanghai, China
2 have shown that drugs demonstrating promising out-
Department of NanoEngineering, University of California San Diego, La Jolla,
CA, USA comes in preclinical settings did not yield significant
Full list of author information is available at the end of the article

© The Author(s) 2024

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction
in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if
changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If
material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain
permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Tang et al. Cell Discovery (2024)10:39 Page 2 of 20

clinical survival benefits, suggesting that current pre- gene expression data, validated using a GBM cell line,
clinical models are suboptimal6,10. Traditional 2D cell bioprinted PDTs, and bioprinted multicellular GBM-
culture fails to recapitulate the heterogeneity of cells and myeloid models. GlioML unveiled unique susceptibility
the extracellular matrix (ECM) of the tumor micro- patterns of glioma patients to a range of compounds,
environment in vivo. In GBM, tumor-associated macro- several of which were confirmed by the bioprinted patient
phages/microglia (TAMs), a highly plastic population of tissues, including lovastatin, dasatinib, and 1S3R-RSL-3
non-neoplastic cells, constitute a substantial part of the (RSL), each operating through distinct mechanisms. Fur-
tumor mass and are implicated in GBM malignancy and ther, this integration discerned discrete angiogenic or
drug resistance11–13. TAMs, originating from two distinct immunosuppressive phenotypic traits in patient-derived
myeloid sources, blood-circulating monocytes, and brain- models and GBM-myeloid models developed using either
resident microglia, form a pro-tumor stroma for GBM microglia or monocytes. Testing of immunotherapy and
growth and progression12. Creating clinically relevant targeted therapy, in conjunction with cytokine profiling
GBM microenvironments requires precise control of and immunofluorescent staining of 3D bioprinted GBM
cellular compositions and ECM properties, have been models, further reinforced GlioML’s discoveries. Our
achieved with recent advancements in 3D bioprinting investigation underscores the synergistic capacity of this
techniques and tissue-relevant biomaterials14,15. The approach to generate clinically pertinent treatment eva-
enabled crosstalk between tumor cells and stromal cells luations for patients and enrich our comprehension of
provides valuable insights into tumor cell dependencies distinct myeloid characteristics in the tumor micro-
and phenotypic features16. However, bioprinting has not environment. This integrative strategy holds potential for
been exploited to develop patient-derived models using further refinement and broad application across diverse
freshly isolated tumor tissues, which might better reca- cancer types, heralding a transformative era in cancer
pitulate the original tissues. research and treatment.
In addition, the advancement of machine learning has
facilitated better pattern recognition from large and Results
complex datasets, such as next-generation sequencing 3D bioprinted PDTs recapitulated genetic features and
data. Machine learning has been successfully deployed in clinical drug responses of patient tumors
various disease management processes and to predict Surgically resected GBM tumor tissues were obtained
cancer drug responses by integrating multi-omics fea- from glioma patients spanning varied demographics
tures, such as gene expression8,17–20. However, current (Supplementary Fig. S1a). PDT cultures were successfully
workflows face several challenges, including potential generated for all obtained patient specimens as drug
overfitting issues due to the abundance of omics data and testing models for 22 adult and 1 pediatric high-grade
relatively scarce drug response data, emphasizing the glioma (HGG) patients, encompassing GBM, astro-
critical importance of feature engineering. Furthermore, cytoma, oligodendroglioma, and pediatric HGG. Patient
current machine learning workflows usually rely on a tumor tissues were rinsed and dissociated using col-
single algorithm for prediction, which is potentially lagenase to generate single-cell solutions. The bioink
unsuitable for different drug compounds operating comprised of gelatin methacrylate (GelMA) and hya-
through distinct mechanisms. luronic acid methacrylate (HAMA), has been formulated
In this study, we reported the first integration of 3D to mimic the glioma ECM with greater fidelity than
bioprinting and machine learning to encompass diverse Matrigel commonly used in organoid cultures16. The
molecular and cellular features to enhance the assessment prepolymer solution and the single-cell suspension were
of glioma treatment responses and microenvironment thoroughly mixed to generate bioink for individual
characteristics. We first leveraged bioprinting to create patients and loaded onto the bioprinter, which had 96
biomimetic, patient-derived tissues (PDTs) that closely independent light sources at the wavelength of 405 nm, to
replicate the original tumor’s genomic characteristics, polymerize the bioink and create PDTs.
gene expression patterns, and cellular compositions, Pairs of patient specimens and their PDT counterparts
maintaining high clinical relevance. The bioprinted were evaluated via RNA sequencing (RNA-seq) and whole
models generated with recurrent patient tissues exhibited exosome sequencing (WES) to determine their molecular
drug resistance to TMZ treatment, which the MGMT statuses (Fig. 1a). RNA-seq showed a high correlation
promoter (pMGMT) methylation assessment failed to between PDTs and their corresponding patient tissues
indicate. We next developed GlioML, a multi-algorithm (Fig. 1b). Gene expression patterns remained mostly
machine learning strategy, to generate independent drug consistent between the original tumor and its PDT
predictors with nine stack-one classic single algorithms counterpart, although variations were observed among
and one stack-two weighted ensemble model. GlioML individuals (Fig. 1c). Principal component analysis (PCA)
produced reliable drug response predictions based on indicates that most tissue samples and corresponding
Tang et al. Cell Discovery (2024)10:39 Page 3 of 20

Fig. 1 Bioprinting of clinically relevant PDTs for drug evaluations. a Schematic of processing workflow for patient specimens, including the
generation of PDTs via bioprinting for drug testing and characterization using RNA-seq, WES, and flow cytometry for both primary tissues and PDTs.
b Pearson correlation graph of the log-transformed gene expression data between primary tissues and PDTs. r: correlation coefficient. c Heatmap
representation of the transcriptome comparison between three pairs of PDTs and their respective tissue samples, demonstrating high similarity,
including a primary GBM patient (4358-HS), a recurrent astrocytoma patient (5256-HS), and a recurrent GBM patient (4089-HS). d Venn diagrams
illustrating genomic concordance between two patient samples and their PDTs, highlighting over 90% overlap in detected single nucleotide variants.
e Concordance of the recurrent clinical status to TMZ sensitivity predicted by pMGMT methylation status or assessed through bioprinted PDTs.
f Comparative efficacy of tumor cell inhibition for TMZ and CCNU in PDTs, showcasing the superior efficacy of CCNU across all patients, especially in
recurrent cases.
Tang et al. Cell Discovery (2024)10:39 Page 4 of 20

PDTs cluster closely, except for an oligodendroglioma, a patients after the CCNU treatment was 68%, compared to
gliosarcoma (GSM), and one GBM specimen showing 115% for TMZ (P < 0.001). PDT responses to clinical
minor deviations (Supplementary Fig. S1b, c). A com- drugs aligned with a meta-analysis identifying CCNU as
parative analysis of RNA-seq data from patient tissues, the most effective chemotherapy for recurrent GBM
bioprinted PDTs, and Matrigel-based patient-derived patients21.
organoids (PDOs) revealed good correlations between
both in vitro models and patient tissues (Supplementary Multi-algorithm workflow identified optimal algorithms for
Fig. S2a). The bioprinted models demonstrate a margin- diverse compounds
ally higher correlation coefficient (0.78 and 0.86 for PDT, We extracted the gene expression data of established
compared to 0.72 and 0.72 for PDO). This might indicate cancer cell lines, as profiled by RNA-seq from the Cancer
a better representation of the in vivo condition by the Cell Line Encyclopedia (CCLE), and the drug response
bioprinted PDTs. Despite identical cellular encapsulation, data from the Cancer Therapeutics Response Portal
the observed differences may stem from the effects of the (CTRP), which included 481 drugs and their corre-
different ECM compositions and properties. Furthermore, sponding drug sensitivity results in well-characterized
the majority of mutations identified in patient tissues were cancer cell lines as area-under-curve (AUC) values23,24.
mirrored in PDT cultures, confirming the genomic fidelity We curated a list of gene sets most relevant to our study,
of the bioprinted PDTs to the original tumors (Fig. 1d). covering canonical pathways, transcription factor-binding
Copy number variations were largely retained across both sites, and oncogenic signature gene sets, as the initial
PDTs and PDOs (Supplementary Fig. S2b). Regarding features for GlioML, the multi-algorithm auto-machine
single nucleotide variants (SNVs) as well as short inser- learning workflow25,26. Expression data of raw protein-
tions and deletions (INDELs), both in vitro models suc- coding genes were converted into single-sample gene set
cessfully conserved the mutational landscape observed in enrichment analysis (ssGSEA) scores of the curated gene
patient tissues (Supplementary Fig. S2c). Nonetheless, sets. These ssGSEA scores provided insight into the level
PDOs exhibited an increased incidence of de novo of regulation within each gene set and offered more
mutations, which was suboptimal for replicating native interpretability than individual gene expression. Of all cell
tissue characteristics. lines in CTRP, 636 were comprehensively characterized in
The drug testing involved only first-passage PDTs the CCLE and therefore extracted for further use. We
derived from freshly dissociated cells to maximize clinical employed the AUC values for each cell–compound pair-
relevance. All PDTs were subjected to the current gold ing as labels, while pairs with missing labels were pre-
standard treatment, TMZ, and lomustine (CCNU), a processed and subsequently excluded. We consolidated
superior option for recurrent GBM21. Additionally, two all this information to generate the raw training data for
platinum-based compounds were tested due to their fre- the workflow (Fig. 2). To mitigate the batch effects in
quent use in pediatric HGG (pHGG) treatment. Viability RNA-seq data resulting from disparate sample handling
assessments were conducted post 48 h treatment with processes, like cell culture methods and library prepara-
100 μM of the compounds (Supplementary Fig. S3). tions, we performed normalization and data cleansing on
Lobaplatin was the only drug among the four tested to the features of both training data and our sample data
have both its median and average tumor inhibition rates before submitting them to GlioML (Supplementary
over 50% in PDT models. PDTs more accurately predicted Fig. S4).
TMZ resistance in most recurrent adult patients who had We incorporated nine supervised base algorithms in our
already undergone TMZ treatment, compared to the workflow, including three variations of Light Gradient-
pMGMT methylation status (Fig. 1e). Although pMGMT Boosting Machine (LightGBM), two of K-Nearest
status (> 10% clinical threshold) indicated TMZ sensitivity Neighbor (KNN), FastAI Neural Network (FastAI),
in 6 out of 9 adult recurrent patients22, drug testing using NeuralNetTorch (NNTorch), Random Forest, and Extra
PDTs revealed that 5 of these 6 patients would not Trees, along with a stack-two weighted ensemble
respond to TMZ treatment. This was evidenced by PDTs model27–32. We chose algorithms representing a broad
demonstrating more than 100% cell viability post- range of learning categories. Each compound underwent
treatment (Table 1). The concordance rate of PDTs in independent training to account for its unique mechan-
reflecting drug resistance among recurrent patients was isms of action and to accommodate the potential for
89%, in contrast to only 33% for pMGMT status. CCNU different algorithms to yield optimal performance. We
exhibited a higher efficacy in inhibiting tumor cells in evaluated model performance via validation scores, with
PDTs, particularly in recurrent GBM patients (Fig. 1f). lower absolute values signifying better performance. The
The median viability of all PDTs after the CCNU treat- tabular data were restructured to incorporate the top n%
ment was 74%, compared to > 100% for TMZ. Moreover, features, and the models were re-trained to improve
the median viability of PDTs generated from recurrent performance further.
Table 1 Individual patient clinical information, pMGMT status, and PDT drug responses to clinical drugs TMZ and CCNU.

ID Age Sex Histologic diagnosis TMZ Treated pMGMT (%) TMZ (%)a CCNU (%)b

2509-HS 43 M Astrocytoma, IV, IDH1 mut (R) Yes 12.0c 121 -

2745-HS 51 M Oligodendroglioma, III, IDH1 mut (R) Yes 31.2c 138 98
c
2766-HS 41 M Oligodendroglioma, III, IDH1 mut No 39.5 106 70
c
3056-HS 35 F Astrocytoma, IV, IDH1 mut (R) Yes 30.7 124 61
c
3199-HS 27 F Giant cell glioblastoma, IV, IDH1/2 WT No 49.0 93 44
3833-HS 21 M Pediatric High-Grade Glioma, IV, IDH1/2 WT (R) No 1.7 99 101
Tang et al. Cell Discovery (2024)10:39

3917-HS 43 F Oligodendroglioma, III (R) Yes 19.0c 114 31

c
4089-HS 45 M Glioblastoma, IV, IDH1/2 WT (R) Yes 38.2 112 57
4240-HS 37 M Gliosarcoma, IV, IDH1/2 WT (R) Yes - 130 39
4358-HS 67 M Glioblastoma, IV, IDH1/2 WT No 1.5 81 82
4468-HS 64 M Glioblastoma, IV, IDH1/2 WT No 1.5 92 36
4671-HS 55 F Glioblastoma, IV, IDH1/2 WT No 7.0 62 29
5256-HS 31 F Astrocytoma, IV, IDH1 mut (R) Yes 47.3c 92 74
0810-HS 66 M Glioblastoma, IV, IDH1/2 WT No 2.0 91 83
1779-HS 52 M Glioblastoma, IV, IDH1/2 WT No 2.25 106 92
c
1869-HS 74 M Glioblastoma, IV, IDH1/2 WT No 29.7 - -
1864-HS 65 M Glioblastoma, IV, IDH1/2 WT No 2.75 106 101
10308-HS 39 M Gliosarcoma, IV, IDH1/2 WT (R) Yes 1.5 104 90
c
0911-HS 58 F Glioblastoma, IV, IDH1/2 WT No 33 84 113
27691-HS 70 F Glioblastoma, IV, IDH1/2 WT No 2.0 102 84
28372-HS 32 F Astrocytoma, IV, IDH1/2 mut (R) Yes 3.0 106 96
29259-HS 67 F Glioblastoma, IV, IDH1/2 WT No - 101 94
c
5371-HS 42 F Astrocytoma, III, IDH1 mut No 21.2 99 102
a
Viability of PDT after TMZ treatment.
b
Viability of PDT after CCNU treatment.
c
pMGMT methylation status indicating sensitivity to TMZ treatment.
Page 5 of 20
Tang et al. Cell Discovery (2024)10:39 Page 6 of 20

Fig. 2 GlioML, a multi-algorithm machine learning approach for drug response prediction. The process consists of four key components: (1)
Data Preprocessing — standardization and labeling of the training dataset, (2) Feature Engineering — derivation of relevant features for robust
glioma drug response prediction, (3) Training Module — a comprehensive process involving initial training, feature reduction, and performance
assessment for individual predictor development, and (4) Downstream Applications — functionalities such as drug response predictions and
microenvironment interpretations.

The training results underscored the superior predictive range. Half maximal inhibitory concentration (IC50)
capacity of neural networks and gradient-boosting algo- values for each drug–compound pair were obtained.
rithms in estimating drug responses from gene expression A second round of feature engineering was conducted to
features. These two categories produced over 98% of the prevent potential overfitting due to the original feature-to-
best predictors from single algorithms, while the KNN sample size ratio, reducing the number of features used in
models did not generate any top predictors (Fig. 3a). training to a quantity less than the total sample number.
Moreover, through an effective combination of base The top 20%, 10%, and 3% highest-ranked features were
models with optimized weights (Supplementary Fig. S5), selected for further training, which resulted in approxi-
the weighted ensemble model outperformed all single mately 1:1, 2:1, and 6:1 sample-to-feature size ratios,
algorithms across all compounds in the training datasets respectively. Reducing the feature size from the original to
(Supplementary Fig. S6). the top 20% (RF20) significantly enhanced the validation
scores for most drugs, except for lovastatin (P = 0.03) (Fig.
GlioML predicted drug response on an independent 3b). The validation scores for six of the nine compounds
patient-derived dataset improved when the feature size was reduced to 3% (RF3)
We next used a line of GBM stem cells (GSCs), CW468, compared to 10% (RF10). However, dasatinib, docetaxel,
to assess the predictive capacity of the workflow on and trametinib experienced diminished performance, sug-
independently generated data. CW468 was profiled by gesting that a 3% cutoff may be excessively stringent for a
RNA-seq, and its gene expression data were processed few compounds. Statistical analysis of validation scores
and analyzed. AUC predictions were obtained for all revealed no significant difference between RF20, RF10, and
compounds from every algorithm and the weighted RF3. Thus, decreasing the feature size to a ratio less than 1:1
ensemble model. The GSCs were then exposed to nine with the sample size was critical, and additional fine-tuning
compounds, including TMZ, the standard drug treatment yielded minor improvements. RF3-trained predictors were
for GBM, with the predicted AUC spanning the entire utilized for downstream investigations.
Tang et al. Cell Discovery (2024)10:39 Page 7 of 20

Fig. 3 Performance evaluation of GlioML for drug response predictions. a The proportion of top-performing predictors generated by different
single machine learning algorithms. b Absolute validation scores for the weighted ensemble model with feature reduction applied to the top 20%,
10%, and 3% of highest-ranked features. Lower scores imply better performance. c Linear regression analysis comparing log-transformed IC50 values
in GSC with predicted drug response AUC by the weighted ensemble model across all four conditions: original, RF20, RF10, and RF3. d–l Individual
linear regression analyses of log-transformed IC50 values and predicted AUC by individual algorithm models: NNTorch (d), FastAI (e), LightGBM (f),
LightGBMXT (g), LGBMLarge (h), Random Forest (i), Extra Trees (j), KNNUinf (k), and KNNDist (l). R2: coefficient of determination.

A strong linear relationship was observed between the linear relationships between the log-transformed IC50
log2 transformed IC50 values in the GSCs and the pre- values and the predicted AUC (Fig. 3d–g). The other
dicted AUC from the weighted ensemble model across all models, including LightGBMLarge (LGBMLarge), Ran-
four conditions, including the original, RF20, RF10, and dom Forest, Extra Trees, KNN uniform (KNNUnif), and
RF3. Simple linear regression analyses on the measure- KNN distance (KNNDist), demonstrated less optimal
ments and predictions produced an R2 value exceeding performance on the GSCs (Fig. 3h–l). These results
0.80 in all groups, with RF3 predictions providing the best indicated that GlioML could be reliably employed for
fit (Fig. 3c). In the single algorithm models, NNTorch, independent datasets and thus can easily accommodate
FastAI, LightGBM and LightGBMXT, exhibited strong other cancer types. Furthermore, it was found that the
Tang et al. Cell Discovery (2024)10:39 Page 8 of 20

weighted ensemble, neural network, and certain variations Nevertheless, upon evaluation with the CW468 dataset,
of LightGBM models provided the most robust predic- the Glioma+ workflow demonstrated reduced perfor-
tions for the evaluated compounds compared to other mance compared to the original GlioML (Supplementary
models, corroborating the training set evaluation. Fig. S7d). This discrepancy of performance in the training
To further enhance the predictive precision of our dataset and independent dataset suggested that the
multi-algorithm workflow, we employed three strategies: improved validation scores may be attributed to over-
Glioma+ (focusing the dataset solely on glioma cases), fitting within a constrained training dataset, which in turn
Metadata+ (adding clinical metadata such as age, gender, could diminish the generalization ability and this was a
and pivotal mutation statuses), and Feature+ (including a crucial aspect of the workflow where predictive accuracy
more extensive array of non-linear conditions). Our ana- across diverse datasets was imperative. Consequently, we
lysis revealed that the Glioma+ approach significantly proceeded with the original GlioML workflow for further
increased accuracy within the training dataset, as indi- investigation.
cated by improved validation scores (Supplementary Fig.
S7a). The Glioma+ model attained a median validation GlioML identified promising compounds and predicted
score of –0.636, markedly better than –1.16 of the base drug efficacies on glioma PDTs
model and –1.14 of the Metadata+ and Feature+ models. We then explored the synergistic potential of integrating
Given the considerable extension in training time and the 3D bioprinted models and GlioML for predicting and
marginal score enhancement from integrating more evaluating individual glioma patient drug responses (Fig.
metadata or features, coupled with the need for time 4a). PCA of the AUC matrix for all patients using GlioML
efficiency, we proceeded with a trial training and evalua- effectively distinguished between World Health Organi-
tion using the Glioma+ condition. This strategic narrow- zation (WHO) grade III and WHO grade IV gliomas,
ing of the dataset also allowed us to incorporate additional suggesting that patients at different stages of glioma
base algorithms, specifically CatBoost and XGBoost, exhibited different drug response patterns (Supplementary
additional leading Boosting algorithms other than the Fig. S8a). Despite the limited number of cells extracted
LightGBM. The remaining methodology paralleled our from primary tissues, several PDTs were exposed to three
preceding multi-algorithmic approach. In the Glioma+ GlioML compounds, along with TMZ. These compounds,
workflow, neural networks and gradient-boosting models including RSL, dasatinib, and lovastatin, demonstrated
still emerged as the most effective, yielding ~97% of the favorable responses in the CW468. Following treatments
top-performing predictors (Supplementary Fig. S7b). The with RSL, dasatinib, and lovastatin, the median viabilities
superior performance of these algorithms compared to were 5.8%, 4.2%, and 50%, respectively. These viabilities
regression or bagging algorithms can likely be ascribed to significantly surpassed those achieved with TMZ. RSL and
their capacity for modeling the intricate patterns present dasatinib, lovastatin, and TMZ emerged in three different
in omic data. Additionally, techniques such as shrinkage clusters based on GlioML’s overall predictions for all
and regularization were pivotal in preventing overfitting, a patient samples (Fig. 4b). Given the effectiveness of RSL
critical concern in datasets where features often out- and dasatinib, the compounds within this cluster were
number samples, as observed in this study. In contrast, identified as potential candidates for glioma treatment and
bagging models, including Random Forests, though gen- could be further explored in the future (Table 2).
erally robust, may fall short of fully capturing the nuanced Although the efficiency of tumor-killing varied among the
complexity of omics datasets and managing the prevalent tested compounds, PDTs treated with CCNU, cisplatin,
noise and outliers. Notably, within the spectrum of lobaplatin, dasatinib, lovastatin, and RSL all displayed
boosting algorithms, LightGBM and XGBoost, which significantly lower tumor viability than untreated controls
utilize asymmetric tree structures, significantly out- (Supplementary Fig. S8b). As PDTs were shown to align
performed CatBoost, which employs symmetric trees. with clinical TMZ and CCNU responses, the superior
When considering the efficiency of neural network models tumor-killing efficacy of GlioML-identified compounds in
compared to gradient-boosting algorithms, the sig- PDTs supports the potential of GlioML to have a clinical
nificantly shorter median fit times for neural network impact on the treatment of glioma.
models such as FastAI and NNTorch (5 s and 46 s, Subsequently, we evaluated GlioML’s capacity to rank
respectively) suggest a computational advantage over the compound efficacy for individual patients by comparing
boosting algorithms, with CatBoost taking 246 s, the measured PDT viability to GlioML predictions.
LightGBM Large taking 167 s, and XGBoost taking 85 s GlioML displayed strong predictive potential for various
for median fit time (Supplementary Fig. S7c). This effi- WHO grade IV gliomas, including GBM, giant cell GBM,
ciency does not seem to come at the cost of performance, IDH mutant astrocytoma, and pHGG. However, it
as the neural network models still deliver a comparable demonstrated suboptimal performance for WHO grade
number of top-performance predictors. III grade oligodendroglioma. Notably, the PDTs treated
Tang et al. Cell Discovery (2024)10:39 Page 9 of 20

Fig. 4 Integration of bioprinted PDTs and GlioML for precision medicine. a Schematic of an integrative workflow that combines a
computational component (gene expression data acquisition and GlioML prediction) with an experimental part (PDT generation and multiple
compound/drug evaluations). b t-SNE plot representing all 481 compounds based on their response profiles to patient samples, with RSL and
dasatinib in the same cluster, while lovastatin and TMZ belong to different clusters. c Linear regression analysis correlating GlioML prediction and log-
transformed cell viability in five cases of high-grade gliomas, including three WHO grade IV adult gliomas, one WHO grade III oligodendroglioma, and
one pHGG. R2: coefficient of determination.

with all four GlioML compounds, giant cell GBM (3199- The oligodendroglioma (3917-HS) yielded a lower R2
HS), recurrent pHGG (3833-HS), recurrent GBM (4089- value of 0.57. The median R2 values were 0.83 for all
HS), and recurrent astrocytoma (5256-HS), yielded R2 WHO grade IV gliomas, suggesting that GlioML is par-
values of 0.83, 0.83, 0.94, and 0.84, respectively (Fig. 4c). ticularly adept at interpreting WHO grade IV gliomas.
Tang et al. Cell Discovery (2024)10:39 Page 10 of 20

Table 2 Promising compounds for glioma treatment based on GlioML predictions for drug response in patient samples.

Bardoxolone methyl Gemcitabine Topotecan Doxorubicin Teniposide

Foretinib NVP-BEZ235 AZD8055 NVP-TAE684 Cucurbitacin I
KPT185 LBH-589 Belinostat Apicidin BRD-A86708339
Obatoclax BRD-K63431240 Pluripotin Dasatinib Afatinib
Trametinib SCH-79797 Methotrexate Daporinad CAY10618
PF-3758309 Omacetaxine mepesuccinate SNS-032 Leptomycin B Narciclasine
SR-II-138A CR-1-31B Dinaciclib Alvocidib AZD7762
AT13387 SNX-2112 BI-2536 GSK461364 Triazolothiadiazine
Vincristine Parbendazole KX2-391 Paclitaxel SB-743921
Rigosertib Docetaxel Brefeldin A Ouabain MLN2238
Avicin D ML210 ML162 1S3R-RSL-3 Neopeltolide
Oligomycin A

Evaluating multimodal therapies through bioprinted PDTs HS). These PDTs were treated with T cells activated by
and multicellular models CD3/CD28 microbeads, at an effector:target (E:T) ratio of
HGGs, particularly GBM, exhibit heterogenous non- 1:1. Initially, 1779-HS comprised 41% CD14+ cells within
neoplastic populations, such as immune cells and endo- 45% of the total CD45+ cell population, whereas 29259-
thelial cells. Patient samples collected in this study HS contained 9% P2RY12+ cells within 9% of the total
exhibited a varied composition of stromal cells, especially CD45+ cells (Fig. 5d). Upon exposure to activated T cells,
the CD45+ immune cells, spanning from < 1% to 44%, and we observed a higher level of apoptosis, as indicated by
the PDTs effectively preserved the stromal cells (Supple- Caspase 3 staining (Fig. 5e), in P2RY12-dominant PDTs
mentary Fig. S9). Comparative analysis between 12 mat- compared to CD14-dominant PDTs, suggesting greater
ched pairs of patient tissues and PDTs showed no susceptibility to T-cell mediated cytotoxicity. Differential
significant variance (P = 0.43), barring a few instances responses were also noted when PDTs were exposed to
where a reduction in CD45+ ratio was observed (Fig. 5a). Bevacizumab, an anti-angiogenic drug. In CD14-
The macrophage populations, including subsets such as dominant PDTs, post-treatment CD31+ cells showed
CD14+ and P2RY12+ cells, which are indicative of TAMs over 50% inhibition relative to the baseline ratio, whereas
of peripheral and cerebral origin, respectively, were found the untreated control maintained the initial CD31+ cell
to constitute a considerable fraction of the tumor mass ratio. Conversely, P2RY12-dominant PDTs exhibited a
and CD45+ cell population, consistent with literature more resistant profile to Bevacizumab, with only a 5%
findings. While the CD14+ cells were well-maintained, inhibition in CD31+ cells post-treatment compared to
the P2RY12+ cells presented challenges in terms of the baseline. Intriguingly, the untreated P2RY12-
extended in vitro preservation (Fig. 5b). Among the dominant PDTs demonstrated a fourfold increase in the
examined patient tissues, CD3+ T cells accounted for less CD31+ cell ratio, suggesting an environment conducive
than 5% of the cell population in all samples, aligning with to angiogenesis.
the immunosuppressive microenvironment of HGGs. Despite their significant clinical relevance, PDTs
Varying amounts of CD31+ endothelial cells were exhibited patient-specific characteristics and presented
observed in patient samples, spanning from < 1% to 28%. inconsistent cellular compositions. This included varying
PDTs also successfully preserved the CD31+ cells, with no quantities of immune cells and cells related to blood
significant differences (P = 0.60) observed between 12 vessels, posing challenges for repeatable mechanistic
pairs of patient tissues and the corresponding PDTs studies. To address this, we employed digital light
(Fig. 5c). processing-based bioprinting to create engineered multi-
We further investigated the utility of bioprinted PDTs cellular GBM models to further investigate the differential
in assessing the efficacy of activated T cells and Bev- behaviors observed in CD14-dominant and P2RY12-
acizumab, in addition to the demonstrated capability of dominant PDTs. We attempted to reconstruct multi-
evaluating small molecule compounds. We selected cellular models that accurately represented the three cell
PDTs representing two distinct subtypes: CD14- populations involved in PDT examinations. These inclu-
dominant (1779-HS) and P2RY12-dominant (29259- ded CD14+ cells representing peripheral-origin TAMs,
Tang et al. Cell Discovery (2024)10:39 Page 11 of 20

Fig. 5 (See legend on next page.)

Tang et al. Cell Discovery (2024)10:39 Page 12 of 20

(see figure on previous page)

Fig. 5 Assessing multimodal treatments in bioprinted PDTs and GBM-myeloid models. a Percentage of CD45+ cells in patient tissue and
bioprinted PDTs. b Percentage of CD45+, CD14+, and P2RY12+ cells in patient tissues and PDTs. c Percentage of CD31+ cells in patient tissue and
bioprinted PDTs. d Immunofluorescent staining of PDTs derived from 1779-HS and 29259-HS. Green: CD45. Red: GFAP. Scale bars, 100 μm.
e Immunofluorescent staining of Caspase 3 in PDTs derived from 1779-HS and 29259-HS after T-cell treatment. Scale bars, 100 μm. f PCA of drug
responses in GSC, Mg, and Mo under different culture conditions. 2D-Mo: 2D cultured monocyte. G-Mo: GSC co-cultured monocyte. 2D-Mg: 2D
cultured microglia. G-Mg: GSC co-cultured microglia. Sphere: traditional sphere cultured GSC. G-GSC: microglia-coculturedd GSC. O-GSC: monocyte
co-cultured GSC. g Cytotoxicity evaluation of T-cell treatment for GSCs in bioprinted GBM-Mg and GBM-Mo models, with E:T ratios of 1:1.
h Cytotoxicity evaluation of bevacizumab treatment (25 μg/mL) for GSCs in bioprinted GBM-Mg and GBM-Mo models. Immunofluorescence staining
of immunosuppressive markers CD163 and CD206 (i), and tight junction marker ZO-1 and glial cell marker GFAP (j) in GBM-Mg and GBM-Mo models.
Scale bars, 100 μm.

P2RY12+ cells representing cerebral-origin TAMs, and and Bevacizumab. Interestingly, although GSCs in the two
CD31+ cells representing endothelial cells. Accordingly, 3D myeloid-co-culture models showed comparable
our assembly included four multicellular combinations: responses to small molecules, they showed distinct sen-
GBM-Monocyte (GBM-Mo), GBM-Microglia (GBM-Mg), sitivities to these other treatment modalities. This finding
GBM-Monocyte-Endothelial Cell (GBM-Mo-EC), and aligned with our observations in the PDTs, where G-GSCs
GBM-Microglia-Endothelial Cell (GBM-Mg-EC). These were more susceptible to T-cell therapy than O-GSCs
models incorporated tumor cells (CW468), myeloid cells (Fig. 5g). Notably, O-GSCs even demonstrated prolifera-
(THP1 monocytes and HMC3 microglia) representing the tion under T-cell treatment at an effector:target (E:T)
peripheral and cerebral TAM origins, respectively, and ratio of 1:1. Post-treatment analysis of the GBM-Mo
endothelial cells (HUVEC), the critical components of models revealed a reduction in CD8+ cytotoxic T cells
blood vessels. The prepolymer solution used for these and an increase in CD4+ T cells. Within the CD4+ T-cell
models maintained the same concentration and light population, a significantly higher proportion of CD25+
intensity as that used for the PDTs. T cells was observed in GBM-Mo models compared to
These engineered multicellular models enabled the GBM-Mg models, suggesting a more immunosuppressive
assessment of various therapeutic approaches, including microenvironment in the former. Conversely, with Bev-
small molecules identified by GlioML, T-cell therapies, and acizumab treatment, O-GSCs demonstrated a higher
targeted treatments like Bevacizumab. An initial evaluation susceptibility than G-GSCs (Fig. 5h). These results were
of the GlioML’s predictive capacity in these bioprinted consistent with those observed in PDTs, suggesting that
multicellular GBM models revealed a linear correlation the differing myeloid cell compositions in these models
between the GlioML predictions and the observed drug might foster distinct GBM microenvironments. These
responses of GSCs in 3D co-cultured models. GlioML environments could be characterized as either immuno-
demonstrated better predictive performance in GBM-Mg suppressive or angiogenic, both of which are recognized
models compared to the GBM-Mo models (Supplementary as key features of GBM.
Fig. S10a, b). The GlioML predictions indicated that tra-
ditionally cultured GSCs had lower drug resistance, Distinct microenvironment characteristics of bioprinted 3D
represented by a smaller AUC value, than 3D myeloid-co- GBM-myeloid models
cultured GSCs (Supplementary Fig. S10c). Moreover, we To elucidate the factors driving the observed diversity in
identified similar drug response patterns in microglia-co- drug sensitivity, we investigated the bioprinted GBM-
cultured GSCs (G-GSCs) and monocyte-co-cultured GSCs myeloid models and GlioML features. We first profiled
(O-GSCs). However, distinct drug susceptibility patterns the gene expression of GSCs under different culture
emerged in GBM-transformed microglia (G-Mg) and conditions, including traditional sphere culture, bio-
GBM-transformed monocytes (G-Mo) (Fig. 5f). G-Mo printed co-culture (GBM-Mg and GBM-Mo), and bio-
displayed significantly different drug response patterns printed triculture (consisting of GSCs, microglia, and
than co-cultured GSCs and G-Mg. Given that our bio- monocytes together). Traditionally cultured GSCs exhib-
printed co-culture models contained a high proportion of ited significant divergence from all the 3D bioprinted
myeloid cells, mirroring clinical scenarios where TAMs cultures. In contrast, the difference between O-GSC and
constitute up to 30% of GBM tissues, the marginally G-GSC was relatively minor, as indicated by the PCA and
reduced accuracy of GlioML in predicting GSC responses differential gene expression analysis (Supplementary
in the GBM-Mo model could likely be ascribed to the Fig. S11a, b). GSC isolated from the triculture model
GlioML’s analysis focusing primarily on tumor cells. (Tri-GSC) also displayed differences from the co-cultured
We subsequently applied the bioprinted GBM-Mg and O-GSC or G-GSC, indicating that combining the two
GBM-Mo models to assess the efficacy of activated T cells myeloid cell types further impacted the GSCs. We noted
Tang et al. Cell Discovery (2024)10:39 Page 13 of 20

distinct GSC morphologies in the GBM-Mg and GBM- The GlioML feature analysis of cells isolated from the
Mo models, with G-GSCs appearing more fibroblastic and bioprinted models aligned with these observations
O-GSCs exhibiting a more rounded shape (Supplemen- (Fig. 6c). Cytokine signaling and inflammatory features
tary Fig. S11c). Immunofluorescent staining of hypoxia were the most enriched in the G-Mo, while angiogenic
marker CA9 suggested a more hypoxic state in GSCs co- features were dominant in the GBM-Mg microenviron-
cultured with myeloid cells than those bioprinted alone ment. Both G-GSC and G-Mg were markedly enriched in
(Supplementary Fig. S11d). angiogenesis-related pathways. Additionally, hypoxia fea-
Variations in responses to T-cell and Bevacizumab tures, hypoxia-induced p53 activation, and glycolysis were
treatments suggested potential roles of immune interac- most enriched in G-Mg. 3D myeloid-co-cultured GSCs
tions and angiogenesis. Upon further analysis, GBM-Mo showed an increased expression in these pathways com-
models showed elevated immunosuppressive CD163 and pared to their traditionally cultured counterparts. The
CD206 expression levels (Fig. 5i). According to reverse GlioML feature analysis aligned with the outcomes of
transcription quantitative polymerase chain reaction (RT- other analysis pipelines, such as gene set enrichment and
qPCR) data, G-Mo had a substantially higher CD163 gene ontology analyses (Fig. 6d–f; Supplementary Fig.
expression than its 2D counterparts, with this elevated S13a–d). We also found these observations consistent
expression maintained in G-Mo for up to a week. In across cell lines, as similar gene expression changes were
contrast, G-Mg maintained a low level (Supplementary observed in GBM models constructed with another GBM
Fig. S12a). Conversely, G-Mg displayed increased cell line, U251 (data not shown).
expression of angiogenesis-related VEGFA based on RT-
qPCR data (Supplementary Fig. S12b). Immuno- Discussion
fluorescent staining revealed a higher ZO-1 tight junction We present a pioneering integration of 3D bioprinting
protein expression in GBM-Mg models, while it was and machine learning in a clinically relevant context,
undetectable in GBM-Mo models (Fig. 5j). Furthermore, advancing experimental and computational methodolo-
the tube formation assay and endothelial growth evalua- gies to predict and evaluate multimodal tumor treatment
tion demonstrated that the conditioned medium of GBM- responses. Moreover, this integrative approach enabled
Mg models facilitated more meshes’ formation and the exploration of the intricate tumor-immune micro-
increased endothelial cell proliferation (Supplementary environment characteristics.
Fig. S12c, d). The bioprinted PDTs using fresh patient samples and a
Cytokines were pivotal in promoting the immunosup- defined hydrogel composition, as the first of their kind in
pressive or angiogenic states in the two GBM-myeloid glioma research, represented an excellent alternative to
models. In GBM-Mg, the most abundant cytokines were Matrigel-cultured tumor organoids. The cost-effectiveness
IL8, IL6, and CCL2, while in GBM-Mo, IL8, CCL5, and and batch-to-batch uniformity of the bioink also support
G-CSF were most prevalent (Fig. 6a). The CCL2/IL6 axis the scalability of tissue model production. They faithfully
was previously associated with TAM-GBM crosstalk, replicated molecular characteristics, cellular composition,
contributing to enhanced GBM invasiveness11. The co- ECM compositions, and clinical drug responses found in
culture supernatant was compared to the mixed super- patient tumors, validating their utility for drug response
natant of monoculture of involved cell types (Fig. 6b; assessment. Notably, the PDTs demonstrated a capacity to
Supplementary Fig. S12e). Both 3D co-culture models accurately predict TMZ resistance in most recurrent adult
exhibited reduced levels of pro-inflammatory molecules, patients when pMGMT status failed to predict. This
such as interferon-gamma (IFN-γ), interferon-gamma- capacity indicates their potential for enhancing treatment
induced protein 10 (CXCL10), and tumor necrosis planning for GBM patients. Furthermore, these models
factor-alpha (TNF-α) compared to the sum of single offer a flexible platform for evaluating patient responses to
cultures, indicating the induction of an immunosup- various therapies, with the ability to differentiate responses
pressive state through cellular interactions and trans- to T-cell and anti-angiogenesis targeted therapies, such as
formations in both models, with a more pronounced Bevacizumab. Despite these advantages, we recognize the
effect in GBM-Mo (Supplementary Fig. S12f)33. Despite limitations of the bioprinted model in forming mature
the increase in IL8 in both co-culture conditions, other vascular structures within the experimental timeframe, a
cytokines manifested different alterations in GBM-Mg limitation we aim to address in future studies.
and GBM-Mo. Most cytokines either did not change In addition, the GlioML workflow demonstrated profi-
significantly or were found at decreased levels in GBM- ciency in predicting small molecule efficacy in tumor cells
Mg. Conversely, in GBM-Mo, several chemokines, across various in vitro systems, including the bioprinted
including G-CSF, CCL3, CCL4, and CCL5, and cyto- PDTs developed in the study. Recognizing the distinct
kines, such as IL4, IL9, IL10, IL12, and IL15, were mechanistic actions of compounds, we hypothesized that
detected at increased levels. a single algorithm might not generate an optimized
Tang et al. Cell Discovery (2024)10:39 Page 14 of 20

Fig. 6 Distinct microenvironment characteristics in GBM-myeloid models. a Absolute cytokine abundance in GBM-Mg and GBM-Mo models.
b Fold-change comparison of cytokine abundance in co-culture supernatants vs supernatants from monocultures of the corresponding cell types.
c GlioML feature analysis of cells isolated from bioprinted GBM-myeloid models and their traditionally cultured counterparts. d Heatmap
representation of the top differentially expressed genes in G-Mg and G-Mo, with related pathways annotated on the side. e GSEA showing greater
enrichment of the angiogenesis pathway in G-Mg. f GSEA demonstrating enhanced chemokine production pathway in G-Mo.
Tang et al. Cell Discovery (2024)10:39 Page 15 of 20

predictor for all drugs. The ensemble model of GlioML Flow cytometry analysis
leveraged the stacking strategy, merging multiple cate- Isolated patient cells were resuspended in a DPBS buffer
gories of algorithms and validated our hypothesis by with 40 μg/mL DNase I and analyzed using CytoFLEX
surpassing the performance of all individual algorithms. flow cytometer (Beckman Coulter). For staining, 1 μL of
This approach underscored the necessity of context- conjugated CD45 antibody (Biolegend, Cat# 304012) was
specific algorithm selection for optimal predictions and added to 100 μL cell suspension (0.1 million cells) and
bypassed the need for arbitrary model selection and incubated for 20 min on ice in the dark. After washing 3
manual parameter tuning typically required in traditional times with DPBS (Gibco), cells were stained with 7-AAD
methodologies. The synergy between GlioML and bio- (Biolegend, Cat# 420403) for 10 min at room temperature
printing yielded a powerful tool for drug assessment and in the dark. A total of 10,000 events were collected for
illuminated differential microenvironments within tumor- each sample. Data were analyzed using FlowJo v10.8. All
immune models. This integration identified a cluster of experiments were performed in triplicates.
promising yet unexpected candidate compounds. This
integrated strategy identified a cohort of unexpected yet GlioML workflow
promising candidate compounds. RSL3, a ferroptosis The GlioML input preparation involved generating
activator, demonstrated strong efficacy on GSCs and bio- features and labels, followed by data cleaning and nor-
printed models. Lovastatin, a drug typically used for malization. Cancer cell line gene expression and ther-
managing high cholesterol and triglyceride levels, displayed apeutic response AUC data were accessed through the
variable efficacy across different PDTs, yet displayed good CCLE23 and the CTRP24,34,35. AUCs were used as labels.
effectiveness in a few samples. GlioML in combination Gene-expression TPM values for clinical samples and cell
with 3D bioprinted multicellular models also revealed line samples generated in this study were calculated from
distinctive drug susceptibility and gene expression patterns raw RNA-seq data and transformed to log2(TPM + 1) to
in GSCs and myeloid cells. These findings, particularly match the CCLE data format. The ssGSEA module on
concerning GBM-transformed microglia and monocytes, GenePattern36 was used to calculate the ssGSEA scores of
enhance our understanding of their unique roles in each sample using a specific list of gene sets curated from
immunosuppression and angiogenesis. In summary, the the Molecular Signatures Database (v7.5)25,37. The cal-
versatility of this integrated computational and experi- culated ssGSEA scores for all curated gene sets were used
mental approach supports its applicability across a range of as the feature values. The data cleaning process involved
cancer types, indicating its promising potential for future filtering out samples that did not have labels. Normal-
advancements in personalized cancer therapeutics. ization was performed by scaling all signature scores to
the range of [0, 1].
Materials and methods AutoGluon 0.5.1 was used to build the GlioML
Patient tissue collection workflow. Nine basic models, including FastAI, Neural
Glioma tissue specimens or blood samples were col- Network Torch, LightGBM, LightGBM_XT, Light-
lected from Huashan Hospital. The research undertaken GBM_Large, Random Forest, Extra Trees, KNN_uni-
was approved by the research ethics committee at Hua- form, KNN_distance, and a weighted ensemble model
shan Hospital, Fudan University, under the ethics were included. AutoGluon Tabular Predictor was used
approval numbers: KY2021-670, KY2021-059, KY2023- with specific parameters, including ‘best_quality’ presets,
846. All participating patients provided written consent bagging and multi-layer stack ensembling, 3 repeats of
for using their samples in research endeavors. Pathologists 5-fold cross-validation, and random search hyperpara-
confirmed the classification of the tumor samples. meter tuning for 512 times for each model. The problem
Detailed patient demographic and clinical data are avail- type was regression; the evaluation metric was mean
able in Table 1 and Supplementary Fig. S1a. squared error.
After initial training, feature importance scores were
Isolation of cells from patient specimens generated using feature_importance(), with sub-
Surgical samples were cut or minced into 1–3 mm sample_size set to 5000 and num_shuffle_sets set to 10 for
pieces and rinsed thoroughly with Dulbecco’s Phosphate- highly accurate importance and P-value estimates. The
Buffered Saline (DPBS, Gibco). The mixture was then scores were presented as a ranking of the most important
treated with collagenase type I (Yeasen) and rocked on a features to the least important features. These scores were
shaker at 37 °C for 1 h to aid digestion. The digested used to interpret the model and identify the most pre-
samples were rinsed with DPBS. The resulting cell sus- dictive features. Then the top 20, top 10, or top 3 percent
pensions were collected by centrifugation at 200× g for of the features were extracted, and the models were re-
5 min. A 5-min red cell lysis procedure was performed if trained based on the newly formulated training data using
excessive red blood cells were present. the same parameters as the initial training.
Tang et al. Cell Discovery (2024)10:39 Page 16 of 20

Cell culture Bioprinting PDTs and GBM-myeloid models

CW468 cells were established by Dr. Jeremy N. Rich lab PDTs were bioprinted using a light-based bioprinter,
(UCSD) and cultured in Neurobasal medium (NBM, Biocube (Cyberiad Biotechnology), with a wavelength of
Gibco, Cat# 21103049) supplemented with 1% B27 minus 405 nm, and multicellular GBM-myeloid models were
vitamin A, 1% L-glutamine, 1% sodium pyruvate, 1% created with a digital light processing bioprinting system,
penicillin/streptomycin (P/S), 10 ng/mL basic human also with a wavelength of 405 nm. Both systems employed
fibroblast growth factor (Peprotech, Cat# AF-100-18B), matching printing parameters, including a power density
and 10 ng/mL human epidermal growth factor (Pepro- (irradiance) of 50 mW/cm2 and a scaffold thickness of
tech, Cat# AF-100-15). THP1 cells (American Type Cul- 0.5 mm. For the PDTs, a cell suspension containing 30
ture Collection, ATCC, Cat# TIB-202) were cultured in million cells/mL was prepared, with 2 μL of the cell-
Roswell Park Memorial Institute (RPMI) 1640 medium material suspension required for each sample. For GBM-
supplemented with 10% FBS and 1% P/S. HMC3 cells myeloid cells, a cell suspension composed of equal parts
(ATCC, Cat# CRL-3304) were cultured in Minimum 40 million/mL GBM cells and 20 million/mL THP1 or
Essential Medium (MEM, Gibco) medium with 10% fetal HMC3 cells was prepared, using 5 μL of cell-material
bovine serum (FBS), 1% P/S, and 1% non-essential amino solution for each printed sample.
acids. HUVEC (Cell Applications, Cat# 200p-05n) were The prepolymer solution used for all bioprinting processes
cultured in endothelial cell growth medium v1 (Cell consisted of 8% GelMA, 2% HAMA (Yuju Technology), and
Applications, Cat# 211-500) supplemented with 1% P/S. 0.2% lithium phenyl-2,4,6-trimethylbenzoylphosphinate
U251 cells were cultured in DMEM medium with 10% (TCI Chemicals). This solution was thoroughly mixed and
FBS and 1% P/S. All cells were maintained at 37 °C with stored at 37 °C in the dark before use. The cell suspension
5% CO2 and passaged with TrypLE (Gibco) or Accutase was combined with the prepolymer solution at a 1:1 ratio
(Stemcell Technologies) every two to three days. Cell line immediately before printing to ensure maximum cell via-
authentication was conducted for immortalized cell lines, bility. The cell-material mixture was exposed to the bio-
including HMC3 and THP1, using short tandem repeat printer for an optimized 15–20 s duration. Afterward, the
(STR) profiling and comparison to known reference bioprinted constructs were rinsed with DPBS and cultured
profiles to ensure data integrity and reliability. in a maintenance medium at 37 °C and 5% CO2.
The maintenance medium for PDTs consisted of DMEM/
Drug response evaluation F12 supplemented with 1% B27, 1% P/S, 10 ng/mL basic
Drug compounds were pre-dissolved in dimethyl sulf- human fibroblast growth factor (Novoprotein), and
oxide (DMSO) or de-ionized water to a stock solution 10 ng/mL human epidermal growth factor (Novoprotein).
based on solubility and stored at –20 °C. Cisplatin was The maintenance medium for GBM-myeloid cells con-
obtained from Qilu Pharmaceutical. Lobaplatin was tained equal parts CW468 medium and THP1/HMC3
obtained from Hainan Changan International Pharma- medium. Bioprinted constructs were cultured 7 days before
ceutical. All other tested compounds were obtained from analytical experiments and 3 days before drug treatments.
Selleck, including RSL (Cat# S8155), Erlotinib (Cat# Model monitoring was conducted regularly using optical
S1023), TMZ (Cat# S1237), JW-55 (Cat# S6745), TGX- microscopy.
221 (Cat# 1169), Dasatinib (Cat# S1021), Lovastatin (Cat#
S2061), Daporinad (Cat# S2799), and Docetaxel (Cat# Cell isolation from bioprinted samples
S1148). The bioprinted samples were incubated with 1 mg/mL
Before drug treatment, the compounds were diluted in collagenase in HBSS (Gibco) for 30–60 min at 37 °C with
NBM to a concentration that ensured DMSO was less constant shaking to digest the hydrogel. The samples were
than 0.1% when treating cells. GSCs were bioprinted at then filtered through a 70-μm cell strainer to remove any
5000 cells per well with 50 μL of medium in 96-well plates remaining clumps. The resulting single-cell suspension
and allowed to culture overnight at 37 °C with 5% CO2. was collected and centrifuged at 200× g for 5 min to
The next day, 50 μL of 2× drug solution was added to the collect cells for subsequent analysis.
wells. Cell viability was assessed at the 72-h time point
after drug treatment using CellTiter-Glo 3D assay (Pro- WES analysis
mega). Briefly, 100 μL of CellTiter-Glo 3D was added, and Patient samples and their corresponding PDTs and
the plates were shaken at 300 rpm for 10 min to ensure PDOs were subjected to WES analysis and germline
uniform mixing. The luminescence was measured on a corrected using matched blood samples. DNA was
plate reader (Tecan) using luciferase reading. The IC50 extracted using the QIAamp FFPE DNA Tissue Kit
values were calculated from dose-response curves gener- (Qiagen). DNA degradation and contamination were
ated using GraphPad Prism 9 software. All experiments assessed on a 1% agarose gel, while concentrations were
were performed as triplicates. determined using the Qubit DNA Assay Kit (Life
Tang et al. Cell Discovery (2024)10:39 Page 17 of 20

Technologies) and the Qubit 2.0 Fluorometer (Life Gene expression analysis
Technologies). DNA was further fragmented to Total RNA was extracted from cell pellets using Trizol and
150–230 bp, and adapters were ligated at both ends of the Microprep kit (Zymo Research). For RNA-seq, paired-end
fragments. The extracted DNA was then PCR-amplified FASTQ sequencing reads were generated and trimmed using
and hybridized to the KAPA HyperExome Probes (Roche) Trim Galore version 0.6.2 and Cutadapt version 2.3,
for enrichment. Libraries were subsequently loaded onto respectively. Transcript quantification was performed using
the Illumina NovaSeq 6000 platform (Honsunbio). Patient Salmon version 0.13.1 in the quasi-mapping mode using
samples, PDTs, and blood DNA samples were respectively transcripts derived from human Gencode release 30
sequenced to average depths of > 400×, > 400×, and > (GRCh38). Salmon quant files were converted using the R
100× in targeted exonic regions. package Tximport, and differential expression analysis was
The obtained reads were processed using the following performed using DESeq2. Genes with an adjusted P-value
pipeline. We aligned the sequencing reads to the human (Padj) less than 0.05 were considered differentially expressed.
reference genome (hg19) using the Burrows-Wheeler To perform gene set enrichment analysis (GSEA), we utilized
Aligner v0.7.17 and then processed them with Picard the GSEA desktop application version 4.1.0 provided by the
v1.119 to remove duplicates. Genome Analysis Toolkit Broad Institute. The gene sets used for enrichment analysis
v4.1.9 was used to conduct INDEL realignment and were obtained from the Molecular Signatures Database
recalibrate base quality scores. Somatic SNVs and INDELs (MSigDB). The enrichment analysis was performed using the
were identified using Mutect2 and Varscan_Indel v2.4.2, gene set permutations option, and gene sets with a false
respectively, and annotated with ANNOVAR. Copy discovery rate (FDR) less than 0.25 were considered sig-
number variations (CNVs) and gene fusions were detected nificantly enriched. The cutoff for gene ontology (GO) ana-
using CNVkit v0.9.9 and Lumpy v0.2.13, respectively. lysis was Padj < 0.01 and log2 fold change > 3. To generate a
Purity and ploidy for each sample were manually calcu- heatmap, the gene expression data is clustered using a
lated by comparing the variant allele frequency (VAF) of hierarchical clustering approach with the “ward.D2” method
somatic SNV and CNV log2 values of multiple variants and “euclidean” distance, both for rows and columns. Pear-
within each sample. The final average ploidy/purity was son’s correlation coefficient I was computed to measure the
taken from the Sequenza or ABSOLUTE estimate that linear correlation between each Tissue–PDT pair. Genes that
most closely matched the manual calculation. No precise were not detected across a pair of samples (i.e., having zero
purity estimates were made for tumor samples with values) were excluded from each comparison.
extremely low purity (< 15%). For RT-qPCR, cDNA was synthesized from total RNA
SNV concordance between tumor–PDT pairs was using the First-Strand cDNA Synthesis Kit (New England
determined by examining the overlap of variant calls Biolabs). After obtaining the cDNA, PowerUp SYBR Green
and variant allelic fractions. For each SNV identified in Master Mix (Thermo Fisher Scientific) was used to carry
the tumor or PDT, SAMTOOLS Pileup (with a mini- out RT-qPCR analysis with specific primers for VEGFA
mum base quality and minimum mapping quality of (NM_001025366.3) and CD163 (NM_203416.3). VEGFA
10) was run at this position for both samples to com- Forward: TTGCCTTGCTGCTCTACCTCCA. VEGFA
pute the variant allele fractions. If read evidence for the Reverse: GATGGCAGTAGCTGCGCTGATA. CD163
SNV was present in both samples (and therefore Forward: AAAAAGCCACAACAGGTCGC. CD163
VAF > 0), the SNV was considered concordant. SNVs Reverse: CTTGAGGAAACTGCAAGCCG. For each
called by two or more variant callers in at least one of sample, the RT-qPCR reactions were conducted in tripli-
the samples were included. CNVs were compared cate, and the expression levels of the target genes were
between PDTs and tumors by plotting CNV log2 values normalized to GAPDH (NM_002046.7) as a reference
across chromosomes. A threshold of –0.235 and 0.2 gene. GAPDH Forward: ACAACTTTGGTATCGTG-
was used to delineate the cutoff for deletions and GAAGG. GAPDH Reverse: GCCATCACGCCA-
amplifications, respectively, based on a diploid sample CAGTTTC. The data were analyzed to determine the fold
with 30% purity. Neutral segments were colored in change in G-Mg and G-Mo gene expression at day 7 and
green, and deletions/amplifications in red. The y-axis day 3 compared to their 2D control. We used a two-tailed
range was smaller for tumor samples to facilitate CNV t-test for statistical analysis, and P-values less than 0.05
identification. SVs were compared between PDTs and were considered statistically significant.
tumors by plotting intra- and inter-chromosome
rearrangements based on LUMPY results. The SNV Tube formation assay and endothelial cell growth
landscape was visualized using the Bioconductor A composite medium was prepared using equal parts
package GenVisR v1.8.0. The CNV per gene heatmap Neurobasal medium, MEM, and RPMI 1640. For tube
was generated using the Bioconductor package Com- formation assessment, conditioned medium from GBM-
plexHeatmap v1.17.1. Mg and GBM-Mo was collected after 72 h of incubation.
Tang et al. Cell Discovery (2024)10:39 Page 18 of 20

In each well of an angiogenesis slide (Ibidi), 10 μL of at a concentration of 25 μg/mL. The bioprinted models
growth factor-reduced Matrigel was added and incubated containing both HUVEC and the cell mixture, or the
at 37 °C for 30 min. Subsequently, 50 μL of HUVEC sus- bioprinted PDTs, were treated with Bevacizumab for
pension (10,000 cells) was introduced into each well. Cells three days. The control group contained the bioprinted
were cultured in conditioned media from GBM-Mg or model without bevacizumab treatment. The assay was
GBM-Mo, with non-supplemented endothelial cell performed in triplicate for all conditions.
growth medium v1 serving as a negative control and The impact of T cells or Bevacizumab on the PDTs was
endothelial cell growth medium containing 50 ng/mL assessed through immunofluorescent staining for Caspase 3,
recombinant VEGF as a positive control. Mesh counts and flow cytometry to quantify CD31+ cell populations. The
were determined using the ImageJ plugin. cytotoxicity of T cells or bevacizumab to GSCs in the
To evaluate endothelial cell growth, 5000 HUVECs were multicellular models was measured using a luciferase assay
seeded and cultured overnight in endothelial cell growth system (Promega). The bioprinted models were digested and
medium v1. Simultaneously, GBM-Mg and GBM-Mo fully lysed, and the substrate was mixed with cell lysates.
models were printed. The following day, conditioned Luminescence was measured using a luminometer (Tecan).
medium from GBM-Mg and GBM-Mo was collected and The assay was performed in triplicate for all conditions.
added to the HUVEC culture, with non-supplemented
endothelial cell growth medium v1 used as a control. Immunofluorescent staining
HUVEC medium was replaced every other day using For immunofluorescent staining, the models were fixed
collected conditioned medium. HUVEC growth was with 1 mL of 4% paraformaldehyde in PBS overnight at
assessed on days 1, 3, and 5 post-incubation in the con- 4 °C. Following PBS washing, the cells were permeabilized
ditioned medium using CellTiter-Glo (Promega), with the with 0.1% Triton X-100 for 20 min and subsequently
experiment conducted in triplicate. blocked with 5% BSA for 90 min. Primary antibodies
against CD163 (1:100 dilution, Proteintech, Cat# 16646-1-
Isolation, culture, and activation of primary human T cells AP), CD206 (1:100 dilution, Abcam, Cat# 91992 S), ZO-1
Human peripheral blood mononuclear cells (PBMCs) (1:100 dilution, Invitrogen, Cat# 339100), Caspase 3 (1:100
were isolated from diluted whole blood (1:1 in DPBS) (San dilution, Proteintech, Cat# 66470-2-lg), CD45 (1:100
Diego Blood Bank) using density gradient centrifugation dilution, Proteintech, Cat# 65109-1-lg), or GFAP (1:100
with lymphocyte separation medium (Corning) or pur- dilution, Sigma, Cat# G3893) were diluted in staining
chased from Hycells Biotechnology. Primary human T cells buffer (BioLegend, Cat# 420201) and incubated at 4 °C
were then isolated from PBMCs using a Pan T-Cell Iso- overnight. Samples were rinsed again with PBS for 3 times
lation Kit (Miltenyi) following the manufacturer’s instruc- and then incubated with secondary antibodies (1:200)
tions using magnetic-activated cell sorting (MACS) conjugated with fluorophores for 2 h at room temperature.
columns. Isolated primary human T cells were cultured DAPI (1:1000) was used for nuclear staining. Samples were
with RPMI 1640 with 10% FBS, 1% P/S, and 100 U/mL imaged using a confocal microscope (Leica SP8), and
recombinant human IL-2 (PeproTech). Primary human image analysis was performed using ImageJ software. A
T cells were activated with Dynabeads Human T-Expander minimum of three samples were analyzed per condition.
CD3/CD28 (Gibco) for 72 h before usage. Cells were cul-
tured at 37 °C with 5% CO2 in a humidified incubator. Luminex for cytokine analysis
Culture supernatants were collected from 3D bioprinted
Activated T cell or bevacizumab cytotoxicity assay models on day 7. The supernatants were centrifuged at
Bioprinted PDTs or GBM-myeloid models were fabri- 3000 rpm for 5 min to remove cellular debris. The assay
cated as described earlier. For the multicellular models, was performed using the Human 27 Cytokine/Chemokine
luciferase-labeled GSCs were used. For the activated Panel. The samples were read using a Luminex 200
T-cell assay, the bioprinted PDTs or GBM-myeloid instrument (Luminex Corporation), and the results were
models were incubated with activated T cells at an analyzed using Milliplex Analyst 5.1 software. Replicates
effector-to-target (E:T) ratio of 1:1. The T cells were were measured for each condition.
added to the bioprinted PDTs or models and incubated at
37 °C with 5% CO2 for 3 days. The control group con- Statistical analysis
tained the bioprinted model without T cells. The assay Statistical analysis was performed using GraphPad Prism
was performed in triplicate for all conditions. 9. The data presented are the mean ± SD from a minimum
For bevacizumab evaluation, HUVEC was mixed with of triplicates. We used a two-tailed t-test to analyze the
the cell mixture at 20 million/mL for printing of the differences between two groups, and for multiple group
multicellular models. The other printing parameters were comparisons, we employed one-way ANOVA followed by
the same as previously described. Bevacizumab was tested Tukey’s post hoc test. We considered a P-value below 0.05
Tang et al. Cell Discovery (2024)10:39 Page 19 of 20

to be statistically significant. To assess the relationships Publisher’s note

between GlioML predictions and measured drug respon- Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.
ses, we fit and calculated the coefficient of determination
(R2) to determine the amount of variance explained by the Supplementary information The online version contains supplementary
models. To perform dimensionality reduction and visua- material available at https://doi.org/10.1038/s41421-024-00650-7.
lize complex datasets, we utilized the scikit-learn library
Received: 27 June 2023 Accepted: 18 January 2024
v1.1.0 in Python to conduct PCA and t-distributed sto-
chastic neighbor embedding (t-SNE).

Acknowledgements
This work was supported by the National Key R&D Program of China References
(2022YFC3401600 to Y.Y.), the National Natural Science Foundation of China 1. Ostrom, Q. T. et al. CBTRUS statistical report: primary brain and other central
(91959112 to Y.Y. and 32301149 to M.T.), the Science and Technology Commission nervous system tumors diagnosed in the United States in 2015–2019. Neuro
of Shanghai Municipality (21140900700 to M.T.), and the Clinical Research Plan of Oncol. 24, v1–v95 (2022).
SHDC (SHDC2020CR2018B to Y.Y.). The authors thank Dr. Chao Tang from 2. Hegi, M. E. et al. MGMT gene silencing and benefit from temozolomide in
Huashan Hospital for assistance in clinical sample collection, Yanbo Wang and glioblastoma. N. Engl. J. Med. 352, 997–1003 (2005).
Wenyang Deng from Cyberiad Biotechnology for assistance in PDT generation 3. Touat, M. et al. Mechanisms and therapeutic implications of hypermutation in
and drug testing, and Dr. Shuo Mu from Honsunbio for WES analysis support. gliomas. Nature 580, 517–523 (2020).
4. Omuro, A. et al. Radiotherapy combined with nivolumab or temozolomide
for newly diagnosed glioblastoma with unmethylated MGMT promoter:
Author details
1 An international randomized phase III trial. Neuro Oncol. 25, 123–134
Shanghai University of Traditional Chinese Medicine, Shanghai, China.
2 (2023).
Department of NanoEngineering, University of California San Diego, La Jolla,
5. Cloughesy, T. F. et al. Neoadjuvant anti-PD-1 immunotherapy promotes a
CA, USA. 3Department of Statistics, University of California Davis, Davis, CA,
survival benefit with intratumoral and systemic immune responses in recur-
USA. 4Department of Neurosurgery, Huashan Hospital, Shanghai Medical
rent glioblastoma. Nat. Med. 25, 477–486 (2019).
College, Fudan University, Shanghai, China. 5National Center for Neurological
6. Reardon, D. A. et al. Effect of nivolumab vs bevacizumab in patients with
Disorders, Shanghai, China. 6Shanghai Key Laboratory of Brain Function and
recurrent glioblastoma: the CheckMate 143 phase 3 randomized clinical trial.
Restoration and Neural Regeneration, Shanghai, China. 7Immunology
JAMA Oncol. 6, 1003–1010 (2020).
Laboratory, Neurosurgical Institute of Fudan University, Shanghai, China.
8 7. Brennan, C. W. et al. The somatic genomic landscape of glioblastoma. Cell 155,
Shanghai Clinical Medical Center of Neurosurgery, Shanghai, China. 9Cyberiad
462–477 (2013).
Biotechnology Ltd., Shanghai, China. 10Department of Human Biology,
8. Sammut, S.-J. et al. Multi-omic machine learning predictor of breast cancer
University of California San Diego, La Jolla, CA, USA. 11Department of
therapy response. Nature 601, 623–629 (2022).
Bioengineering, University of California Berkeley, Berkeley, CA, USA.
12 9. Varn, F. S. et al. Glioma progression is shaped by genetic evolution and
Department of Bioengineering, University of California San Diego, La Jolla,
microenvironment interactions. Cell 185, 2184–2199.e16 (2022).
CA, USA. 13Alfred E. Mann Department of Biomedical Engineering, University of
10. Wick, W. et al. Lomustine and bevacizumab in progressive glioblastoma. N.
Southern California, Los Angeles, CA, USA. 14Department of Orthopaedic
Engl. J. Med. 377, 1954–1963 (2017).
Surgery, University of California San Diego, La Jolla, CA, USA
11. Hambardzumyan, D., Gutmann, D. H. & Kettenmann, H. The role of microglia
and macrophages in glioma maintenance and progression. Nat. Neurosci. 19,
Author contributions 20–27 (2016).
M.T., S.J., and Y.Y. conceptualized the study. M.T. and S.J. led the methodology 12. Pombo Antunes, A. R. et al. Single-cell profiling of myeloid cells in glio-
development. M.T., S.J., X.H., C.J., Y.G., Y. Qi, Y.X., E.Y., N.Z., E.B., D.Y., Y. Qu., L.L., blastoma across species and disease stage reveals macrophage competition
and D.B. conducted the experiments. M.T. and S.J. handled the visualization of and specialization. Nat. Neurosci. 24, 595–610 (2021).
data. Funding was acquired by M.T. and Y.Y. Project administration was carried 13. Sa, J. K. et al. Transcriptional regulatory networks of tumor-associated mac-
out by X.H., C.J., and Y.Y. Supervision was provided by M.T. and Y.Y. The original rophages that drive malignancy in mesenchymal glioblastoma. Genome Biol.
draft was written by M.T., and the review and editing were performed by M.T., 21, 216 (2020).
S.J., and Y.Y. All authors read and approved the final manuscript, contributed 14. Yu, C. et al. Photopolymerizable biomaterials and light-based 3D printing
feedback throughout the development process, and agreed to the final strategies for biomedical applications. Chem. Rev. 120, 10695–10743
submitted version. (2020).
15. Akbari, M. & Khademhosseini, A. Tissue bioprinting for biology and medicine.
Data availability Cell 185, 2644–2648 (2022).
All raw sequencing data has been deposited in the Sequence Read Archive 16. Tang, M. et al. Three-dimensional bioprinted glioblastoma microenvironments
(BioProject accession number: PRJNA952850) and National Omics Data model cellular dependencies and immune interactions. Cell Res. 30, 833–853
Encyclopedia (project number: OEP004093). All other data can be obtained (2020).
upon request by contacting the corresponding authors. All biological materials 17. Kong, J. et al. Network-based machine learning in colorectal and bladder
employed in this manuscript will be provided upon request. organoid models predicts anti-cancer drug efficacy in patients. Nat. Commun.
11, 5485 (2020).
Material availability 18. Geeleher, P., Cox, N. J. & Huang, R. S. Clinical drug response can be predicted
All biological materials employed in this manuscript will be provided upon using baseline gene expression levels and in vitro drug sensitivity in cell lines.
request. Genome Biol. 15, R47 (2014).
19. Vanguri, R. S. et al. Multimodal integration of radiology, pathology and
genomics for prediction of response to PD-(L)1 blockade in patients with non-
Code availability small cell lung cancer. Nat. Cancer 3, 1151–1164 (2022).
All computational algorithms used in the manuscript are cited in the 20. Boehm, K. M. et al. Multimodal data integration using machine learning
respective figure legends and detailed in the Materials and methods section. improves risk stratification of high-grade serous ovarian cancer. Nat. Cancer 3,
Further information can be provided upon request. 723–733 (2022).
21. McBain, C. et al. Treatment options for progression or recurrence of glio-
Conflict of interest blastoma: a network meta‐analysis. Cochrane Database Syst. Rev. 2021,
The authors declare no competing interests. CD013579 (2021).
Tang et al. Cell Discovery (2024)10:39 Page 20 of 20

22. Esteller, M. et al. Inactivation of the DNA-repair gene MGMT and the clinical 30. Geurts, P., Ernst, D. & Wehenkel, L. Extremely randomized trees. Mach. Learn.
response of gliomas to alkylating agents. N. Engl. J. Med. 343, 1350–1354 63, 3–42 (2006).
(2000). 31. Altman, N. S. An introduction to kernel and nearest-neighbor nonparametric
23. Ghandi, M. et al. Next-generation characterization of the cancer cell line regression. Am. Statistician 46, 175–185 (1992).
encyclopedia. Nature 569, 503–508 (2019). 32. Caruana, R., Niculescu-Mizil, A., Crew, G. & Ksikes, A. Ensemble selection from
24. Rees, M. G. et al. Correlating chemical sensitivity and basal gene libraries of models. ICML '04: Proceedings of the twenty-first international
expression reveals mechanism of action. Nat. Chem. Biol. 12, 109–116 conference on Machine learning. https://doi.org/10.1145/1015330.1015432
(2016). (2004).
25. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based 33. Castro, F., Cardoso, A. P., Gonçalves, R. M., Serre, K. & Oliveira, M. J. Interferon-
approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. gamma at the crossroads of tumor immune surveillance or evasion. Front.
Sci. USA 102, 15545–15550 (2005). Immunol. 9, 847 (2018).
26. Liberzon, A. et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics 34. Basu, A. et al. An interactive resource to identify cancer genetic and
27, 1739–1740 (2011). lineage dependencies targeted by small molecules. Cell 154, 1151–1161
27. Howard, J. & Gugger, S. Fastai: a layered API for deep learning. Information 11, (2013).
108 (2020). 35. Seashore-Ludlow, B. et al. Harnessing connectivity in a large-scale small-
28. Ke, G. et al. LightGBM: a highly efficient gradient boosting decision tree. molecule sensitivity dataset. Cancer Discov. 5, 1210–1223 (2015).
NIPS’17: Proceedings of the 31st International Conference on Neural Infor- 36. Reich, M. et al. GenePattern 2.0. Nat. Genet. 38, 500–501 (2006).
mation Processing Systems. 3149–3157 (2017). 37. Subramanian, A. et al. A next generation connectivity map: l1000 platform and
29. Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001). the first 1,000,000 profiles. Cell 171, 1437–1452.e17 (2017).