Gene Set Variation Analysis (GSVA)

Pathway activity was estimated using GSVA with custom gene sets (33). GSVA is an unsupervised non-parametric alternative to ORA and GSEA as it does not depend on prior selection or ranking of the genes from group comparisons. The score produced by GSVA can be interpreted as the coordinated activation of genes from the gene set, summarizing the expression profile within individual samples. The gene sets used herein were compiled from literature and functional annotation databases. The ‘granulomatosis’ gene set was assembled with genes sourced from the Human Phenotype Ontology (HP:0002955), Gene Ontology (GO:0002432), and DisGeneNet (granulomatous diseases with score ≥ 0.1) totaling 30 genes ( Table S7 ). The macrophage polarization signatures are mainly from (34), where ‘M1’ (n genes = 25), ‘M2’ (n=20), ‘M2a’ (n=12), ‘M2b’ (n=9) and ‘M2c’ (n=12). ‘Pro-M1’ and ‘Pro-M2’ gene sets were built from literature (35, 36). Autophagy genes were retrieved from Reactome v.75 accessions: ‘Macroautophagy’ (R-HSA-1632852.8, n=137), ‘Chaperone Mediated Autophagy’ (R-HSA-9613829.3, n=22), and ‘Late endosomal microautophagy’ (R-HSA-9615710.3, n=34). Wilcoxon signed-rank test was used to test the differences between mean ranks between infected and mock macrophages signatures. Spearman’s rank correlation coefficient was calculated alongside 95% confidence intervals using DescTools R package v.0.99.40 (37). Principal component analysis (PCA) was applied to further explore the macroautophagy genes in the dataset. Mean-centered and variance standardized expression matrix with the 137 genes from macroautophagy gene set was subjected to PCA computed with FactoMineR v.2.4 (24). PCA scatter plot and contribution bar plots were graphed with factoextra v.1.0.6.