v0.1.0

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
 

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "bronko"
-version = "0.0.3"
+version = "0.1.0"
 edition = "2021"
 authors = [
     "Ryan Doughty <rdd4@rice.edu>",

src/call.rs

@@ -100,6 +100,34 @@ 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 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 {
+        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,39 +1055,56 @@ 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(); 
 
-                    // 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;
+                        // filter out if greater than strand odds ratio (default 8)
+                        if sor > args.strand_odds_max {
+                            continue;
+                        }
 
-                    if c_k < n_kmer_per_strand && d_k < n_kmer_per_strand {
-                        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;
+
+                        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
                     }
                 }
 
                 //Get minor af (filter out if below reporting threshold)
                 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 +1115,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;
                 }
 

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")]

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";