Genome-based features are increasingly used for viral identification and classification, yet the extent to which viral taxa are intrinsically diagnosable from genomic sequence has not been systematically evaluated across taxonomic ranks. Here, we quantify the diagnosability of viral taxa at family, genus, and species levels using interpretable diagnostic rules derived from k-mer sequence features. Kmer presence/absence (PA) characters were determined using a kmer sweep approach. Diagnostic rule sets were discovered using the PA characters and evaluated using precision, recall, confusion matrices, and receiver operating characteristic (ROC) analysis. At the family level, diagnosability across ~360 viral families is limited, reflecting extensive sequence heterogeneity and weakly conserved diagnostic features at deep evolutionary scales. Partitioning families by Baltimore classification reduces this heterogeneity but does not substantially increase the number of diagnosable families. This result likely reflects both the extensive evolutionary divergence among viral families and the fact that Baltimore groups are defined by genome type and replication strategy rather than shared evolutionary history, which may limit the extent to which such partitions capture consistent diagnostic signal. In contrast, diagnosability increases substantially at finer taxonomic levels. At the genus and species levels, diagnostic rules are readily identified and achieve strong classification performance metrics (mean AUC ≈ 0.98 for diagnosable taxa). However, single k-mer rules are often insufficient for optimal discrimination. Instead, higher-order rule structures—particularly paired and composite k-mer diagnostics—substantially improve performance, indicating that combinatorial sequence features are required to resolve closely related taxa. These results demonstrate that diagnostic signal becomes more accessible at finer taxonomic resolution but requires increasing rule complexity to fully capture taxonomic differentiation. Importantly, we show that diagnosability and classification are distinct properties: diagnostic features may exist even when accurate classification is not achievable under competitive frameworks. This study provides a quantitative and interpretable framework for evaluating genome-based viral diagnostics and shows that diagnostic rule complexity scales predictably with taxonomic resolution.

Materials and Methods

Overview 

For this archive we have grouped the scripts and readme files into the following larger directories:

• Family_script_final_2.zip: holds the scripts described below and in the document README_Family.docx
• Family_Script_Library.docx: An archival list of the scripts used in this study
• README_Family.docx: A fully documented README document that describes all of the scripts in this archive.

For each analytical step, we describe first the goal of the method, second the methodological approach, and third the scripts used to implement the analysis. All scripts are available via DataDryad using the filenames listed below.

Viral genome datasets and taxonomic scope: The goal of this step was to define the taxonomic scope of the analysis and identify viral families suitable for evaluating genome-based diagnosability. Because diagnosability depends on the presence of consistent sequence features within taxa, analyses were restricted to families with sufficient genomic sampling to support independent rule discovery and validation.

For family-level analyses, families were retained if at least four complete genomes were available, allowing partitioning into minimal training (n ≥ 2) and test (n ≥ 2) sets. These rules decreased the number of viral families available for examination from >350 to about 250.  This reduction means that over 100 families have too few representative genomes to be useful in a study of diagnostic character utility.

For genus- and species-level analyses, additional criteria were applied to ensure robust representation of taxonomic structure. First, families were retained at this level only if the number of available genomes exceeded the number of recognized species by approximately fourfold, ensuring adequate sampling within taxa for both rule discovery and evaluation. Second, families were required to contain at least two genera, enabling assessment of diagnosability within a hierarchical taxonomic framework.

Together, these criteria define minimal conditions for independent validation and meaningful taxonomic comparison within a family. Application of these filters reduced the dataset from approximately 360 viral families to 79 families suitable for downstream analyses. Although these restrictions represent a substantial reduction in taxonomic scope, the retained families span a wide range of genome architectures, sizes, host associations, and evolutionary divergence levels (Supplemental Table 1). Genome sizes range from small circular genomes (~2 kb) to large double-stranded DNA genomes (>100 kb), providing a representative cross-section of viral genomic diversity.

Candidate families and genome sampling statistics were identified using summary tables derived from NCBI viral genome downloads and associated taxonomy. For each family-level FASTA file, a corresponding “truth” table was generated containing accession numbers and standardized taxonomic assignments based on ICTV classifications. The following scripts were used in this step: taxonomy_of_taxonomy.sh; qc_taxonomy_of_taxonomy.sh; build_family_master_table.py.  These scripts generated summary tables describing genome counts, taxonomic structure, and sampling ratios for each family and are available as all scripts are on DataDryad.

Genome acquisition and dataset curation: The goal of this step was to obtain complete viral genome sequences and construct curated taxonomic reference tables linking each genome accession to family, genus, and species assignments. Genome sequences were downloaded from NCBI nucleotide and viral genome databases using family-level queries restricted to complete genomes. For each accession, associated taxonomic metadata were extracted and mapped to three hierarchical levels (family, genus, species). Because public viral genome annotations frequently contain inconsistencies, a series of curation steps was applied to ensure taxonomic coherence. These included removal of incomplete records, normalization of taxonomic strings, standardization of genus- and species-level naming conventions, exclusion of ambiguous or placeholder taxa (e.g., “sp.”, “uncultured”), correction or removal of inconsistent taxonomic assignments, and filtering of sparsely sampled taxa that could not support downstream analyses. These curated truth tables serve as authoritative mappings between genome accessions and taxonomic labels and form the foundation for all subsequent analyses. Ensuring taxonomic consistency at this stage is critical for evaluating diagnosability, as errors or ambiguities in taxonomic labels would directly confound assessment of diagnostic rule performance.** Genome acquisition and curation were performed using the following scripts: download_family_genomes.sh; build_family_truth_table.sh; build_truth_batch.sh; qc_taxonomy_of_taxonomy.sh; normalize_taxonomy_strings.py; build_family_kmer_cache.py; discover_ladder_rules.py; aggregate_family_results.py; download_from_accession_list.py.  Outputs from this step included FASTA genome files and curated taxonomic truth tables.

Train–test dataset partitioning: The goal of this step was to generate independent datasets for diagnostic rule discovery and evaluation. For each viral family, genomes were partitioned into training (rule discovery) and test (evaluation) sets using stratified sampling to preserve taxonomic representation across hierarchical levels (family, genus, species). Within each family, 50–70% of genomes were assigned to the training set and 30–50% to the test set. This partitioning ensures strict separation between rule discovery and evaluation, preventing overfitting and enabling unbiased assessment of diagnostic performance. To maintain computational tractability while preserving taxonomic diversity, a maximum of 10,000 genomes per family was imposed. For families exceeding this limit, stratified subsampling was applied to retain representative diversity across taxa without introducing sampling bias.** Partitioning was implemented using the following scripts: select_family_genomes.py; stratified_subsample_genomes.py; build_train_test_split.py.  Outputs included selected genome tables and train–test split tables defining the composition of each dataset.

Diagnostic rule types and rule hierarchy: The goal of this step was to define diagnostic sequence features and their hierarchical organization. Diagnostic rules were constructed from k-mer sequence features, where a k-mer is a contiguous nucleotide substring of length k. For each genome, k-mer presence/absence was recorded and compared across taxa to identify features that uniquely diagnose taxonomic groups.

Three classes of diagnostic rules were defined:

1. Single k-mer rules (S): A k-mer present in all genomes of a target taxon and absent from all other taxa.

2. Paired k-mer rules (P): Pairs of k-mers whose joint presence/absence uniquely identifies a taxon when single k-mers are insufficient.

3. Composite rules (C): Higher-order combinations of k-mers capturing conditional or multi-feature relationships required to diagnose taxa not resolvable by single or paired rules.

These rule types form a hierarchical “rule ladder” of increasing complexity (S → P → C), in which progressively more complex combinations of sequence features are used to achieve diagnosability (Figure 1). Importantly, because k-mer features are represented as binary presence/absence characters, classical sequence diagnostic categories based on multi-state nucleotide patterns (e.g., single-pattern or paired-pattern diagnostics) are not directly represented in this framework. Instead, such patterns are implicitly captured through combinations of binary features within the S, P, and C rule hierarchy. Rule discovery was implemented using the following scripts: scripts/build_family_kmer_cache.py; scripts/discover_ladder_rules.py 

Outputs included rule tables for each taxon and rule type.

Figure 1. Discovery of diagnostic k-mer rule sets and analysis ladder. Datasets were partitioned into a rule discovery dataset (70%; red) and an independent test dataset (30%; green). K-mers were generated from genomes in the rule dataset, and a presence/absence matrix was constructed for all taxa.(1) K-mer generation produces libraries of sequence features of length k.  (2) Presence/absence matrices are used to identify diagnostic rules for each taxon. Single k-mers may uniquely diagnose taxa (S rules), while combinations of k-mers may be required when single features are insufficient (P rules). More complex taxa require composite combinations of multiple k-mers (C rules). (3) Discovered rules are applied to genomes in the independent test dataset.  (4) Predictions are summarized in a confusion matrix.  (5) Diagnostic performance is evaluated using standard classification metrics (accuracy, recall, precision, F1 score, specificity).The rule ladder (S → P → C) represents increasing levels of combinatorial complexity required to achieve taxonomic diagnosability.

Diagnostic performance metrics and diagnosability: The goal of this step was to quantify diagnostic performance and to distinguish between the existence of diagnostic signal (diagnosability) and classification accuracy. In this framework, diagnosability refers to the ability of sequence features to produce unambiguous, taxon-specific rules, whereas accuracy reflects the correctness of taxonomic assignments when such rules are applied. Rule sets derived from training data were applied to independent test datasets. Classification outcomes were summarized using confusion matrices, from which standard performance metrics were calculated, including precision, recall (sensitivity), and F1 score. In addition, two metrics were used to capture properties specific to diagnostic systems. Coverage was defined as the proportion of genomes receiving at least one diagnostic assignment, reflecting the extent to which taxa are diagnosable under the rule framework. Ambiguity rate was defined as the proportion of genomes assigned to multiple taxa, indicating cases in which diagnostic signal is present but not uniquely resolving. Together these metrics distinguish between (i) the presence of diagnostic signal, (ii) the ability to assign genomes, and (iii) the accuracy of those assignments, allowing diagnosability to be evaluated independently of classification performance.** Scripts used in this step include: score_family_diag_rules.py; eval_family_diag_results.py 

Diagnostic scoring and threshold calibration: To enable quantitative evaluation of rule-based diagnostics, a scoring framework was implemented to summarize rule support for each genome–taxon assignment. For each genome, the presence of diagnostic rules was evaluated across all candidate taxa, and a composite score was computed based on the number and strength of supporting rules. Scores were normalized to allow comparison across taxa and rule types. In cases where multiple taxa received support, the taxon with the highest score was selected as the predicted assignment; ties were treated as ambiguous classifications. To assess the relationship between score and classification reliability, thresholds were applied to exclude low-confidence predictions. By varying the score threshold, a trade-off between coverage (proportion of genomes classified) and accuracy was evaluated. This threshold-based analysis enables identification of high-confidence diagnostic assignments and provides a quantitative framework for calibrating rule-based classification.  Scripts used in this phase of the work are: score_family_diag_rules.py; eval_family_diag_results.py; sweep_thresholds.py; plot_score_vs_correctness.py

ROC and pseudo-ROC analysis: The goal of this step was to evaluate the separability of taxa under different rule types and levels of rule complexity. Receiver operating characteristic (ROC) curves were used to assess classification performance by examining the trade-off between true positive and false positive rates.In rule-based systems, continuous classification scores are not inherently produced. To approximate ROC behavior, thresholds were applied to rule support, defined as the proportion of training genomes within a taxon satisfying a given rule. By varying this threshold, true positive and false positive rates were calculated, and the area under the curve (AUC) was used as a summary measure of classification performance.In addition, pseudo-ROC curves were constructed to evaluate the accumulation of diagnostic signal. In this approach, diagnostic rules were ranked by support, and classification performance was evaluated as rules were incrementally added. This procedure provides an analog to feature-ranking approaches in standard classification frameworks, allowing assessment of how rapidly diagnostic power increases as additional rules are incorporated. The resulting pseudo-AUC quantifies the rate at which diagnostic signal accumulates, with AUC ≈ 1.0 indicating that a small number of high-support rules are sufficient to achieve strong classification performance. ROC and pseudo-ROC analyses were performed across rule types and taxonomic levels to compare the diagnostic properties of different feature classes. Scripts used in this step include:  plot_family_diag_metrics.py; generate_heatmaps.py. 

Comparative analysis across viral families: The goal of this step was to compare diagnosability across viral families (and genera and species). Performance metrics were aggregated across families, and rule sets were evaluated across different rule types and parameter settings. For each family, the best-performing rule set was identified based on AUC, evaluated in conjunction with coverage and ambiguity to ensure that strong classification performance reflected robust and usable diagnostic signal. Summary tables were generated containing AUC, precision, recall, F1 score, and coverage metrics for each family.Heatmaps were constructed to visualize variation in diagnostic performance across families and taxonomic levels. These visualizations highlight differences in both classification accuracy and diagnosability. Families lacking reliable diagnostic signal were classified as “unsolved,” defined as those failing to achieve adequate coverage and/or producing ambiguous or low-confidence assignments under the rule framework. These families represent cases in which genomic features do not provide sufficient discriminatory power for robust taxonomic diagnosis. Scripts used in this step include: aggregate_family_results.py; generate_heatmaps.py; plot_family_comparison.py. 

Benchmark evaluation on ICTV test dataset: To evaluate performance on a large and heterogeneous dataset, diagnostic rule sets were applied to an independent benchmark dataset derived from ICTV-associated viral genomes. Predicted taxonomic assignments were compared to curated truth labels using accession-based matching. Performance was evaluated using standard classification metrics, including accuracy, precision, recall, and F1 score. In addition, coverage and ambiguity were assessed to quantify the proportion of genomes receiving high-confidence assignments. This benchmark analysis provides an external validation of diagnostic performance across a broad spectrum of viral diversity. Scripts used in this step include: eval_family_diag_results.py; aggregate_family_results.py; plot_score_vs_correctness.py.

Computational implementation: All analyses were implemented using Python scripts and Unix shell workflows organized into a modular pipeline. The pipeline comprises genome preprocessing, k-mer generation, rule discovery, rule scoring, performance evaluation, and cross-family aggregation. Analyses were conducted for k-mer lengths ranging from k = 8 to k = 15 to evaluate the effect of feature resolution on diagnostic performance. Detailed analyses focused on k = 9–12, with k = 11 selected for downstream analyses based on optimal performance across ROC/AUC metrics, coverage, and classification accuracy. Scripts used for these summaries are:  build_family_kmer_cache.py; discover_ladder_rules.py; aggregate_family_results.py; kmer_sweep_analysis.py. 

 Output and reproducibility: All outputs, including rule tables, classification results, ROC and pseudo-ROC statistics, and summary performance metrics, were stored as tab-delimited files and graphical figures. Intermediate datasets (e.g., train–test partitions, k-mer presence/absence matrices, and rule sets) were retained to ensure full traceability of all analytical steps. This organization enables complete reproducibility of the pipeline from raw genome acquisition through final evaluation. All scripts (Supplemental Table 2) and associated data products are available via DataDryad, allowing independent replication and extension of the analyses.