Study design and patient cohort

Patients were enrolled between August 2018 and September 2021 at the University Hospital Cologne, division of infectious diseases. EPTB diagnosis was made in patients with the detection of M. tuberculosis complex (culture positivity and/or PCR positivity) from bodily secretions or tissue. In three culture-negative patients, a clinical diagnosis was made based on epidemiological exposure combined with physical findings, radiographic findings, and appropriate histopathological findings (e.g., necrotizing granuloma, acid-fast staining positivity). Except for one patient with missing data, all patients showed positive results in TB-specific IFN gamma release assays (IGRA), confirming previous infection with M. tuberculosis complex. All patients underwent systematic baseline imaging, including chest X-ray, abdominal ultrasound, and lymph node ultrasound. For specific clinical findings, suspected disseminated disease, or inconclusive chest X-ray results, additional computed tomography (CT) or magnetic resonance imaging (MRI) was performed to assess the extent of the disease within a staging protocol. Radiological data were centrally reviewed to document organ involvement, lesion characteristics, and evidence of dissemination, defined as involvement of two or more non-contiguous sites. Baseline assessments included symptom history, physical examination, and measurement of inflammatory markers such as CRP, leukocyte count, and ESR. Patients suffering exclusively from central nervous system (CNS) TB were excluded from the study, as CNS TB primarily affects the neuroimmune system.

Detailed information on all patients included in the study are provided in Supplementary Table 1. Sample size and epidemiological imbalances of EPTB did not allow for matching according to the sex of patients. IGRA negative healthy volunteers were used as controls.

Peripheral blood mononuclear cell (PBMC), plasma and serum isolation of whole blood

PBMCs were purified using density gradient centrifugation using Cytiva Ficoll-Paque® (Fresenius, Bad Homburg, Germany) and Leucosep^TM tubes (Greiner, Kremsmünster, Austria). At this step plasma was taken off and stored at − 80 °C. Afterwards, PBMCs were washed several times with PBS and stored in freezing media (FBS + 10% DMSO) at − 150 °C.

Whole blood RNA isolation

PAXgene Blood RNA Kit (Preanalytix, Hombrechtikon, Germany) was used according to the manual provided by the vendor. Briefly, PAX tubes were thawed at room temperature, centrifuged and washed with RNase-Free Water. Next, pellet is resuspended in resuspension buffer and binding buffer and proteinkinase K. Incubation times were performed as recommended and the spin column procedure was performed as suggested by the manual. Upon elution, RNA was incubated for 5 min at 65 °C and RNA was stored at − 80 °C.

Bulk RNA sequencing

Sequencing library prep was performed with in total input of 100 ng RNA according to the NEBNext Ultra RNA library prep protocol (New England Biolabs, Ipswich, MA, USA). Next, the RNA libraries were validated with Tape Station 4200 (Agilent Technologies, Santa Clara, CA, USA). Library quantification was performed using the KAPA Library Quantification Kit (VWR International, Radnor, PA, USA) and the 7900HT Sequence Detection System (Applied Biosystems, Foster City, CA, USA). Sequencing was done with NovaSeq6000 sequencers (Illumina, San Diego, CA, USA) with a PE100bp read length, aiming at 50 M reads/sample. Quality control of samples was performed using fastQC (v0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and multiQC⁴⁴ (v1.14) (https://multiqc.info/). Sequenced reads were aligned against the Gencode human reference genome v33 using STAR⁴⁵ (v2.7.10b) (https://github.com/alexdobin/STAR, https://doi.org/10.1093/bioinformatics/bts635) and reads were summarized using no strandness information. The alignment pipeline was executed using SnakeMake (v7.20.0) (https://snakemake.github.io/, https://doi.org/10.12688/f1000research.29032.2).

Bulk RNA-seq data preprocessing

The following analysis steps were performed in R⁴⁶ (v4.1.0) and R Studio⁴⁷ (v1.4.1717). After importing the count matrix, gene with less than 10 counts in less than 4 samples were excluded, resulting in 20,206 genes kept in the analysis. Normalization of the count matrix was computed with DESeq2⁴⁸ (v1.32.0), a variance stabilizing transformation applied using the DESeq2 vst function at default settings and batch-corrected for different sequencing time points with limma⁴⁹ (v3.48.3).

Differential expression analysis in bulkRNA-seq data

Differential expression analysis based on DESeq2 was performed by adjusting p-values according to independent hypothesis weighting from IHW⁵⁰ (v1.20.0) at default settings and applying empirical Bayes shrinkage estimators from apeglm⁵¹ (v1.14.0). Genes with an abs(fold change) > 2 and an adjusted p-value 

Gene set enrichment analysis

Gene Ontology (GO)^52,53, Kyoto Encyclopedia of Genes and Genomes (KEGG)^54,55,56 and Molecular Signature Database Hallmark gene set^57,58 enrichment was performed the clusterProfiler⁵⁹ (v4.4.4) enricher function. Gene-Term associations for the respective databases were downloaded from the Molecular Signature Database⁵⁷ (v2023.1).

Gene co-expression network analysis

Gene co-expression network analysis was performed using hcocena^31,60 (v1.0.0). The 9010 most variable genes of the normalized, batch-corrected data were used as input and gene-gene correlations were computed using Pearson’s correlation. Gene pairs with a Pearson’s correlation coefficient lower than 0.751 were excluded from the network, resulting in a network with 5042 genes, 178,105 edges and an R²-value of 0.8. Leiden clustering with a resolution of 1 identified 10 modules with a minimum size of 50, containing a total of 4965 genes.

Cell type deconvolution

Cell type abundances of the batch-corrected data were determined by CIBERSORTx⁶¹ (https://cibersortx.stanford.edu/) with default parameters using BD Rhapsody datasets of cohort 2 from EGAS00001004571³² filtered for whole-blood samples and healthy controls as reference. From each cell type, 1000 cells were sampled if available, and the normalized data was used to generate the single-cell reference matrix.

Classifier set up and signature assessment

To create a balanced dataset of control and EPTB patients, 13 additional control samples were pre-processed. Genes with less than 10 reads in less than 2 samples were excluded, resulting in a count matrix of 19,979 genes. Further pre-processing steps were performed as previously described. The combined datasets was filtered for the intersection of genes included in both datasets, consisting of 18,663 genes from 22 control and 29 EPTB patients and genes were ranked per patient. To assess the classification performance, a random forest classifier was constructed with caret⁶² (v7.0-1) and randomForest (v4.7-1.2)⁶³. The expression values of the combined datasets were rank transformed and the data split into a training and test datasets in a 2:1 ratio. RF hyper-parameters were set to default and the ‘mtry’ parameter was set to \(\root{{2}}\of{n-1}\) where n is the signature length. The performance was estimated using a leave-one-out cross validation and was evaluated with MLeval (v0.3)⁶⁴ and by calculating the area under the receiver operating characteristic curve (AUROC) generated with plotROC⁶⁵ (v2.3.1) and pROC (v1.18.8)⁶⁶. Gene set variation analysis was performed with GSVA⁶⁷ (v1.40.1) without kernel estimation of the cumulative density function on the ranked combined datasets to assess the classification performance for each inflammatory stage. Significance of the results was calculated with a Wilcoxon test followed by a Benjamini-Hochberg adjustment.

GSVA-based gene prioritization approach

Group-specific EPTB signatures were generated by either using the intersection of identified DEGs between the different EPTB inflammatory stages and control samples or between a selected EPTB inflammation immunotype and the other inflammation immunotypes following an adapted optimization strategy published in Knoll et al.⁶⁸. In short, the intersection of DEGs was ordered by the (average) log2 fold change, and ties were solved based on mean expression level. An iterative GSVA was then performed, starting with the gene with the highest (average) log2 fold change, adding the gene with the next highest log2 fold change with each iteration. For each iteration, a Wilcoxon test was performed between the different EPTB inflammatory stages and control samples or between a selected EPTB inflammatory stage and the other inflammatory stages, and the average p-value was calculated. Genes which increased the p-value were excluded from the signature. The final signature consisted of the gene set with the lowest p-value.

Enrichment of signatures in other diseases

To assess the disease specificity of the core signature, bulk RNA-seq and microarray datasets of other viral, bacterial or pulmonary diseases were analysed. Datasets were selected based upon the following criteria: (i) The data was generated from whole blood, (ii) sample collection protocols were as similar as possible to the described above (i.e., blood collection tubes), and (iii) included healthy controls as a reference. The following datasets were considered: Bloom et al.³⁹ (GSE42834 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42834]), Zhai et al.³⁸ (GSE68310 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68310]), Aschenbrenner et al.³⁷ (EGAS00001004503 [https://ega-archive.org/studies/EGAS00001004503]) and Herwanto et al.³⁶ (GSE154918 [https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE154918]) together with the EPTB datasets. Each datasets was pre-processed as previously described and transformed into a log2-space while genes not included in all studies were excluded from the datasets. The final datasets consisted of 494 samples and 11,778 genes. All datasets were integrated with the sva⁶⁹ (v3.40.0) ComBat function prior to performing a classification benchmark. GSVA between the different disease states and controls was performed without kernel estimation of the cumulative density function. Significance of the results was calculated with a Wilcoxon test followed by a Benjamini-Hochberg adjustment. For the EPTB core signature and for each PTB signature, an RF classifier was constructed based on the rank-transformed integrated datasets using the architecture described above. Classification performance was evaluated based on an overall score reflecting the sum of the AUROC, specificity, precision and recall.

BD Rhapsody blood single-cell RNA-seq and AbSeq

Whole transcriptome analyses and surface protein expression, using the BD Rhapsody Single-Cell Multiomics System (BD, Biosciences) were performed on PBMCs of 3 selected donors of each stratified EPTB group and 3 healthy controls.

Cells from four samples per condition were labeled with sample tags (BD Human Single-Cell Multiplexing Kit) following the manufacturer’s protocol. Briefly, a total number of 1 × 10⁶ cells were resuspended in 90 μl of Stain Buffer (FBS) (BD PharMingen). The sample tags were added to the respective samples, and the samples were incubated for 20 min at room temperature. After incubation, 200 μl stain buffer was added to each sample and centrifuged for 5 min at 300 × g and 4 °C. After washing, cells were resuspended in 500 μl of cold Sample Buffer (BD, Biosciences) and counted using a Neubauer Hemocytometer. Based on the counting, labeled samples were pooled equally.

After pooling, the cell suspension was stained with the BD AbSeq Immune Discovery Panel (BD, Biosciences, Cat#25970) to assess surface marker expression. Reconstitution and cell staining were performed according to the manufacturer’s protocol. After staining, the cells were resuspended in 300 μl of cold Sample Buffer, and a second cell counting step was conducted. The concentration was adjusted to achieve 120,000 cells in 1300 μl of cold Sample Buffer. The pooled samples were then split in two for super-loading two BD Rhapsody cartridges, with approximately 60,000 cells each.

Single cells were isolated using Single-Cell Capture and cDNA Synthesis with the BD Rhapsody Express Single-Cell Analysis System according to the manufacturer’s recommendations (BD Biosciences). cDNA libraries were prepared using the mRNA Whole Transcriptome Analysis (WTA), AbSeq and Sample Tag Protocol (BD Biosciences). The final libraries were quantified using a Qubit Fluorometer with the Qubit dsDNA HS Kit (ThermoFisher), and the size-distribution was measured using the Agilent high-sensitivity D5000 assay on a TapeStation 4200 system (Agilent Technologies). Sequencing was performed in paired-end mode (R1 85 bp and R2 215) on a NovaSeq 6000 with NovaSeq 6000 S2 Reagent Kit (300 cycles) chemistry. After demultiplexing of bcl files using bcl2fastq2 (v2.20) from Illumina and quality control, a barcode whitelist provided by BD Biosciences was used to filter the paired-end scRNA-seq reads for valid cell barcodes. Cutadapt⁷⁰ (v1.16) package (Martin, 2011) was used to trim adapter sequences and to filter reads for a PHRED quality score of 20 or above. Next, STAR⁴⁵ (v2.6.1b) was used for alignment against the Gencode v33 human reference genome and Dropseq-tools (v2.0.0) were used to quantify gene expression and collapse to UMI count data (https://github.com/broadinstitute/Drop-seq/). For SampleTag oligo-based demultiplexing of single-cell transcriptomes and subsequent assignment of cell barcodes to their sample of origin, the respective multiplexing tag sequences and AB-seq sequences were added to the reference genome and quantified, as well⁷¹.

scRNA-seq data pre-processing

The following analysis steps were performed in R⁴⁶ (v4.1.0) and R Studio⁴⁷ (v1.4.1717). scRNA-seq analysis was performed using Seurat^72,73 (v4.0.4). Cells of each cartridge were demultiplexed in a two-step process. First, cells that were considered as doublets based on the Seurat HTODemux function at default settings were discarded prior to genotype-based demultiplexing using Vireo⁷⁴ (v0.5.6) and cellsnp-lite⁷⁵ (v1.2.2). Only cells defined as singlets by the Seurat HTODemux function or cells that could be assigned with a 90 % probability to its respective donor based on the single-nucleotide polymorphism (SNP) information were kept. In addition, cells were excluded if they had 1) more than 35 % mitochondrial genes approximately equal to three mean absolute deviations above the median according to scater⁷⁶ (v1.20.1), 2) more than 5 % hemoglobin genes, 3) less than 465 and 338 UMIs, as well as 412 and 328 genes based on the inflection point of the cells ranked by either the number of UMIs or genes of each cartridge, respectively, and 4) a prominent multi-lineage marker gene expression. Genes expressed in less than in 5 cells were excluded. After QC, a total of 61,537 PBMCs expressing 28,604 genes were analysed. Normalization, scaling and dimensionality reduction of the WTA and AbSeq data were computed with Seurat functions. The gene expression values were normalized by multiplying the total UMI counts per cell by 10,000 (TP10K) and applying a natural log transformation by log(TP10K + 1). Subsequently, the normalized data was scaled as well as centered and a regression between each gene and donor was performed to account for the intra-individual heterogeneity of each donor. Following a variance stabilization using the ‘vst’ method, PCA of the 2000 most variable features was performed to reduce the dimensionality of the data. Protein expression data of each cell was normalized by applying a centered log-ratio transformation prior to scaling and centering. A two-dimensional, multimodal representation of the data was calculated by generating a weighted nearest neighbor (WNN) graph of the first 30 PCs of the gene expression data and the normalized protein expression data with the shared nearest neighbor (SNN) pruning cut-off set to 1/25. Nearest neighbors were then used to calculate a UMAP. Cells were clustered based on the WNN graph using the Louvain algorithm at a resolution of 0.3. Cluster identities were determined using the Seurat FindAllMarkers function using the Wilcoxon rank sum test, defining cluster marker genes as genes expressed in at least 25 % of the cells of either of the two compared populations, with a difference of at least 10 % between the populations, and a log fold change > 0.25. During log fold change calculation, the mean was computed as implemented in Seurat⁷⁷ (v5):

where \(x\) is defined as the respective gene and \(n\) as the number of cells. Subsequently, cells were annotated in a three-step process by overlapping cluster marker genes and protein marker expression with literature-known cell-type marker genes as CD14 Monocytes (CD14, S100A8, S100A9), CD16 Monocytes (CD16, CDKN1C, FCGR3A), mDCs (CD1C, FCER1A), pDCs (CLEC4C, IRF7, LILRA4), CD4 T cells (CD3, CD4, TCF7), CD8 T cells (CD3, CD8, CD8A), NK cells (CD56, NKG7, GZMB), B cells (PAX5, IGHD, CD22, CD19), Plasmablasts (CD27, IGHA1, IGHG1, IGHA2), Megakaryocytes (TUBB1, ITGA2B, PPBP), Progenitors (CD34, GATA2), Proliferating cells (TYMS), and Mixed (CD19, CD14).

Selection and annotation of monocytes

Monocytes were selected and annotated in an iterative process. Genes expressed in less than 5 cells were excluded from the analysis. In the first iteration, the data was subset to all cells defined as either CD14 monocytes, CD16 monocytes or mDCs. Then, the subset was normalized, scaled and centered as described above. For dimensionality reduction of the gene expression data, a PCA was performed of the 2000 most variable genes following a variance stabilization using the ‘vst’ method. To create a two-dimensional, multimodal representation of the data, a WNN graph of the first 15 PCs and the normalized and scaled protein expression data was generated with an SNN pruning cut-off set to 1/25. Nearest neighbors were then used to calculate a UMAP. Cells were clustered based on the WNN graph using the Louvain algorithm at a resolution of 0.94. Clusters with high expression of literature-known marker genes of non-monocytes (such as NK cells, T cells or mDCs) in addition to the typical monocyte markers were excluded from further analysis, resulting in 20,395 total monocytes. After the clean-up, a second iteration of the previous pre-processing steps was performed at final clustering resolution of 0.9 and cluster identities were determined using the Seurat FindAllMarkers function using the Wilcoxon rank sum test defining cluster marker genes as genes expressed in at least 25 % of the cells of either of the two compared populations, with a difference of at least 10 % between the populations, a log fold change > 0.25. During log fold change calculation, the mean was computed as described above. Subsequently, cells were annotated in a two-step process based on cluster marker genes and protein marker expression as CD14⁺HLA-DR^hi (HLA-DQA1, HLA-DPB1, HLA-DQB1), CD14⁺S100⁺ (S100A12), CD14⁺CCL⁺ (CCL3, CCL3L3, CXCL8), CD14⁺CD163⁺ (CD163, F13A1), CD14⁺THBD⁺ (THBD, THBS1, RGCC), CD14⁺IFI⁺ (IFI44, IFIT3, IFIT2), CD14⁺CD16^hiC1Q⁺ (C1QA, C1QB, C1QC), and CD14^–CD16^hi (CDKN1C, FCGR3A).

Selection and annotation of T and NK cells

T and NK cells were selected and annotated in an iterative process. Genes expressed in less than at least 5 cells were excluded. In the first iteration, the data was subset to all cells defined as either CD4 T cells, CD8 T cells or NK cells. Then, the subset was normalized, scaled and centered as described above. For dimensionality reduction of the gene expression data, a PCA was performed of the 2000 most variable genes following a variance stabilization using the ‘vst’ method. To create a two-dimensional, multimodal representation of the data, a WNN graph of the first 15 PCs and the normalized and scaled protein expression data was generated with an SNN pruning cut-off set to 1/25. Nearest neighbors were then used to calculate a UMAP. Cells were clustered based on the WNN graph using the Louvain algorithm at a resolution of 1.1. Clusters with high expression of literature-known marker genes of non-T/ NK cells or multiple cell types were excluded from further analysis, resulting in 29,662 total T and NK cells. After the clean-up, a second iteration of the previous pre-processing steps was performed with a final clustering resolution of 0.7 and cluster identities were determined using the Seurat FindAllMarkers function using the Wilcoxon rank sum test defining cluster marker genes as genes expressed in at least 25% of the cells of either of the two compared populations, with a difference of at least 10% between the populations, a log fold change > 0.25. During log fold change calculation, the mean was computed as described above. Subsequently, cells were annotated in a two-step process based on cluster marker genes and protein marker expression as CD4 naïve (CCR7, LEF1, CD4, CD45RA, CCR7, CD62L), CD4 memory (AQP3, CD4, CD25), Activated CD4 memory (AQP3, CXCR3, ANXA5, CD4, CD28, HLA-DR), GZMH⁺ CD4 memory (GZMH, FGFBP2, CD4), Treg (FOXP3, IL2RA, IKZF2, CD4, CD25), CD8 naïve (CD8B, CCR7, LEF1, CD8, CD45RA, CCR7, CD62L), GMZK⁺ CD8 memory (CD8B, GZMK, CD8), GMZH⁺ CD8 memory (CD8B, GZMH, FGFBP2, CD8), MAIT (KLRB1, SLC4A10, CD8, CD161), gamma-delta (gd) (TRGC1, TRGC2, TRDC, CD3), CD4-CD8-double-negative (dn) (CD3, CD27, CD28), CD56^dim NK (KLRF1, IL2RB, CD56, CD16), CD56^bright NK (KLRF1, XCL2, NCAM1, CD56).

For differential expression analysis, cluster annotations of memory T cells were summarized into broader categories to ensure sufficient cell numbers per compared group and cluster ( ≥ 100 cells), as follows: CD4 memory (CD4 memory, Activated CD4 memory, GZMH⁺ CD4 memory), CD8 memory (GZMK⁺ CD8 memory, GZMH⁺CD8 memory).

Annotation and integration of other disease scRNA-seq datasets

The datasets from Schulte-Schrepping et al. (COVID-19, EGAS00001004571)³², Zhang et al. (Influenza, GSE243629)⁷⁸, and Yanagihara et al. (lung cancer, GSE215219)⁷⁹ were processed based on the steps used for processing the T and NK cell subset of our scRNA-seq EPTB dataset. Based on the respective dataset’s distribution, the lung cancer dataset’s cells were additionally filtered for > 250 genes and 10 % mitochondrial counts, while the Influenza dataset’s cells were filtered for > 250 genes, 20 % mitochondrial and 10 % hemoglobin counts. For both the lung cancer and Influenza datasets, the harmony⁸⁰ (v0.1.0) algorithm was applied in downstream processing to account for donor-specific batch effects. All datasets were subset for T and NK cells based on literature-known marker gene expression in an iterative clustering and filtering approach, as described before and using the following parameters: COVID-19 dataset: first 10 PCs, resolution 1.2; Influenza dataset: first 20 corrected PCs, resolution 0.4, then first 13 corrected PCs, resolution 0.7, then first 12 corrected PCs, resolution 0.7; lung cancer dataset: first 15 PCs, resolution 0.7, then first 13 corrected PCs, resolution 0.3. Final cell type annotation was performed separately for each dataset using Seurat reference mapping (FindTransferAnchors and MapQuery functions) from the EPTB dataset based on the first 20 PCs. Integration of all datasets was done for a combined UMAP visualization based on the first 10 PCs using reciprocal PCA. The final dataset contained 63,019 cells from the COVID-19 dataset (n = 13 controls, 8 mild COVID, 9 severe COVID patients), 22,670 cells from the Influenza dataset (n = 2 controls, 2 pregnant controls, 2 Influenza-infected children, 5 Influenza-infected adults), 64,559 cells from the lung cancer dataset (n = 4 patients) and the 29,662 cells from the EPTB dataset. Pregnant controls were not included in subsequent analyses. One control and a mild COVID donor each were excluded from cell state quantification due to a cell number

Correlation analysis of cell type frequencies

Correlation between measured (flow data and scRNA-seq) and estimated (bulk RNA-seq) cell type frequencies of patients present in all data modalities was calculated using Pearson’s correlation coefficient r.

Differential expression analysis in scRNA-seq data

Differential expression analysis was performed based on the Seurat FindMarkers with a Wilcoxon rank sum test. Differentially expressed genes were defined as genes with an abs(log2 fold change) > 0.25, were expressed in > 10 % of cells in the group of interest, and were below the 0.05 threshold for the bonferroni-adjusted p-values. Overrepresentation of DEGs in functional terms of the GO data base^52,53, KEGG^54,55,56 data base and Molecular Signature Database Hallmark gene sets ^57,58 was assessed using the clusterProfiler⁵⁹ (v4.0.5) enricher function. Terms with a bonferroni-adjusted p-value 

Module score calculation

Module scores were calculated using the R/Seurat AddModuleScore function at default settings for each respective signature. For the enrichment of the C1Q⁺ intermediate monocyte signature, the top 20 upregulated markers from Hillman et al.³⁵ were enriched per monocyte state. Statistics were computed using a Kruskal-Wallis test on donor level.

Pathway activity analysis

TNFα pathway activity was inferred using the PROGENy⁸¹ (v1.22.0) and decoupleR⁸² (v2.6.0) R packages in R (v4.3.0) and R Studio (v2023.3.0.386). Based on the PROGENy curated collection of pathways and target genes (n = 500 per pathway), a multivariate linear model was fit for each cell using the normalized gene expression data matrix as input. The model predicted gene expression based on pathway-gene interaction weights, obtaining t-values that were interpreted as pathway activity scores.

Luminex cytokine array

Luminex Discovery Assay (R&D Biotechne, Minneapolis, MN, USA) was performed according the manufacturer’s instructions. Samples were thawed on ice and centrifuged at 1000 × g for 5 min. Prior to analysis, samples were diluted 1:4 with Calibrator Dilutent (R&D Biotechne). All analytes were measured with Luminex 200 xMAP system (Thermo Fisher Scientific, Waltham, MA, USA) in technical duplicates. Standard curves were performed for all analytes, and xPONENT software (Thermo Fisher Scientific) was used for data collection and analysis. GraphPad Prism 9.5.1 (GraphPad) was used to analyse the data.

Flow cytometry and cluster analysis

PBMCs were thawed and washed (centrifugation at 400 x g for 5 min, for all washing steps) with DPBS. Next, PBMCs were blocked 20 min on ice with DPBS containing 10% FBS and 5 mM EDTA. PBMCs were washed with DPBS containing 1% FBS and 5 mM EDTA and stained with previously titrated antibodies (all antibodies are listed in Supplementary Table 3) in DPBS containing 1% FBS and 5 mM EDTA for 30 min at 4 °C. Afterwards, cells were washed and directly acquired using a S3 Symphony Cytometer (BD, Heidelberg, Germany). Cytometer settings were kept the same for all experiments, and data was analysed using FlowJo software (BD) (Version 10.8.1).

To verify identified cell types and visualize changes, unsupervised clustering was carried out. For this, the FlowSOM clustering algorithm⁸³ (version 4.0.0, available at https://www.flowjo.com/exchange/#/plugin/profile?id=7) was applied to the cells based on a chosen set of lineage-defining markers for each cell type. The cells were clustered into 100 different groups, which were visualized in a minimal spanning tree. From these, a biologically appropriate number of metaclusters for the respective cell types was generated. Metaclusters were named based on the included cell’s relative expression of defining surface markers. To circumvent bias caused by the variable number of patients within the different groups, the number of included cells was equalized before clustering using the FlowJo DownSample plugin (Version 3.3.1, available at https://www.flowjo.com/exchange/#/plugin/profile?id=25). In order to better visualize group-specific changes in individual clusters, the cell data was reduced to two dimensions using the t-distributed stochastic neighbor embedding (t-SNE) algorithm in FlowJo⁸⁴. The acquired FlowSOM metaclusters were then projected onto the respective t-SNE blots.

Statistical analysis

Statistical analysis was performed using GraphPad Prism 9.5.1 software (GraphPad) and the R package rstatix⁸⁵ (v0.7.0). For all experiments, full statistical parameters (value of n, standard deviation and statistical tests) are provided within the corresponding figure legends. For the whole study, p-values less than or equal to 0.05 were considered statistically significant and numerical p-values are shown within the figures. All data points represent biologically independent data points. Graphical illustrations were created with BioRender.com.

Ethics

The study is part of the Cologne EX-TB project (ClinicalTrials.gov Identifier: NCT06875336). All procedures involving human participants were conducted in accordance with the ethical standards of the institutional research committee, the 1964 Declaration of Helsinki and its later amendments, and the principles outlined in the Belmont Report. Written informed consent was obtained from all study participants prior to inclusion. Individuals gave consent that pseudonymised data can be provided in the manuscript. The study protocol, methodology, and all associated documents (e.g., informed consent form) were approved by the ethics committee of the University Hospital Cologne (identifier: 18-079).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.