aardvark-bio 0.10.5

Aardvark - A tool for sniffing out the differences in vari-Ants
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247


use anyhow::ensure;
use clap::Args;
use log::info;
use serde::Serialize;
use std::path::PathBuf;

use crate::cli::core::{check_optional_filename, check_required_filename, AFTER_HELP, FULL_VERSION};
use crate::parsing::noodles_helper::get_vcf_sample_name;

#[derive(Args, Clone, Default, Serialize)]
#[clap(author, about, 
    after_help = &**AFTER_HELP
)]
pub struct CompareSettings {
    #[clap(default_value = "")]
    #[clap(hide = true)]
    aardvark_version: String,

    /// Reference FASTA file
    #[clap(required = true)]
    #[clap(short = 'r')]
    #[clap(long = "reference")]
    #[clap(value_name = "FASTA")]
    #[clap(help_heading = Some("Input/Output"))]
    pub reference_fn: PathBuf,

    /// Truth variant call file (VCF)
    #[clap(required = true)]
    #[clap(short = 't')]
    #[clap(long = "truth-vcf")]
    #[clap(value_name = "VCF")]
    #[clap(help_heading = Some("Input/Output"))]
    pub truth_vcf_filename: PathBuf,

    /// Query variant call file (VCF)
    #[clap(required = true)]
    #[clap(short = 'q')]
    #[clap(long = "query-vcf")]
    #[clap(value_name = "VCF")]
    #[clap(help_heading = Some("Input/Output"))]
    pub query_vcf_filename: PathBuf,

    /// Confidence regions (BED)
    #[clap(short = 'b')]
    #[clap(long = "regions")]
    #[clap(value_name = "BED")]
    #[clap(help_heading = Some("Input/Output"))]
    pub regions: Option<PathBuf>,

    /// Stratifications, specifically the root file-of-filenames TSV
    #[clap(short = 's')]
    #[clap(long = "stratification")]
    #[clap(value_name = "TSV")]
    #[clap(help_heading = Some("Input/Output"))]
    pub stratifications: Option<PathBuf>,

    /// Output directory containing summary and VCFs
    #[clap(required = true)]
    #[clap(short = 'o')]
    #[clap(long = "output-dir")]
    #[clap(value_name = "DIR")]
    #[clap(help_heading = Some("Input/Output"))]
    pub output_folder: PathBuf,

    /// Optional output debug folder
    #[clap(long = "output-debug")]
    #[clap(value_name = "DIR")]
    #[clap(help_heading = Some("Input/Output"))]
    pub debug_folder: Option<PathBuf>,

    /// Optional comparison label for the summary output
    #[clap(long = "compare-label")]
    #[clap(value_name = "LABEL")]
    #[clap(help_heading = Some("Input/Output"))]
    #[clap(default_value = "compare")]
    pub compare_label: String,

    /// The sample name to use in the truth VCF [default: first sample]
    #[clap(long = "truth-sample")]
    #[clap(value_name = "SAMPLE")]
    #[clap(help_heading = Some("Input/Output"))]
    #[clap(default_value = "", hide_default_value = true)]
    pub truth_sample: String,

    /// The sample name to use in the query VCF [default: first sample]
    #[clap(long = "query-sample")]
    #[clap(value_name = "SAMPLE")]
    #[clap(help_heading = Some("Input/Output"))]
    #[clap(default_value = "", hide_default_value = true)]
    pub query_sample: String,

    /// The minimum gap (bp) between variants to split into separate sub-regions
    #[clap(long = "min-variant-gap")]
    #[clap(value_name = "BP")]
    #[clap(help_heading = Some("Region generation"))]
    #[clap(default_value = "50")]
    pub min_variant_gap: usize,

    /// Disables variant trimming, which may have a negative impact on accuracy
    #[clap(long = "disable-variant-trimming")]
    #[clap(help_heading = Some("Region generation"))]
    pub disable_variant_trimming: bool,

    /// Maximum edit distance in the WFA comparison before quitting
    #[clap(long = "max-edit-distance")]
    #[clap(value_name = "INT")]
    #[clap(help_heading = Some("Compare parameters"))]
    #[clap(default_value = "5000")]
    #[clap(hide = true)] // if you remove this, make sure you re-enable the CLI outputs
    pub max_edit_distance: usize,

    /// Maximum branch factor in the query optimizer; limits work on dense variant regions
    #[clap(long = "max-branch-factor")]
    #[clap(value_name = "INT")]
    #[clap(help_heading = Some("Compare parameters"))]
    #[clap(default_value = "50")]
    pub max_branch_factor: usize,

    /// Enables an exact-match compute shortcut at the cost of variant-level assessment accuracy
    #[clap(long = "enable-exact-shortcut")]
    #[clap(help_heading = Some("Compare parameters"))]
    #[clap(hide = true)] // if you remove this, make sure you re-enable the CLI outputs
    pub enable_exact_shortcut: bool,

    /// Enables the haplotype scoring metrics
    #[clap(long = "enable-haplotype-metrics")]
    #[clap(help_heading = Some("Optional metrics"))]
    #[clap(default_value = "false")]
    pub enable_haplotype_scoring: bool,

    /// Enables the weighted haplotype scoring metrics
    #[clap(long = "enable-weighted-haplotype-metrics")]
    #[clap(help_heading = Some("Optional metrics"))]
    #[clap(default_value = "false")]
    pub enable_weighted_haplotype_scoring: bool,

    /// Enables the record-basepair scoring metrics
    #[clap(long = "enable-record-basepair-metrics")]
    #[clap(help_heading = Some("Optional metrics"))]
    #[clap(default_value = "false")]
    pub enable_record_basepair_scoring: bool,

    /// Number of threads to use in the benchmarking step
    #[clap(long = "threads")]
    #[clap(value_name = "THREADS")]
    #[clap(default_value = "1")]
    pub threads: usize,

    /// Enable verbose output.
    #[clap(short = 'v')]
    #[clap(long = "verbose")]
    #[clap(action = clap::ArgAction::Count)]
    pub verbosity: u8,

    // Debug options that are generally hidden and just for quick testing
    /// Skips a number of sub-problems (debug only); non-0 values may create partial values in output
    #[clap(hide = true)]
    #[clap(long = "skip")]
    #[clap(default_value = "0")]
    pub skip_blocks: usize,

    /// Takes a number of sub-problems (debug only); non-0 values may create partial values in output
    #[clap(hide = true)]
    #[clap(long = "take")]
    #[clap(default_value = "0")]
    pub take_blocks: usize,
}

pub fn check_compare_settings(mut settings: CompareSettings) -> anyhow::Result<CompareSettings> {
    // hard code the version in
    settings.aardvark_version = FULL_VERSION.clone();
    info!("Aardvark version: {:?}", &settings.aardvark_version);
    info!("Sub-command: compare");
    info!("Inputs:");

    // check for all the required input files
    check_required_filename(&settings.reference_fn, "Reference FASTA")?;
    check_required_filename(&settings.truth_vcf_filename, "Truth VCF")?;
    check_required_filename(&settings.query_vcf_filename, "Query VCF")?;
    check_optional_filename(settings.regions.as_deref(), "Regions")?;
    check_optional_filename(settings.stratifications.as_deref(), "Stratifications")?;
    
    // dump stuff to the logger
    info!("\tReference: {:?}", &settings.reference_fn);
    info!("\tTruth VCF: {:?}", &settings.truth_vcf_filename);
    if settings.truth_sample.is_empty() {
        settings.truth_sample = get_vcf_sample_name(&settings.truth_vcf_filename, 0)?;
    }
    info!("\tTruth sample: {:?}", &settings.truth_sample);
    info!("\tQuery VCF: {:?}", &settings.query_vcf_filename);
    if settings.query_sample.is_empty() {
        settings.query_sample = get_vcf_sample_name(&settings.query_vcf_filename, 0)?;
    }
    info!("\tQuery sample: {:?}", &settings.query_sample);
    if let Some(hcr_fn) = settings.regions.as_deref() {
        info!("\tRegions: {hcr_fn:?}");
    } else {
        info!("\tRegions: None");
    }
    if let Some(filename) = settings.stratifications.as_deref() {
        info!("\tStratifications: {filename:?}");
    } else {
        info!("\tStratifications: None");
    }

    // check for the haplotype scoring metrics
    if settings.enable_haplotype_scoring || settings.enable_weighted_haplotype_scoring || settings.enable_record_basepair_scoring {
        info!("Optional metrics:");
        info!("\tHaplotype scoring: {}", if settings.enable_haplotype_scoring { "ENABLED" } else { "DISABLED" });
        info!("\tWeighted haplotype scoring: {}", if settings.enable_weighted_haplotype_scoring { "ENABLED" } else { "DISABLED" });
        info!("\tRecord-basepair scoring: {}", if settings.enable_record_basepair_scoring { "ENABLED" } else { "DISABLED" });
    }

    // 0 is just a sentinel for everything
    if settings.take_blocks == 0 {
        settings.take_blocks = usize::MAX;
    }
    
    // outputs
    info!("Outputs:");
    info!("\tCompare label: {:?}", &settings.compare_label);
    info!("\tOutput folder: {:?}", &settings.output_folder);
    if let Some(debug_folder) = settings.debug_folder.as_ref() {
        info!("\tDebug folder: {debug_folder:?}");
    }

    // other misc parameters
    info!("Region generation parameters:");
    ensure!(settings.min_variant_gap > 0, "--min-variant-gap must be >0");
    info!("\tMinimum variant gap: {}", settings.min_variant_gap);
    info!("\tVariant trimming: {}", if settings.disable_variant_trimming { "DISABLED "} else { "ENABLED" });

    info!("Compare parameters:");
    ensure!(settings.max_branch_factor > 0, "--max-branch-factor must be >0");
    info!("\tMax branch factor: {}", settings.max_branch_factor);
    if settings.enable_exact_shortcut {
        info!("\tExact match shortcut: {}", if settings.enable_exact_shortcut { "ENABLED" } else { "DISABLED" });
    }

    if settings.threads == 0 {
        settings.threads = 1;
    }
    info!("Processing threads: {}", settings.threads);

    Ok(settings)
}