Experimental and computational workflow.

Data resources for the project included peripheral blood (PB) and peripheral blood leukocyte (PBL) RNA-seq data from four studies; Two studies from (31) (OGR25-BTB) and (35) (MCL14-BTB) investigating transcriptional differences between non-infected animals (bTB−) and cattle naturally infected with M. bovis (bTB+), and two studies consisting of bTB− and bTB+ cattle experimentally infected with M. bovis and longitudinally sampled across an experimental time course from (36) (MCL21-BTB) and (38) (WIA20-BTB). Data analysis procedures included: (1) an intra-dataset differential expression (DE) and functional enrichment analysis; (2) development and tuning of machine-learning models in the training set via five-fold cross-validation using differentially expressed genes identified in the training set (see “Methods”); (3) evaluation of hyperparameter-tuned models in the testing set and assessing performance by calculating the area under the receiver operating characteristic curve (AUROC), sensitivity and specificity, respectively; and (4) evaluation of hyperparameter-tuned models in non-bTB datasets comprised of non-infected animals and cattle infected with M. avium spp. paratuberculosis (MAP), bovine herpes virus (BoHV−1), or bovine respiratory syncytial virus (BRSV), respectively.

Overview of datasets used in this manuscript.

Intra-dataset differential expression and functional enrichment analysis.

(A) Upset plot illustrating the number of significantly (Padj. < 0.05) differentially expressed genes (DEGs) identified in the four bovine tube culosis (bTB) datasets and the overlap of these DEGs among the studies. (B) Bi-directional jitter plots for significantly impact pathways (SIPs) perturbed by DEGs exhibiting increased and separately, decreased expression in each dataset, respectively for four databases: (1) gene ontology (GO) biological processes (GO:BP); (2) GO cellular component (GO:CC); (3) Kyoto Encyclopaedia of Genes and Genomes (KEGG); and (4) Reactome (REAC). Dashed dotted lines indicate the −log10Padj. threshold of 0.05 for characterising a SIP. The total number of input DEGs and identified SIPs are also detailed for each analysis of increased and decreased DEGs in each dataset, respectively. The jitter plots are arranged from left to right according to their documented level of innate immune response pathway activation. (C) Heatmap showing the log2 fold-change (LFC) in effect size estimates between control non-infected cattle and animals infected with M. bovis for 50 significant DEGs between bTB− and bTB+ cattle identified in all naturally infected datasets and in at least one time point in each of the time course datasets, respectively. Red colours indicate increased expression, and blue colours indicate decreased expression in the bTB+ group relative to the bTB− group. The * sign indicates that the LFC value was significantly greater than or less than 0 in the corresponding comparison between bTB+ and bTB− cattle.

Machine learning (ML) model training and evaluation.

(A) Volcano plot illustrating significantly differentially expressed genes (DEGs) identified in the training set for the bTB+ (n = 93) vs bTB− (n = 69) contrast with thresholds determined by FDR-Padj. < 0.05 and an absolute log2 fold-change (LFC) > 0. Genes are coloured to indicate increased (red) or decreased (blue) expression in the bTB+ group relative to the bTB− group, respectively, and triangular datapoints indicate the genes retained for subsequent ML analysis with the criterion of mean median-of-ratios normalised expression count >100 in the training dataset. (B) Bi-directional jitter plot for significantly impacted pathways (SIPs) perturbed by DEGs exhibiting increased and separately, decreased expression in the training set for four databases: (1) gene ontology (GO) biological processes (GO:BP); (2) GO cellular component (GO:CC); (3) Kyoto Encyclopaedia of Genes and Genomes (KEGG); and (4) Reactome (REAC). Dashed dotted lines indicate the −log10Padj. threshold of 0.05 for characterising a SIP. The total number of input DEGs and identified SIPs are also detailed for each analysis of increased and decreased genes in each dataset, respectively. (C) Area under the receiver operating characteristic curve (AUROC) values for eight hyperparameter-tuned ML models in each of the five folds of the training set. (D) An AUROC depicting the performance of each of the hyperparameter-tuned models in the testing set. (E) Boxplot showing the distribution of AUROC estimates of animals in the testing dataset depending on whether the dataset comprised of animals experimentally or naturally infected with M. bovis. (F) Heatmap showing the predicted class of animals in the testing set when selecting an optimum threshold to maximise the sensitivity (> 0.90), where possible, of each classifier.

Greedy-forward search strategy evaluation in training and testing sets.

(A) Line plot showing the average area under the receiver operating characteristic curve (AUROC) across the five folds in the training set for each combination of genes identified in the first pass or second pass of the greedy-forward search algorithm, respectively. (B) Line plot showing the 95% confidence interval of AUROC estimates for the 13-gene set, 17-gene set and combined 30-gene set in each of the five cross-validation folds in the training set. (C) Boxplots showing the distribution of the scaled infection z-score derived from each of the three gene sets for all bTB− (blue; circle) and bTB+ (red; triangle) animals in the training set, separated by study of origin. (D) Boxplots showing the distribution of the scaled infection z-score derived from each of the three gene sets for all bTB− (blue) and bTB+ (red) animals in the testing set, separated by study of origin. (E) An AUROC depicting the performance of each of the three gene sets in the testing set. (F) Heatmap showing the predicted class of animals in the testing set when selecting an optimum threshold to maximise the sensitivity (> 0.85) for each of the three gene sets.

Evaluation of machine learning and greedy forward search strategies in external datasets.

(A) An area under the receiver operating characteristic curve (AUROC) showing the ability of the hyperparameter-tuned models to discriminate between control non-infected animals and cattle infected with M. avium spp. paratuberculosis (MAP) from (39) (ALO19-MAP) dataset. (B) An AUROC showing the ability of the hyperparameter-tuned models to discriminate between control non-infected animals and cattle infected with bovine herpes virus (BoHV-1) from (41) (ODO23-BRD) dataset. (C) An AUROC showing the ability of the hyperparameter-tuned models to discriminate between control non-infected animals and cattle infected with bovine respiratory syncytial virus (BRSV) from (40) (JOH21-BRD) dataset. (D) Boxplots showing the distribution of the scaled infection z-score derived from each of the three gene sets for all control non-infected animals and cattle infected with either MAP, BoHV-1, or BRSV, respectively. (E) Three AUROCs illustrating the performance of the three gene sets in terms of discriminating between control non-infected animals and cattle infected with either MAP, BoHV−1, or BRSV, respectively. (F) An AUROC curve showing the discriminatory power of all hyperparameter-tuned ML models and gene sets when comparing bTB+ cattle to MAP+, BoHV-1+ and BRSV+ cattle, respectively.