Assess Alzheimer’s Disease via Plasma Extracellular Vesicle-derived mRNA

Integrated transcriptome and pathway analysis

As shown in Figure 1A, we isolated and purified small EVs from 0.4 mL of plasma specimens, collected from a cohort of 82 individuals. The specimen information is summarized in Table S1. Following this, mRNA was extracted from these EV samples for sequencing library preparation. Subsequent RNA sequencing and gene alignment led to the identification of 19,495 genes. Categorization by biotype revealed the presence of 15,664 protein-coding genes (80.34%), 2,391 long non-coding RNAs (12.26%), and 1,440 genes of other biotypes (7.38%) within the identified dataset. This dataset, notably enriched with protein-coding genes, forms the foundation for the analysis conducted in our study. The volcano plots of the differential expression analysis within the MCI group (Figure 1B) revealed 1,097 significantly expressed protein-coding genes when compared to the Control group, including 897 upregulated and 200 downregulated genes. In the AD group (Figure 1C), 1,701 protein-coding genes were found to be significantly differentially expressed in comparison to the Control group, with 1,421 upregulated and 280 downregulated. We analyzed differentially expressed genes (DEGs) in each group to assess their involvement in biological pathways, utilizing g:Profiler against the KEGG database. The analysis indicated significant enrichment in the AD pathway. Specifically, in the MCI group, we observed an enrichment with a p-adjusted value of 1.36e-7, involving 27 AD-related genes, which have been identified in g:Profiler. In the AD group, the enrichment was more pronounced, with a p-adjusted value of 1.25e-13 and 44 AD-related genes. These genes are labeled in Figures 1B and 1C. Furthermore, we visualized the expression levels of these AD-related genes through heatmaps in Figures 1D and 1E. The heatmaps display distinct expression patterns in comparison to Control samples for both MCI and AD cases, offering insightful visual representation of the data.

• Download figure
• Open in new tab

Figure 1. Analysis of mRNA from plasma-derived EVs highlighting key AD-related genes in MCI and AD groups.

(A) Illustration of sample-preparation workflow. (B-C) Volcano plots showing differentially expressed protein-coding genes in MCI and AD groups compared to the Control group (p-adjusted < 0.05 and log2 fold change > |0.5|). The labeled genes are AD-related. (D-E) Heatmaps display expression patterns of the AD-related DEGs in MCI and AD groups compared to the Control group. Gene clustering in the heatmaps was based on hierarchical clustering using the Euclidean distance metric. (F-G) GO annotations for the AD-related DEGs in MCI and AD groups. (H-I) KEGG pathway enrichment analysis for AD-related DEGs in MCI and AD groups.

We conducted a comparison analysis of GO annotations for AD-related DEGs in the MCI (Figure 1F) and AD (Figure 1G) groups, focusing on the top 10 most significant annotations in each category. In the Biological Processes (BP) category, DEGs in MCI group predominantly engage with the Wnt signaling pathway, with 8 genes in both canonical and general Wnt pathways, and 10 genes in cellular response to oxygen-containing compounds. On the other hand, in the AD group, DEGs are mainly involved in energy metabolism pathways, such as oxidation of organic compounds (13 genes) and generation of metabolites and energy (14 genes). In the Cellular Components (CC) domain, MCI DEGs are primarily associated with protein complexes, including the catalytic complex (13 genes) and organelle membrane (14 genes). For AD group, the emphasis shifted to cellular structures like organelle membranes (24 genes) and catalytic complexes (23 genes). Within the Molecular Functions (MF) category, MCI genes were involved in binding activities, including frizzled binding (6 genes) and G protein-coupled receptor binding (6 genes). In contrast, AD genes focused on enzymatic functions, such as transmembrane transporter (11 genes), kinase (9 genes), and protein kinase activities (9 genes). This transition reflects the shift in molecular characteristics and functional emphasis from MCI to AD.

We subsequently conducted a KEGG pathway analysis of AD-related DEGs in both groups, revealing similarities and distinctions between the two conditions. The top 20 significant pathways for MCI (Figure 1H) and AD (Figure 1I) provide a detailed insight into the genetic landscape and underlying molecular mechanisms of both states. A key observation was the overlap in several pathways, including Alzheimer’s disease, neurodegeneration-related pathways, and other neurologically-linked diseases such as Spinocerebellar ataxia and Parkinson’s disease, common to both MCI and AD. Notably, the enrichment data indicated a stronger association with the AD group, especially in the Alzheimer’s disease pathway, where enrichment scores were significantly higher at 57.69 for AD group compared to 34.50 for MCI group. The difference suggests a progressive genetic intensification from MCI to AD state. While these shared genetic pathways revealed their interconnected nature, distinct pathways unique to each condition were also identified. MCI showed a tendency towards pathways typically linked to cellular processes and carcinomas, such as the mTOR signaling pathway, pluripotency of stem cells, and Melanogenesis. In contrast, AD was more aligned with pathways influencing metabolic and systemic processes, including Diabetic cardiomyopathy, Oxidative phosphorylation, and the Calcium signaling pathway. A further analysis of gene counts within these pathways revealed a consistent increase in gene numbers in AD group. For example, in the Huntington disease pathway shared by both conditions, AD group had a higher gene count with 22 genes compared to 12 genes in MCI group. This pattern indicates a potential accumulation of genes in AD, thereby amplifying its genetic complexity relative to MCI.

Identification of potential AD biomarkers

In the study to identify biomarkers for AD, we focused on a subset of DEGs shared by MCI and AD groups (Figure 2A), which made up 22.5% of the total AD-related gene pool. This indicates a notable genetic overlap between MCI and AD, suggesting these genes are potentially critical in AD progression. The expression patterns of these genes, analyzed through variance-stabilizing transformation (VST), showed distinct trends: 10 genes were upregulated and 3 were downregulated in MCI and AD stages compared to the Control group. These trends were quantified statistically as shown in Figure 2B. Generally, the upregulated genes showed higher median VST values in diseased states, while downregulated genes had lower median VST values, indicating a reduction in expression as the disease advances. An exception was PSMA2, which did not exhibit a clear up or downregulation trend when MCI and AD were compared with the Control group. The statistical significance of these differential expressions was marked with asterisks above the boxplots (Figure 2B).

• Download figure
• Open in new tab

Figure 2. Integrative analysis of 13 shared DEGs.

(A) Venn diagram showing the overlap of differentially expressed AD-related genes between the MCI and AD groups. (B) Box plot showing variations in gene expression across the 13 shared AD genes through different disease stages: Control, MCI, and AD. (*p < 0.05, **p < 0.005, ****p < 0.00005, ******p < 0.0000005) (C) Enrichment bubble plot illustrating the GO annotations for the 13 genes. (D) Enrichment bubble plot illustrating the KEGG pathways enriched by the 13 genes. (E) Summary of pathway enrichment listing the number of genes involved in each pathway. (F) Sankey diagram mapping the 13 genes to their associated biological processes. (G) Chord plot showing the relationships between the 13 genes and their linked KEGG pathways.

Functional enrichment and pathway analysis of common genes

We further conducted GO annotations for the 13 AD-related DEGs shared between MCI and AD groups (Figure 2C). The analysis ranked these annotations based on their enrichment scores relative to the background dataset. The significance of these findings was further confirmed by corresponding p-adjusted values. Within the Biological Process (BP) category, there was a significant enrichment in processes linked to cellular metabolism, including energy derivation through oxidation of organic compounds and regulation of precursor metabolites and energy generation, engaging 5 and 4 genes respectively. Notably, pathways related to cellular response to oxygen-containing compounds and response to these compounds involved 7 genes each, suggesting a potential connection to oxidative stress responses in neurodegenerative pathogenesis. In the Cellular Component (CC) domain, the protein-containing complex was prominent, with 11 genes involved, suggesting a major role of these molecules in AD-related cellular changes. The catalytic complex, enriched with 8 genes, emphasized the importance of enzymatic activities in the AD cellular environment. Additionally, the protein kinase complex and endoribonuclease complex were also highlighted, indicating a complex interaction of various protein complexes in the disease development. For Molecular Function (MF) annotations, insulin binding emerged as highly significant, involving 2 genes, indicating potential disruptions in insulin signaling pathways in AD. This is accompanied by nucleotide-binding functions such as ATP binding, adenyl ribonucleotide binding, and adenyl nucleotide binding, each implicating 4 genes. These findings suggest alterations in energy dynamics and signaling processes in AD.

Figure 2F shows a deeper insight into the distribution and roles of the 13 AD-related DEGs in Biological Process (BP), particularly in the pathways critical to cellular metabolism and the body’s response to oxidative stress. The Sankey diagram mapped the connection of these genes to specific biological processes, emphasizing both their gene ratio and p-adjusted significance. The analysis revealed a significant enrichment in the pathway for energy derivation through the oxidation of organic compounds, involving 5 genes: AKT2, PPIF, IDE, NDUFV1, and INSR. This pathway shows a notable p-adjusted value of 0.0013, confirming its significance. Additionally, the regulation of the generation of precursor metabolites and energy, which includes 4 genes (AKT2, PPIF, IDE, and INSR) with a p-adjusted value of 0.0017, further highlights the crucial role of these genes in key metabolic pathways that are often disrupted in AD.^{21, 22} Moreover, the pathways related to cellular responses to oxygen-containing compounds are notable as well, each involving 7 genes - AKT2, WNT2, PPIF, IDE, CHRM1, INSR, and ERN1. With a gene ratio of 0.538 and identical p-adjusted values of approximately 0.0021, these pathways emphasize the collective role of these genes in the oxidative stress response, which is crucial in the pathogenesis of neurodegenerative diseases like AD.^23–25

The KEGG pathway analysis of 13 AD-related DEGs is shown in Figure 2D, ranked by their enrichment scores. Notably, the Alzheimer’s disease pathway, with a significant enrichment score of 15.625 and a remarkably low p-adjusted value of 2.37e-16, emerged at the top of the list. This pathway, involving all 13 genes, emphasized the common genetic basis from MCI to AD. The detailed visualization of this pathway is illustrated in Figure S1, showing the involvement of these genes in a range of AD-associated biological processes. The ’Pathways of Neurodegeneration - Multiple Diseases’ follows with an enrichment score of 8.265, involving 10 genes, emphasizing their broader impact on various neurodegenerative diseases beyond Alzheimer’s. The Spinocerebellar ataxia and Huntington’s disease pathways also showed notable enrichment scores of 7.437 and 3.693, respectively, indicating genetic overlaps with these specific neurodegenerative disorders. Further, pathways like Diabetic Cardiomyopathy, Parkinson’s Disease, and Non-alcoholic Fatty Liver Disease, with enrichment scores of 3.257, 2.688, and 2.371, respectively, suggested the AD-related genes’ relevance to metabolic and systemic diseases, linking Alzheimer’s pathology with wider metabolic issues. Additionally, signaling pathways like the cGMP-PKG signaling pathway, key for vascular and neuronal functions, and the ’Longevity Regulating Pathway’, associated with aging, are also highlighted in the analysis with enrichment scores of 2.244 and 1.734. Even lesser-scored pathways, such as Amyotrophic Lateral Sclerosis and Prion Disease, demonstrated the implicated genes’ involvement in various disease mechanisms. Categorization of the pathways in Figure 2E indicated the significance of these genes in human diseases, particularly neurodegenerative disorders. The distribution of genes across these pathways, as reflected in the gene counts, aligns with the level of enrichment, painting a comprehensive picture of these genes’ roles in disease-related biological processes.

Figure 2G presents a chord diagram that maps the complex relationships between genes and pathways, providing a visualization of their regulatory impact. Notably, gene ATP2A2, interacting with various pathways including Diabetic cardiomyopathy and Spinocerebellar ataxia, indicated the possible wide-ranging effects of certain genes on both neurodegenerative and metabolic diseases. This suggests that the dysregulations observed in Alzheimer’s could have extensive systemic effects, potentially influencing or being influenced by peripheral tissue dysfunction.^{26, 27} NDUFV1, PSMA2, and SEM1 showed broad engagement with pathways such as Huntington disease, Amyotrophic lateral sclerosis, and Prion disease. Their upregulated expression in these diverse neurodegenerative conditions might suggest a common pathological mechanism or a unified response to neurodegenerative stress.²⁸ AKT2 and PPIF also exhibited wide-ranging associations with multiple pathways. AKT2’s involvement in these pathways reflected its dual role in the central nervous system’s pathology and peripheral metabolic disturbances related to AD.²⁶ PPIF’s widespread links suggested its role in mitochondrial dysfunction, a hallmark of neurodegenerative diseases.²⁹ The genes INSR, ATG13, and ERN1 display a diverse array of connections as well. INSR’s links to these pathways suggested that insulin resistance, a feature often seen in metabolic syndrome, may also play a role in neuropathology of AD.^{30, 31} On the other hand, the ties of ATG13 and ERN1 to longevity and mitochondrial pathways emphasized the significance of cellular aging and energy regulation in the progression of AD.^{32, 33} These correlations are enlightening as they reveal the intersection of metabolic regulation, aging, and neurodegeneration, suggesting that genes involved in AD could be crucial in these broader biological processes.

Biomarker selection through correlation analysis and clustering

We conducted a refinement of the 13 AD-related DEGs through biweight midcorrelation analysis to characterize the similarity in expression among the biomarker candidates. This approach identified gene pairs with a high degree of correlation, specifically those with a correlation coefficient above |0.65|. For instance, ATP2A2 and SEM1, as shown in Figure 3A, were closely aligned, suggesting potential exclusion to maintain dataset integrity and avoid redundancy. Additionally, the dataset was assessed through a comprehensive clustering analysis using Euclidean distance, resulting in a heatmap shown in Figure 3B. This heatmap clustered genes based on their expression patterns across all diagnostic categories: AD, MCI, and Control. A dendrogram on the heatmap shows the hierarchical relationships between genes, aiding in identifying closely related pairs. To minimize dataset redundancy, gene pairs that clustered closely, such as ATG13 with INSR, PSMA2 with HSD17B10, NDUFV1 with PPIF, WNT2 with CHRM1, and SEM1 with ATP2A2, were marked for potential exclusion. The goal was to select one gene from each pair for the final biomarker panel, ensuring a unique and non-overlapping set for further validation.

• Download figure
• Open in new tab

Figure 3. Biomarker selection through correlation analysis and clustering.

(A) Biweight correlation matrix showing genes with strong correlations. Color intensity represents the correlation strength with darker hues indicating stronger correlations. (B) Heatmap with a dendrogram constructed using Euclidean distance, showing gene expression similarities among the 13 AD genes. (C-E) Scatter plots showing the median expression distributions of the Control and AD/MCI groups across three selected gene sets.

Following the analysis, we derived three potential biomarker panels for further examination. As illustrated in Figures 3C-E, each panel is characterized with a scatter plot showing the median expression values for the two diagnostic categories, the Control group and a combined group of MCI and AD. Panel 1 (Figure 3C), the unrefined group, includes all the 13 genes. Panel 2 (Figure 3D), consisting of 8 genes, was created by removing one gene from each redundant pair identified in the clustering assessment. An alternative to panel 2 (Figure S2A), also containing 8 genes, was formed by choosing the alternate gene from each pair that was omitted in panel 2, serving as a comparison. The last panel (Figure 3E) was compiled by selecting one gene from the dendrogram branch of downregulated genes and three from the initial splits of upregulated genes, ensuring diversity of gene expression patterns in the biomarker panel. The scatter plots showed that the selected biomarker panels maintained a similar range of gene expression, allowing comprehensive coverage of potential biomarkers while minimizing redundancy, thereby enhancing their overall applicability in the subsequent validation stage.

Development and characterization of AD prediction model

For each biomarker panel, we developed an AD prediction model by training a linear Support Vector Machine (SVM) classifier with the dataset. A nested cross-validation (5-fold inner and outer cross-validation) was first employed to determine the optimal hyperparameters for the SVM model. Subsequently, both the AD/MCI and control groups were randomly divided into a training set (80%) and a testing set (20%) each time for 500 trials of model training and testing, in order to characterize model performance in terms of Receiver Operating Characteristic (ROC) Area Under Curve (AUC) value. Figure 4A-C illustrate the mean ROC curves for each biomarker panel.

• Download figure
• Open in new tab

Figure 4. Model characterization for AD evaluation.

(A-C) Training and testing of an SVM classification model for each biomarker panel with the performance characterized through ROC curves, highlighting the AUC and 95% confidence interval. Inset figures present the distribution of AUC values over 500 iterations, illustrating the model’s performance consistency. (D-F) Comparison of AD scores for each sample as generated by the prediction models, differentiating between AD/MCI and control groups. The significance of differences was characterized using the Wilcoxon test. (G-I) Evaluation of correlation between AD scores and CDR values. Statistical differences among groups were assessed using the Kruskal-Wallis test.

The shading of the curves indicates the standard deviations across the 500 training and testing iterations, reflecting the consistency of the models’ performance. Additionally, the distribution of AUC values (inset of Figure 4A-C) provides insights into the variation of the models’ prediction outcomes over these iterations.

The first biomarker panel (Figure 4A), incorporating all the 13 genes, displayed outstanding training performance with an AUC of 1.00. This high level of accuracy continued in the testing phase, where it obtained a mean AUC of 0.9983, with a tight 95% confidence interval (CI) between 0.9978 and 0.9989, indicating strong predictive reliability. The second biomarker panel (Figure 4B), a refined cohort consisting of 8 biomarkers, recorded a mean training AUC of 0.9992, slightly lower than that of the 13-gene set. In the testing, it reached a mean AUC of 0.9864, with a 95% CI stretching from 0.9845 to 0.9882, suggesting a slight reduction in classification capability compared to the first panel, but it still retained high diagnostic precision. An alternative version of the second panel (Figure S2B), also with 8 biomarkers, had a comparable but lower AUC for both training and testing to the original second biomarker panel. The final biomarker panel (Figure 4C), which included only 4 genes, exhibited a training AUC of 0.9611. Although lower compared to the other panels, it still maintained a strong capacity to differentiate between conditions. For this biomarker panel, the testing phase yielded a mean AUC of 0.9397, accompanied by a wider 95% CI ranging from 0.9349 to 0.9444, reflecting greater variability in the model’s predictive accuracy.

The AD classifier for each marker panel was developed by averaging the weights and bias terms from the linear SVM model across 500 iterations, as detailed in Supplementary Information. This process yielded a final model that generated an AD score for each sample. These scores were then compared between two groups (AD/MCI and Control), as illustrated in Figures 4D-F and Figure S2C. Generally, it was observed that as the number of genes in the biomarker panels reduced from 13 to 4, the overlap in the interquartile ranges of the AD/MCI and Control groups became more significant. This trend was confirmed by an increase in the p-value from the Wilcoxon test, indicating less distinction between the groups. However, the model retained significant discriminative ability even with only 4 gene markers (Figure 4F), as evidenced by a statistically significant p-value of 3.2e-11.

The prediction model not only demonstrated the ability to discriminate between the Control and AD/MCI groups but also exhibited a strong correlation with the clinical dementia rating (CDR). We subdivided the samples into various groups based on their CDR values to compare their AD scores, as shown in Figures 4G-I and Figure S2D. As anticipated, there was a significant increase in the median AD scores, progressing from the control group (CDR 0) to the AD/MCI group (CDR 0.5-3). More importantly, the median scores showed a gradual decrease as the CDR value increased from 0.5 to 3. This trend was consistently observed across the three biomarker panels, indicating the model’s ability to differentiate among AD stages.

Clinical specimens

We analyzed a total of 82 plasma specimens, including 25 samples from healthy controls, 13 from individuals with MCI, and 44 from patients diagnosed with AD. Demographics and clinicopathological characteristics of these participants are provided in Table S1. These specimens were from three independent patient cohorts: Washington University in St. Louis (Knight ADRC), Indiana University (NCRAD), and PrecisionMed. The study was approved by the institutional review boards of all participating institutions with written informed consent provided by all participants involved in the study.

EV isolation, library preparation, and RNA sequencing

For isolating EVs, 350 μL of each plasma sample was first diluted with 1×PBS to a final volume of 15 mL. The diluted plasma was then passed through a 0.22 μm syringe filter to remove large particles and aggregates. The isolation of EVs was conducted using the EXODUS platform, which utilized 25 mm-diameter exosome isolation devices (EID) according to the previously published method.⁴⁶ For comparison, the exoEasy Maxi Kit (Qiagen) was employed. Once isolated, the EVs were resuspended in 400 μL of PBS. From this suspension, a 100 μL aliquot of sample was further diluted 4 times for nanoparticle tracking analysis (NTA) using NanoSight (Malvern Panalytical, LM14). The concentration and size distribution of EVs are shown in Figure S3.

Another 300 μL of EVs in PBS were then mixed with an equal volume of Dynabeads M-270 Streptavidin beads (ThermoFisher Scientific, 65305), previously conjugated with a biotinylated oligonucleotide that included a PCR handle and a Poly(dT) sequence. These beads were suspended in a lysis buffer containing 0.4% Sarkosyl to facilitate the capture of mRNA from lysed EVs. Following this, solid-phase reverse transcription was performed using Maxima H Minus Reverse Transcriptase (ThermoFisher Scientific, EP0752). The cDNA synthesized on the beads underwent amplification through SMART PCR, using Kapa HiFi Hotstart Readymix (Roche, KK2601), and the resulting DNA was subsequently purified using AMPure XP beads (Beckman Coulter, A63881). The purified DNA was further subject to a second round of PCR and purification. The PCR products then underwent tagmentation using a Nextera XT Library Preparation Kit (Illumina, FC-131-1096). For fragmented mRNA within EVs, NEBNext Small RNA Library Prep Set (New England BioLabs, E7330L) was employed for library preparation according to manufacturer’s instructions. Quality control of the NGS library preparation was accessed via a TapeStation 4200 with High Sensitivity D5000 ScreenTape (Agilent). Finally, RNA sequencing was performed on an Illumina NovaSeq 6000 platform with at least 20M of reads for each library.

Genome alignment and gene annotation

The sequencing data were first trimmed using Trimmomatic (version 0.4)⁴⁷. Subsequently, these trimmed reads were aligned to the Ensembl reference genome for Homo sapiens (GRCh38.109) using the STAR aligner (version 2.7.10a)⁴⁸. Duplicate reads were identified and eliminated using Picard tools (version 3.0). The raw reads were extracted using the gene transfer format from the Ensembl database corresponding to the Homo sapiens GRCh38.109 genome. This process was facilitated by several R libraries, including GenomicFeatures (version 1.46.5)⁴⁹, GenomicAlignments (version 1.30.0)⁴⁹, and Rsamtools (version 2.10.0)⁵⁰. The reads obtained were then annotated with external gene names and biotypes using the biomaRt package (version 2.50.3)^{51, 52}. Low-read counts were filtered out, retaining only those present in at least 20% of individual samples within each of the study groups (i.e., Control, MCI, and AD), with a minimum of 2 read counts per sample. Moreover, any genes that lacked annotations of external gene names were excluded from the analysis.

Statistical analysis

Differential gene expression analysis was conducted using DESeq2 (version 1.40.2)⁵³ supplemented by the sva package (version 3.48.0)⁵⁴ to account for potential hidden batch effect. Criteria for significance were established by requiring a log2 fold chance > |0.5| and adjusted p-value (p-adjusted) < 0.05. The preprocessed data were normalized using DESeq2, followed by a filtering process to retain only those entries where at least 2 samples exhibited normalized values >10 before passing through the sva package and subsequent processes. Following the differential expression analysis, data subjected to variance stabilizing transformation (VST) were retrieved from DESeq2, processed further with the sva and limma packages (version 3.56.2)⁵⁵ for downstream analyses per recommendation outlined in the DESeq2’s vignette.

Functional profiling was then performed using g:Profiler (version e110_eg57_p18_4b54a898)⁵⁶, employing the g:SCS significant threshold with p-adjusted < 0.05. In identifying KEGG AD-related genes from the pool of differentially expressed genes, the analysis utilized statistical domain scope specified as ‘All known genes.’ In contrast, for more focused profiling of these determined AD-related genes, the statistical domain scope was narrowed to ‘Only annotated genes.’

The linear SVM model, implemented in the e1071 package (version 1.7-13)⁵⁷, was employed to develop a diagnostic classifier. This model was paired with a nested cross-validation strategy with inner and outer folds set to 5, facilitated by the caret package (version 6.0-94)⁵⁸, to address class imbalances within the dataset. The classifier’s performance was assessed through ROC analysis, utilizing the pROC package⁵⁹.

For correlation analysis, the biweight midcorrelation from the WGCNA package ^{60, 61} was utilized. Clustering on the heatmap was performed using the ‘complete’ method with Euclidean distance as the measure. Statistical differences in decision values between AD/MCI and Control groups were determined using the Wilcoxon test. Furthermore, the Kruskal-Wallis test was applied to compare decision values across various CDR groups from CDR 0 to 3.

Visualizations

Volcano plots were produced with GraphBio⁶². Boxplots were constructed using BoxPlotR⁶³. The mean ROC curves were generated by RStudio. For the other plots, SRPlot, a web tool for data analysis and visualization, was employed.