From e36a35d4893b05e33268532f272b12a59f719bd5 Mon Sep 17 00:00:00 2001 From: Ryan Doughty Date: Tue, 11 Nov 2025 11:12:54 -0600 Subject: [PATCH 1/3] parameter updates and new variant calling arguments --- Cargo.toml | 2 +- src/call.rs | 80 +++++++++++++++++++++++++++++++++++++++------------ src/cli.rs | 22 ++++++++++++-- src/consts.rs | 11 ++++--- 4 files changed, 89 insertions(+), 26 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 4ac55b4..5bb2b2e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "bronko" -version = "0.0.3" +version = "0.1.0" edition = "2021" authors = [ "Ryan Doughty ", diff --git a/src/call.rs b/src/call.rs index c143e44..9e986f4 100644 --- a/src/call.rs +++ b/src/call.rs @@ -100,6 +100,30 @@ fn check_args(args: &CallArgs) { warn!("Number of kmers per strand set very high, only strongly supported variants will be returned") } + + if args.strand_balance_ratio < 0.0 { + error!("Strand balance ratio is set to below 0, must be between 0.0 and 1.0"); + std::process::exit(1); + } else if args.strand_balance_ratio > 1.0 { + error!("Strand balance ratio is set above 1, must be between 0.0 and 1.0"); + std::process::exit(1); + } else if args.strand_balance_ratio == 1.0 { + warn!("Strand balance ratio is set to 1, all variants will pass this filter"); + } + + if args.min_variant_depth < 0 { + warn!("Minimum variant depth set below 0, all variants will be returned"); + } + + if args.variant_multiplier < 1.0 { + error!("Noise multiplier for variant detection is set to below 1.0, must be greater than 1.0 (recommended between 1.3-2.0)"); + std::process::exit(1); + } else if args.variant_multiplier > 2.0 { + warn!("Strand balance ratio is set above 2, may experience a drop in recall (we recommend ~1.5)"); + } else if args.variant_multiplier == 1.0 { + warn!("Noise multiplier for variant detection set to 1.0, all variants will pass this filter"); + } + if args.first_pairs.len() != args.second_pairs.len() { error!("Number of paired end sequences do not match, exiting."); std::process::exit(1); @@ -1027,31 +1051,43 @@ pub fn call_variants( } // NEW Strand filter logic - let mut sor = args.strand_odds_max + 1.0; + let mut sor = args.strand_odds_max + 1.0; //default is above the given so filtered if not able to calculate if strand_filter { let a = row[ref_base as usize] as f64 + 1.0; //ref fwd let b = row_rev[ref_base as usize] as f64 + 1.0; //ref rev let c = row[alt_base as usize] as f64 + 1.0; //alt fwd let d = row_rev[alt_base as usize] as f64 + 1.0; //alt rev - //Using GATK strand odds ratio - let r = (a*d)/(b*c); - let ref_ratio = (a.min(b)) / (a.max(b)); - let alt_ratio = (c.min(d)) / (c.max(d)); - - sor = (r + (1.0 / r)).ln() + ref_ratio.ln() - alt_ratio.ln(); + //calculate the difference between the strands + let ref_total = a + b + c + d; + let min_strand_depth = (a+c).min(b+d); + let min_strand_percent = min_strand_depth / ref_total; - // filter out if greater than strand odds ratio (default 2) - if sor > args.strand_odds_max { - continue; - } + //if the strand balance filter is on (default), then always do SOR filter + //if the strand balance filter is off, then only do SOR filter when 1 strand is > strand_balance_ratio of total depth (by default 10% of total depth, aka all cases where 1 strand is less than 10% of the total depth are ignored) + if (!args.no_strand_balance_filter) | ((args.no_strand_balance_filter) & (min_strand_percent >= args.strand_balance_ratio)) { + + //Using GATK strand odds ratio + let r = (a*d)/(b*c); + let ref_ratio = (a.min(b)) / (a.max(b)); + let alt_ratio = (c.min(d)) / (c.max(d)); + + sor = (r + (1.0 / r)).ln() + ref_ratio.ln() - alt_ratio.ln(); + + // filter out if greater than strand odds ratio (default 8) + if sor > args.strand_odds_max { + continue; + } - // additional filtering for low kmer support across both strands - let c_k = count[alt_base as usize] as usize; - let d_k = count_rev[alt_base as usize] as usize; + // additional filtering for low kmer support across both strands + let c_k = count[alt_base as usize] as usize; + let d_k = count_rev[alt_base as usize] as usize; - if c_k < n_kmer_per_strand && d_k < n_kmer_per_strand { - continue; + if c_k < n_kmer_per_strand && d_k < n_kmer_per_strand { + continue; + } + } else { + sor = -1.0; //output SOR == -1.0 if the strand is unbalanced so not tested } } @@ -1059,7 +1095,12 @@ pub fn call_variants( let alt_count = row_total[alt_base as usize]; let af = alt_count as f64 / total_depth as f64; - if af < min_af || af < (2.0*baseline_noise[i].max) { + let y0: f64 = args.variant_multiplier; + let p0: f64 = 0.5; + let a: f64 = 0.03; + let factor = y0 + p0 * a.powf(100.0 * af); + + if af < min_af || af < (factor.max(y0) * baseline_noise[i].max) { continue; } @@ -1070,7 +1111,10 @@ pub fn call_variants( } else { if total_depth < args.min_depth as u64 { // filter out minor variants when the total depth is too low continue; - } + } + if alt_count < args.min_variant_depth as u64 { //filter out minor variants when variant depth is below threshold + continue; + } num_minor_variants += 1; } diff --git a/src/cli.rs b/src/cli.rs index 93626b8..9120de1 100644 --- a/src/cli.rs +++ b/src/cli.rs @@ -93,7 +93,7 @@ pub struct CallArgs { pub use_full_kmer: bool, //number of buckets to ignore on ends of kmers - #[clap(long="n-fixed", default_value_t = DEFAULT_N_FIXED, help_heading="ALGORITHM", help="Number of fixed positions at the end of each kmer")] + #[clap(long="n-fixed", default_value_t = DEFAULT_N_FIXED, help_heading="ALGORITHM", help="Number of fixed positions at the end of each kmer that cannot contribute to pileup")] pub n_fixed: usize, //VARIANT CALLING PARAMETERS @@ -106,9 +106,17 @@ pub struct CallArgs { pub no_end_filter: bool, //do not do strand filtering - #[clap(long="no-strand-filter", default_value_t = DEFAULT_NO_STRAND_FILTER, help_heading="VARIANT CALLING PARAMETERS", help="Do not filter variants that are present on one strand but not the other")] + #[clap(long="no-strand-filter", default_value_t = DEFAULT_NO_STRAND_FILTER, help_heading="VARIANT CALLING PARAMETERS", help="Do not utilize SOR test to filter variants that are present on one strand but not the other")] pub no_strand_filter: bool, + //do not filter variants with extreme strand balance issues + #[clap(long="no-strand-balance-filter", default_value_t = DEFAULT_NO_STRAND_BALANCE_FILTER, help_heading="VARIANT CALLING PARAMETERS", help="Allow variants with extreme strand disbalance pass without SOR check (will not matter if --no-strand-filter present)")] + pub no_strand_balance_filter: bool, + + //strand balance ratio for strand to be ignored + #[clap(long="balance-ratio", default_value_t = DEFAULT_STRAND_BALANCE_RATIO, help_heading="VARIANT CALLING PARAMETERS", help="Percent of total depth that one strand must be under to be considered unbalanced (must be [0.0-1.0])")] + pub strand_balance_ratio: f64, + //the number of kmers per strand that are required to call a variant #[clap(long="n-per-strand", default_value_t = DEFAULT_N_KMERS_PER_STRAND, help_heading="VARIANT CALLING PARAMETERS", help="Min number of unique kmers to observe to call a variant at any site (needed on both strands if strand filter active)")] pub n_per_strand: usize, @@ -117,9 +125,17 @@ pub struct CallArgs { pub strand_odds_max: f64, //min depth to call a variant - #[clap(long="min-depth", default_value_t = DEFAULT_MIN_DEPTH, help_heading="VARIANT CALLING PARAMETERS", help="Minimum total depth at an allele to call a minor variant (default=50*min_kmers)")] + #[clap(long="min-depth", default_value_t = DEFAULT_MIN_DEPTH, help_heading="VARIANT CALLING PARAMETERS", help="Minimum total depth at an allele to call a minor variant (default=100*min_kmers)")] pub min_depth: usize, + //minimum reads/kmers supporting a minor variant + #[clap(long="min-variant-depth", default_value_t = MIN_KMER_COUNT, help_heading = "VARIANT CALLING PARAMETERS", help="Minimum depth of a minor variant to be called present (default=min_kmers)")] + pub min_variant_depth: usize, + + //noise multiplier + #[clap(long="noise-multiplier", default_value_t = DEFAULT_NOISE_MULTIPLIER, help_heading = "VARIANT CALLING PARAMETERS", help="How much greater (1x, 1.5x, etc) the minor allele frequency of a variant must be above estimated baseline noise in that region (must be > 1.0x). Note that for variants under 1%, multiplier will be increased exponentially up to +0.5 more")] + pub variant_multiplier: f64, + //OUTPUT PARAMETERS //todo add output locations, output formats #[clap(short, long="output", help_heading="OUTPUT", default_value = DEFAULT_OUT_FOLDER, help="Folder to output all resulting files")] diff --git a/src/consts.rs b/src/consts.rs index fb3c46f..6ca8ed6 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -2,17 +2,20 @@ pub const BRONKO_VERSION: &str = env!("CARGO_PKG_VERSION"); pub const MIN_KMER_SIZE: usize = 15; pub const DEFAULT_KMER_SIZE: usize = 21; //kmer size pub const MAX_KMER_SIZE: usize = 31; -pub const MIN_KMER_COUNT: usize = 10; // minimum number of times a kmer must occur to be run through algorithm +pub const MIN_KMER_COUNT: usize = 3; // minimum number of times a kmer must occur to be run through algorithm pub const DEFAULT_MIN_AF: f64 = 0.03; // minimum allele frequency to be reported pub const DEFAULT_NO_FILTER_ENDS: bool = false; //flag to not filter variants on the ends of sequences pub const DEFAULT_NO_STRAND_FILTER: bool = false; // do not use a strand filter +pub const DEFAULT_NO_STRAND_BALANCE_FILTER: bool = false; //do not filter variants with extremely unbalanced strands +pub const DEFAULT_STRAND_BALANCE_RATIO: f64 = 0.1; //percent of total depth that one strand must be under to be called extremely unbalanced pub const DEFAULT_N_KMERS_PER_STRAND: usize = 2; //number of kmers that are needed to support each strand to call a variant pub const DEFAULT_MAX_STRAND_ODDS: f64 = 6.0; // Strand odds ratio that must be under +pub const DEFAULT_NOISE_MULTIPLIER: f64 = 1.5; //How much variant must be above baseline noise (with exponential increase below 1% (going up to + 0.3 higher)) pub const DEFAULT_TSV_PILEUP: bool = false; //print a pileup as well -pub const DEFAULT_ALIGNMENT: bool = false; -pub const DEFAULT_KEEP_KMER_INFO: bool = false; +pub const DEFAULT_ALIGNMENT: bool = false; // return an alignment +pub const DEFAULT_KEEP_KMER_INFO: bool = false; //also print out kmer counts (deleted by default) pub const DEFAULT_N_FIXED: usize = 2; //number of fixed bases in the kmer pub const DEFAULT_USE_FULL_KMER: bool = false; //use the full length of the kmer -pub const DEFAULT_MIN_DEPTH: usize = MIN_KMER_COUNT*50; +pub const DEFAULT_MIN_DEPTH: usize = MIN_KMER_COUNT*100; pub const DEFAULT_INDEX_OUTPUT: &str = "bronko"; pub const DEFAULT_OUT_FOLDER: &str = "bronko_output"; \ No newline at end of file From 6bcbd4135f23d8ef0bb04afa5c541674cd16587b Mon Sep 17 00:00:00 2001 From: Ryan Doughty Date: Tue, 11 Nov 2025 11:16:07 -0600 Subject: [PATCH 2/3] small warning for depth argument --- src/call.rs | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/call.rs b/src/call.rs index 9e986f4..2730ac9 100644 --- a/src/call.rs +++ b/src/call.rs @@ -112,7 +112,11 @@ fn check_args(args: &CallArgs) { } if args.min_variant_depth < 0 { - warn!("Minimum variant depth set below 0, all variants will be returned"); + warn!("Minimum variant depth set below 0, all variants will be returned if passing other thresholds"); + } + + if args.min_depth < 0 { + warn!("Minimum total depth for minor variant calling set below 0, all variants will be returned if passing other thresholds"); } if args.variant_multiplier < 1.0 { From 5b41fb906cbe3acf6fd06ad2b8dc3fc09d7253a0 Mon Sep 17 00:00:00 2001 From: Ryan Doughty Date: Tue, 11 Nov 2025 11:19:12 -0600 Subject: [PATCH 3/3] Changelog for v0.1.0 --- CHANGELOG | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index d2573b4..cea290f 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -2,6 +2,11 @@ All major changes to bronko will be documented in this file +## [v0.1.0] - 11/11/2025 +- First major release of bronko, corresponding to first manuscript +- Updated defaults for multiple variant calling parameters, introduced multiple others to allow more control over variant calls +- By default now significant strand balance issues will be filtered + ## [v0.0.3] - 10/28/2025 - Introduced baseline noise filter using streaming version of modified thompson tau filter