Study participants and samples

Forty adult CL patients participated in this study, and were recruited in our outpatient clinics, one in the city of Tumaco (South Pacific coast of Colombia), and another in the urban center of Cali. One patient was excluded due to progression to mucosal leishmaniasis, in accordance with the clinical protocol designed for this study (Fig. 1A). Of the remaining participants, 35 completed clinical follow-up. One participant did not donate blood or tissue samples and was therefore not included in the transcriptomic analyses. Most of the 34 enrolled participants were males of Afrocolombian descent, presented with ulcerated skin lesions, and most were infected with L.V. panamensis, the remaining with L. V. braziliensis. Five patients were treated with miltefosine (MLF), all cured. Twenty-nine patients were treated with GLUC, 19 cured, and 10 experienced TF. Due to the small sample size of MLF-treated patients and the absence of TF in this group, only samples from GLUC-treated patients (29 in total) were carried forward for transcriptomic analyses (Fig. 1B). Adherence to treatment in GLUC-treated patients was > 95% (measured as the proportion of ampoules administered from those prescribed). Significant differences were noted in the age and ethnicity of patients who cured vs. patients with TF in the aggregated data of patients recruited in Cali and Tumaco (Supplementary Table 1a). These two variables (age and ethnicity) were collinear (Supplementary Fig. 1A, B), and the magnitude of their effect on the outcome of treatment was marginal (Supplementary Fig. 1C, D). No further evidence of the influence of any of these variables on the outcome of treatment was identified in a multivariable logistic regression analysis (Supplementary Fig. 1E for the full cohort, or Supplementary Fig. 1F for patients recruited in Tumaco only).

Flow chart (A) showing the recruitment of participants in the study ending with a total of 34 enrolled participants. Study design (B) detailing participant enrollment, treatment, sample collection and analysis, and determination of treatment outcome. Created in BioRender. Gomez, M. (2025) https://BioRender.com/g83q151.

Leishmania species, drug resistance, and virulence of the infecting strain may contribute to TF. In L. V. braziliensis and L. V. guyanenesis infections, virulence has been associated with the presence of the Leishmania RNA virus (LRV)^16,17. To assess the possible impact of these variables, we evaluated the Leishmania isolates from study participants for susceptibility to pentavalent antimony and the presence of LRV. Leishmania strains were isolated from 24 patients treated with GLUC (14 cures and 10 TF); isolates from the remaining 5 could not be propagated. The predominant species isolated was L.V. panamensis (n = 19), followed by L.V. braziliensis (n = 4) and L.V. guyanensis (n = 1). No statistical difference was observed between the Leishmania species and the outcome of treatment in this patient cohort (Supplementary Table 1). As intracellular amastigotes, 10 strains were susceptible to GLUC, 12 were resistant (of which 5 were isolated from patients who were cured, and 7 from TF), and data was unavailable for 2 (Supplementary Fig. 2A). No statistically significant difference was found (Fisher´s P = 0.23). All strains evaluated were negative for LRV.

To explore the possible contribution of genomic/transcriptomic differences of the parasite to the outcome of treatment, Leishmania transcriptomes were analyzed from the dual RNA-Seq data collected from lesion biopsies, Mo, Nφ, and Eφ samples from CL patients before and over the course of treatment. Only samples that passed a threshold of single mapped parasite reads ≥ 30000, and ≥ 3000 observed genes were used for the sufficient representation of the parasite genome, resulting in a total of 15 samples analyzed (corresponding to 8 biopsy samples, 1 Eφ, 1 Mo, and 5 Nφ samples from 11 patients; 6 cures and 5 with TF) (Supplementary Data 1). No significant difference in the parasite transcriptome profiles could be discerned by PCA using either uncorrected (Supplementary Fig. 2B) or surrogate variable analysis (SVA) adjusted data (Supplementary Fig. 2C). Collectively, LRV content, drug susceptibility testing, Leishmania species typing, and transcriptomic analyzes from infected human samples rule out a significant contribution of the parasite to the outcome of treatment in this cohort of patients.

Global assessment of samples and human transcriptomic profiles

Lesion biopsy samples were collected before initiation of treatment, and peripheral blood samples obtained pre-treatment (Pre-Tx), mid-way through treatment (Mid-Tx, day 8), and at the end of treatment (End-Tx, day 20). Mo, Nφ, and Eφ were isolated from all blood samples. The cell populations were evaluated by flow cytometry and light microscopy, and purity confirmed > 95% (Supplementary Fig. 3). cDNA libraries were constructed and sequenced from a total of 186 samples (Supplementary Data 1 and Fig. 2A). Following low coverage filtering (Supplementary Fig. 4), two samples were removed. We used principal component analysis (PCA) and a correlation heatmap to visualize the relationship between samples (Fig. 2B, C). The PCA resulted in the expected grouping by cell type (Fig. 2B). A similar clustering was observed in the hierarchical clustering analysis (Fig. 2C). These data support the quality, reproducibility, and specificity of transcriptomes derived from the isolated blood cells from our patient cohort.

RNA-seq was carried out on biopsies collected pre-treatment and on 3 types of leukocytes (Mo, Nφ, and Eφ) collected pre-treatment (Pre-Tx), mid-treatment (Mid-Tx), and at the end of treatment (End-Tx). C = cures, F = treatment failure (A). A principal component analysis (PCA) plot (B) and heatmap of hierarchical clustering analysis using pairwise correlations (C) are shown for all samples. The analyses were performed using all annotated protein-coding genes following filtering for low counts, CPM, quantile normalization, and log₂ transformation. In the PCA plot, the first two principal components are shown on the X and Y-axes, respectively, with the proportion of total variance attributable to that PC indicated. Each sample is represented as a single point with color indicating cell type and shape indicating treatment outcome. Colors along the top of the heatmap indicate the cell type. Colors within the inset of the heatmap represent correlation values. The green line represents the frequency of correlation values.

The use of two clinics for patient recruitment (one in Tumaco and another in Cali) necessitated the evaluation of the data to account for possible clinic-associated batch effects. When all samples were colored by clinic on the same PCA plot, no grouping by clinic was discernible within cell types (Supplementary Fig. 5A – all samples); however, a grouping of samples by cell type revealed a significant amount of variance which separated Tumaco from Cali, confirming the presence of a batch effect attributable to the clinic (Supplementary Fig. 5B–D). Similar analyses were performed evaluating the effect of ethnicity (which was also cross-correlated with age, without finding any substantial contribution of this variable to the overall variability of the transcriptomic data (Supplementary Fig. 5F–H). To highlight the variability between clinics, we modeled the ‘clinic’ variable in SVA and used the surrogate variables-modified counts to generate the PCA plots, further supporting the hypothesis of a strong batch effect associated with the clinic (Fig. 3A–D).

Principal Component Analysis (PCA) plots of the global gene expression profiles of all samples combined (A) and individual blood cells (eosinophils, monocytes, and neutrophils) and lesion biopsies collected from study participants (B–E) at two different clinics (one in Tumaco and another in Cali), and three visits (pre-treatment and twice during the course of treatment). Here, and based on the results shown in Supplementary Fig. 5, we modeled the ‘clinic’ variable in SVA and used the surrogate variables-modified counts to generate the PCA plots. The analyses were performed using all annotated protein-coding genes following filtering for low counts, CPM normalization, and log₂ transformation. In the PCA plots, the first two principal components are shown on the X and Y-axes, respectively, with the proportion of total variance attributable to that PC indicated. Each sample is represented as a single point with color indicating either the clinic (A), or clinic, treatment and outcome (B–E). The inner and outer colored ellipses represent the 90% and 95% confidence intervals, respectively. These are only shown for classes containing more than 3 samples.

The relative contribution of metadata factors to the transcriptional profile was further examined via a series of principal component (PC) and surrogate variable (SV) loading analyses. This was performed by taking the first 5 PCs of the normalized (CPM and filtered) data and collecting their F-statistics and associated P-values with respect to a series of empirically observed important metadata factors: clinic, donor, and visit number. PCs and factors with high F-statistics and low P-values were deemed noteworthy. The normalized data was then passed to SVA with cure/TF status as the variable of interest, and the F and P-values were collected. We found that the most prominent SV in the expression data was ‘clinic’ (SV4 and SV5), followed by modest F-statistic/P-values for ‘donor/participant’ (SV4) and, to a much lesser extent, ‘visit number’ (SV1 and SV2). The same values were then collected from the SVA-adjusted transcriptional profile (Supplementary Data 2). Indeed, the effect of ‘clinic’ increased significantly relative to the other two factors (as shown by higher F and P-values in the Post-SVA PC1), indicating a persistent batch effect associated with the clinic.

Based on the significant batch effect introduced by the patient recruitment clinic and the skewed representation of cured patients recruited in Cali (9 of 10 CL patients recruited in this site cured), we excluded all samples obtained from patients in the Cali clinical site from the initial transcriptomic analyses and biomarker discovery. Therefore, the transcriptomic variance associated with therapeutic outcome was analyzed on SVA-adjusted data sets from samples derived from patients recruited in Tumaco. After subsetting to include only samples collected from Tumaco and carrying out similar loadings analyses as described above (Supplementary Data 2), the effect of ‘visit’ continued as minor, showing that the expression changes in the samples are modest between samples collected Pre-Tx, Mid-Tx or at End-Tx. The largest effect remained attributable to the ‘donor’.

Transcriptomes from lesion biopsies obtained pre-treatment corroborate heightened tissue cytolytic activity in treatment failure patients

A recent comparative transcriptome analysis of Pre-Tx lesion biopsies from CL patients infected with L.V. braziliensis showed a gene signature of heightened cytolytic activity in lesions from TF patients, compared to those who cured¹⁵. To explore the congruity of those findings with infections with other L. Viannia species (L.V. panamensis), we examined CL lesion transcriptomes in our study cohort. Overall, the transcriptomic profiles of lesion biopsies from patients who cured were indistinguishable from those of patients with TF (Supplementary Fig. 5E), even following SVA (Fig. 3E). As expected, only few genes (n = 28) were differentially expressed (DE) (P < 0.05; |log₂FC | ≥ 1), 17 of which were up-regulated and 11 down-regulated in patients with TF, compared to cures (Supplementary Data 3a). Notably, among up-regulated genes, a signature of increased cytolytic activity (GZMB, NCR1, SH2D1B, PRF1, KLRC1, GNLY, FGFBP2, KIR2DL4, CCL3 and CCL4) was found in tissue samples from TF patients (Supplementary Data 3b), consistent with previous findings from TF patients infected with L.V. braziliensis¹⁵. The minimal difference in the global transcriptomes of skin lesions from cured and TF patients was expected since bulk transcriptomes from complex multicellular tissues are often skewed to reflect the most abundant or transcriptionally active cells within the sample. Nevertheless, the fact that a clear transcriptomic signature of enhanced cytolytic activity was detected in TF suggests that systemic differences leading to the activation and/or recruitment of cytotoxic Natural Killer (NK) and CD8⁺ T cells could be contributing to this enhanced inflammatory state.

Transcriptomic profiles of innate immune cells do not change over the course of treatment

Our understanding of the participation of innate immune responses in the outcome of antileishmanial therapy is almost exclusively limited to the role of macrophages as primary host cells for the parasite. However, mounting evidence shows that other innate cells, including Nφ, and more recently Eφ and NK cells, participate in the inflammatory responses that contribute to CL immunopathology, and thus their functions are relevant to the outcome of treatment¹⁸.

Among the hallmarks of innate immune cell functions are the velocity and robustness of their elicited responses, which in turn require tight mechanisms of control to avoid host injury. To explore the dynamics of these responses and their participation in therapeutic responsiveness, we analyzed samples collected from CL patients Pre-Tx, Mid-Tx, and at End-Tx. PCA plots of the transcriptomes of individual cell types did not reveal any significant clustering of samples based on visit (i.e., over the course of treatment) (Fig. 4A–C). Furthermore, a correlation analysis of DE genes (DEGs) derived from SVA-adjusted data using two different models, one which explicitly included visits in the DE model and the other which did not, showed strong and significant correlations for the log₂FC values in Mo, Eφ, and Nφ (r > 0.9), as well as for the respective P-values (ρ > 0.82) (Fig. 4D–F). These data rule out a substantial effect of including visits in the DE model for cure vs. TF. Based on the above, we opted to group transcriptome samples from all time points from each cell type for the identification of biomarkers and transcriptional signatures of cure and TF.

PCA for all samples were performed using log₂ transformed, quantile and CPM-normalized, and filtered counts of annotated protein-coding genes (A–C). The first two principal components are shown on the X and Y-axes, with the proportion of total variance attributable to that PC indicated. D–F show correlation plots of the log₂FC values for DEGs identified using DESeq2 with the default model including cure/TF as factors with SVA-adjustment (x-axis), and SVA-adjusted data where cure/TF and visit (Pre-Tx, Mid-Tx, and End-Tx) were considered as factors. Pearson correlation coefficient (r) is shown for DEGs and Spearman correlation coefficient (ρ) for associated P-values. (G–I) Correlation plots of the log₂FC values for DEGs using Dream (x-axis) and DESeq2 using SVA-adjusted data. Pearson correlation coefficient (r) is shown for DEGs and R² of the linear regression.

Because we sampled the same patients before treatment and twice during the course of several weeks post-treatment, we wanted to rule out that the multiple data points were not amplifying the effects observed. To that end, we included “patient” as a random effect variable in a series of linear mixed models (LMM) using the dream¹⁹ functionality from variancePartition²⁰. Our LMM included visit, final outcome, and cell type as fixed effects; and donor as a random effect. We then compared the DE outputs of dream (log₂FC and P-values in the contrast of Cure vs. TF; Supplementary Data 4) against the DESeq2 data. High correlations of log₂FC values (≥ 0.84) were observed for all cell types (Fig. 4G–I), indicating a strong similarity among DEG lists obtained with the two analytical approaches. Up to 40% of genes deemed significant by DESeq2 in any of the three cell types were also significant using dream. Network analyses using as input the DEGs resulting from dreams also yielded consistent results (genes and pathways) to those derived from the DESeq2 data (Supplementary Fig. 6).

Monocyte, neutrophil and eosinophil transcriptomes from CL patients who cure differ from those with TF

Examination of Nφ, Mo, and Eφ transcriptomes showed a clear clustering of samples by treatment outcome, in each of the three cell types, revealing a distinct separation of samples along the first principal component (PC1) (Fig. 5A). On average, PC1 explained 18% of the variance across the three cell types. A DE analysis in each of the three cell type transcriptomes revealed 3 to 4 times more significant DEGs between cures and TF in Mo, when compared to Eφ and Nφ (Supplementary Data 5), consistent with higher transcriptional activity of Mo^21,22. Those DE profiles are represented in the form of volcano plots (Fig. 5B). When a |log₂FC | ≥ 1 threshold was established, a comparable number of DEGs was observed between cures and TF among the three cell types: 191 in Eφ, 160 in Nφ and 112 in Mo (Supplementary Data 5 and Fig. 5C).

PCA plots, after batch correction estimated by SVA, of transcriptomes collected at all-time points from monocytes (Mo), eosinophils (Eφ), and neutrophils (Nφ) (A). Volcano plots of DEGs using DEseq2 (Wald test and its Benjamini-Hochberg adjusted P-value < 0.05 and |log₂FC | ≥ 1) for each cell type, with genes up-regulated in TF vs. cures colored pink, and those down-regulated in TF vs. cures labeled purple. The top 10 up- and down-regulated genes are labeled (B). Venn diagram showing the intersection of common DEGs (P < 0.05 and |log₂FC | ≥ 1) between cell types (C). The inner and outer colored ellipses represent the 90% and 95% confidence intervals, respectively. These are only shown for classes containing more than 3 samples.

In Eφ, a robust induction of IFNα/β signaling was revealed by gene up-regulation of innate immune receptors, signaling molecules, transcription factors, and regulators of signaling pathways (Supplementary Data 5). The gene encoding IRF7 was up-regulated in TF Eφ. Together with IRF3, IRF7 is the canonical transcriptional regulator of Type I IFNs²³. IFIH1 (gene encoding MDA5), a RIG-I-like receptor and cytoplasmic sensor of dsRNAs, was also up-regulated in TF. MDA5 has been demonstrated to be an amplifier of innate immune responses and associated with autoinflammation²⁴. Genes encoding downstream effectors (including OAS1, OAS2, OAS3, OASL, BST2, MX1, MX2, IFI6, XAF1, GBP2, and IFI27), and pathway regulators (IFIT5, ISG15, RSAD2, USP18, HLA-G, DHX58, and DDX60) were found in enriched gene categories in Eφ (Fig. 6A, B), substantiating elements of a Type I IFN signature in TF. Type I IFN response pathway genes were also enriched in Nφ transcriptomes from TF (Supplementary Data 6 and Supplementary Data 7). Notably, genes encoding HERC6, IFI44L, ISG15, USP18, IFI27 and DDX60 (Fig. 6C, D) were also expressed at higher levels in TF Nφ. Consistently, Mo transcriptomes from TF patients were also enriched in mRNAs encoding Type I IFN-related genes (Fig. 6E, F), suggesting synergistic innate cell functions towards hightened type I IFN inflammation in TF. Downregulation of IL1R1 and IL1R2 (decoy) receptors was observed in Mo from TF patients, and downregulation of molecules related to MHCII antigen presentation was found in all cell types from TF patients (Supplementary Data 5, Supplementary Data 6 and Fig. 6A, C, E). Genes involved in wound healing and cell proliferation (HBGEF, EGR1, and EGR3) were down-regulated in Eφ of TF patients. In Mo from TF, down-regulated antimicrobial peptide genes (including CTSG, LTF, CAMP, DEFA3, and LCN2) and immune receptor activity genes (including IL2RB, IL1R2, HLA-DQA1, HLA-DQB1, IL1R1, and CXCR4) were significantly enriched. Altogether, these results suggest that mechanisms underlying TF are associated with an enhancement of Type I IFN signaling, a dampening of antimicrobial effectors (antimicrobial peptides), and functions linking the innate and adaptive immune systems (antigen presentation).

Panels A, C, E, and G. Significantly DEGs (P_adj < 0.05 and |log₂FC | ≥ 1) were used as input for network analyses. STRING V12.0 was used to construct networks based on co-expression, databases, and experiment terms, with line thickness representing the confidence of the interaction. Enriched categories shown overlaying the networks are those selected from KEGG, GO, or Reactome (one-sided Fisher test with g:SCS-corrected P-values) terms with FDR < 0.05 and strength > 1. The blue border in nodes represents up-regulated genes in TF compared to cures; red borders depict downregulated genes. The intensity of the border color reflects the magnitude of the DE in our dataset. Genes belonging to the Type I IFN cluster are grouped under a blue shade, MHCII in red, and others in gray. Panels B, D, F, H. The same significant DEGs were used for over-representation analyses by gProfiler2. Terms with adjusted P-values ≤ 0.05 (g:SCS method in gProfiler) from the GO biological processes were clustered by term similarity, and the resulting trees were plotted. Groups of nodes are colored by their parent GO term, the size of terminal nodes denotes the number of genes observed in the category, and the color of each node represents the adjusted significance of the observed over-representation.

Twelve DEGs were common among all cell types: IFI44L, IFI27, PRR5, PRR5-ARHGAP8, RHCE, FBXO39, RSAD2, SMTNL1, USP18, AFAP1, SIRPG and OTOF (Fig. 5C). Although more than 70% of significant DEGs with |log₂FC | > 1 were unique to each cell type (underscoring a cell-specific response) gene enrichment analyses showed common significantly enriched features between the different innate cells (Supplementary Data 6, Fig. 6B, D, F). Among up-regulated genes, enrichment of type I IFN responses was common to Mo, Eφ, and Nφ from TF patients, and this was also supported by GSVA (Supplementary Data 7). No similarities in gene enrichment analysis were found for down-regulated genes among cell types, with the exception of MHCII-related genes (Supplementary Data 6).

Considering the overlap of DEGs among cell types, we evaluated whether an IFN gene signature could be detected in a combined analysis of transcriptomes derived from all innate cell types. As observed in the PCA before SVA (Supplementary Fig. 7A), most of the variance in the data (81% for PC1 and PC2 combined) could be attributed to the cell type. When SVA was used (Supplementary Fig. 7B), a partial separation of transcriptomes was observed between cells isolated from CL patients who cured versus those with TF. Differential expression analyses revealed 210 genes ( | log₂FC | ≥ 1; 137 upregulated in TF compared to cures and 73 downregulated (Supplementary Data 5). Consistently, the IFN and MHC signatures identified in the analyses of each cell population were also observed in the combined analysis (Fig. 6G, H).

In a complementary approach to DE analyses, we employed weighted gene co-expression network analysis (WGCNA) to identify co-expressed gene clusters associated with therapeutic response. The Mo data resulted in four gene modules with a significant association. The four modules comprised 478 genes and were enriched in genes encoding ribosomal, mitochondrial and DNA nuclear activity and IFNα/β-inducible proteins (Supplementary Data 8and Supplementary Fig. 8). In Nφ, three modules comprising 581 genes were associated with therapeutic response and were enriched in cellular pathways similar to those found in our DE analyses: Type I IFN, MHC-II, DNA nuclear activity, glycolytic metabolism, cell migration, vesicular transport, cell death and the immunoproteasome. For Eφ, both associated modules comprised 200 genes, with an enrichment of genes related to Type I IFN responses, consistent with DE analysis, and Nφ and Mo WGCNA. Of the total genes contained within the significant modules from the three cell types, 36 were shared, and all were up-regulated in TF. Among those, 24 were related to the IFN responses (Supplementary Fig. 8C, D), and this was consistent with the DE analyses. Prominent in WGCNA were STAT1 and STAT2, two known transcription factors involved in IFN signaling and expression of interferon-stimulated genes (ISGs). Both transcription factors were up-regulated in TF, suggesting coordinated and sustained IFN effector functions elicited downstream of IFN receptor ligation. Together, WGCNA and DE analyses support the participation of a systemic pro-inflammatory environment sustained in TF, mediated in part by Type I IFNs.

To explore whether these IFN signatures observed in innate cells were also detected at the lesion site, we compared CPM values of any differentially expressed ISGs identified in any of the innate cell populations analyzed, including those common to all cell types. Significantly higher expression of OAS1, OASL, GBP2, MX1, DDX60, IFIT5 and XAF1 was observed in lesion biopsies from TF patients (Supplementary Fig. 9A). Higher expression of MX1 and OASL was corroborated in an independent transcriptomic data set from lesion biopsies of other 11 patients previously published by our group²⁵ (Supplementary Fig. 9B). In a combined analysis of all lesion biopsy transcriptomes (Supplementary Fig. 9C; cures n = 14, TF n = 11), OAS1, OAS3, OASL, RSAD2 and MX1 were significantly higher in TF biopsies, suggesting involvement of local and systemic IFN responses in treatment failure.

IFNα/β stimulated genes constitute a hallmark signature of TF in CL

Based on the functional commonalities between the cell-specific transcriptional profiles, we recognized a common innate gene signature, from which we were able to identify biomarkers that predict TF. From the gene lists derived from our WGCNA and DE analyses, we selected the 12 common DE genes in Mo, Nφ and Eφ: IFI44L, IFI27, PRR5, PRR5-ARHGAP8, RHCE, FBXO39, RSAD2, SMTNL1, USP18, AFAP1, OTOF and SIRPG. Each of those common genes was up-regulated in TF, with the exception of AFAP1, which was down-regulated. We, therefore, constructed a composite score for each patient based on this innate gene signature of up-regulated genes (which excluded AFAP1). RPKM values were extracted from all Pre-Tx samples for each cell type. Those were selected for their predictive value in guiding early therapeutic interventions. RPKM values were used to construct raw and normalized scores (Z-score), the latter to account for possible outliers within groups. Both raw and Z-scores were higher in TF patients, and this difference was statistically significant for the scores derived from Nφ, and for the Mo + Nφ composite scores (Fig. 7A–D and Supplementary Fig. 10). Consistently, raw and Z-scores from Pre-Tx Nφ samples and the Mo + Nφ scores, were significantly predictive of the therapeutic response (Fig. 7E, F).

RPKM data of the 11 up-regulated innate signature genes from Pre-Tx samples was used to construct neutrophil (A, C) and multi-cell (monocyte + neutrophil panels B, D) composite scores. Raw (A, B) and normalized scores (C, D) are shown. Statistical significance was evaluated by analysis of variance. P-values from two-tailed analyses are shown. Representation of receiver operating characteristic (ROC) curves for raw (filled line and black circles) and normalized (dashed line and open circles) scores from neutrophils (E) and a combined score of monocytes plus neutrophils (F). ROC curves represent the area under the curve (AUC) of the false positive rate versus the true positive rate. P-values, AUC, and cutoff values based on Youden’s J statistic are shown in the graph plots. Mono+Neutro: combined scores of monocytes and neutrophils.

We next applied a machine learning approach to carry out a comprehensive analysis of all data on hand and complement our focused DE analyses. The expression matrices of all samples obtained from Tumaco or both clinics (Tumaco and Cali) were split into 10 rounds of training (40%) and testing sets (60%), through random partitioning and cross-validation. The training sets were used to generate K-nearest neighbors (KNN), logistic regression (GLM), gradient boost (XGBoosted GLM), and random forest models²⁶. Following an evaluation of the resulting models by comparing the predictions of the training data to the known clinical outcomes (Supplementary Data 9), we used those models to predict clinical outcomes in the test partition and evaluated their performance with an emphasis on a dataset that only included Pre-Tx samples as these would be the most translatable for clinical and therapeutic decision making. Overall, the GLM models showed better performance. Using data from all blood cell types collected from all visits in the Tumaco clinic only, specificity, sensitivity, and accuracy, were 0.83, 0.84, and 0.82, respectively (Supplementary Data 9). When the models were restricted to include only Pre-Tx samples, the metrics were impacted (as an example, GLM metrics changed to specificity = 0.78, sensitivity = 0.70, accuracy = 0.69). A second round of the dual analysis described above included samples from both clinics, substantially improving the performance of all models with accuracy up to 0.89, sensitivity 0.87, specificity 0.95, and AUC 0.85 (for GLM) (Supplementary Data 9). As expected, the addition of samples from both study sites, and from all subsequent visits enhanced the performance of all models (Supplementary Data 9), possibly overcoming the clinic-based batch effect observed for DE analyses.

To explore the relationship between the predictive features identified by ML and DE, we conducted a “pseudo-bulk” DE analysis of all innate cell samples obtained at the three time points from patients from Cali and Tumaco (Fig. 8A, B). A total of 282 genes from the combined DE analysis in the Cure/TF contrast were observed, among which were the 12 innate biosignature genes (Supplementary Data 10, Fig. 8C). The top 300 most variable importance genes for each of the four ML algorithms used was contrasted against the DE gene set (Supplementary Data 10). Interestingly, the most significant overlap between ML and DE analyses resulted from the KNN model, sharing 18 genes, five of which also feature in the innate biosignature (IFI27, USP18, OTOF, SIRPG, and RSAD2). Although the three other algorithms performed better than KNN in terms of specificity, sensitivity, and accuracy metrics, they fared worse in their overlap with the DE analyses. This likely results from the requirement of significantly larger datasets for ML compared to DE analyses.

PCA for all samples collected from patients in Cali and Tumaco (Mo, Nϕ, and Eϕ, collected Pre-Tx, Mid-Tx, and at End-Tx) were performed using log₂ transformed, quantile and CPM-normalized, and filtered counts of annotated protein-coding genes (A), SVA-adjusted data considering cell-type for variance adjustment (B). The first two PCs are shown on the X and Y-axes, with the proportion of total variance attributable to that PC indicated. The inner and outer colored ellipses represent the 90% and 95% confidence intervals, respectively. C Volcano plot of DEGs, using DEseq (Wald test and its Benjamini-Hochberg adjusted P-value < 0.05 and |log₂FC | ≥ 1), highlighting in labels the top 10 up- or down-regulated genes in TF compared to cures.