|
|
@@ -22,12 +22,25 @@
|
|
|
//! - 12 columns (extended, optional):
|
|
|
//! 12. fb_invs (comma-separated per-base counts)
|
|
|
//!
|
|
|
-//! Parsers accept either 11 or 12 columns; if the 12th column is absent, `fb_invs` is empty.
|
|
|
+//! - 13 columns (extended, optional):
|
|
|
+//! 13. depths_sens_ratio (comma-separated per-base fwd/(fwd+rev) ratios, `NaN` where total depth is 0)
|
|
|
+//!
|
|
|
+//! - 16 columns (extended, current):
|
|
|
+//! 14. ratio_lowq (sum(low_qualities) / sum(depths) for the bin)
|
|
|
+//! 15. ratio_fb (max per-base fb_invs count / depth at that same position)
|
|
|
+//! 16. ratio_sens (pooled sum(fwd) / sum(fwd+rev) for the bin)
|
|
|
+//!
|
|
|
+//! Parsers accept 11, 12, 13, or 16 columns. Column counts of 14 or 15 are rejected as invalid —
|
|
|
+//! fields 13-16 are only ever written together, so a partial set indicates a corrupt or
|
|
|
+//! hand-edited file rather than an older schema version.
|
|
|
//!
|
|
|
//! ## Semantics (bitwise-identical mode)
|
|
|
//!
|
|
|
//! - **Per-base vectors** (`depths`, `low_qualities`, `fb_invs`) are accumulated **per alignment**
|
|
|
-//! for each bin it overlaps.
|
|
|
+//! for each bin it overlaps. `depths` is derived from separately tracked forward/reverse
|
|
|
+//! per-base counts (see `depths_sens_ratio` below); it is not accumulated directly.
|
|
|
+//! - **`depths_sens_ratio`** is a per-base vector: for each position, `fwd / (fwd + rev)` where
|
|
|
+//! `fwd`/`rev` are the forward/reverse-strand depth at that position, or `NaN` if depth is 0.
|
|
|
//! - **`n_reads` and ratios** are computed using a per-bin **qname-deduplicated representative record**:
|
|
|
//! the representative is the **last** alignment encountered for a given qname within that bin
|
|
|
//! (“last wins”), matching the previous `HashMap` overwrite behavior.
|
|
|
@@ -35,6 +48,10 @@
|
|
|
//! - `ratio_sa` = (number of representative records with an SA tag) / `n_reads`
|
|
|
//! - `ratio_start` = (maximum representative-start pileup at a single position) / `n_reads`
|
|
|
//! - `ratio_end` = (maximum representative-end pileup at a single position) / `n_reads`
|
|
|
+//! - `ratio_lowq` = sum(`low_qualities`) / sum(`depths`) over the bin, pooled (not per-position mean)
|
|
|
+//! - `ratio_fb` = max per-base `fb_invs` count / `depths` at that same position
|
|
|
+//! - `ratio_sens` = pooled sum(fwd) / sum(fwd + rev) over the bin (not a mean of
|
|
|
+//! `depths_sens_ratio`, to avoid over-weighting low-depth positions)
|
|
|
//!
|
|
|
//! ## Performance notes
|
|
|
//! - The largest win comes from avoiding thousands of small random indexed reads.
|
|
|
@@ -69,7 +86,6 @@ use crate::config::Config;
|
|
|
use crate::pipes::{Initialize, ShouldRun};
|
|
|
use crate::positions::GenomeRange;
|
|
|
use crate::runners::Run;
|
|
|
-use crate::scan::bin::{Bin, BinStats};
|
|
|
use crate::variant::vcf_variant::Label;
|
|
|
|
|
|
/// Represents a count of reads in a genomic bin, including various metrics and outlier information.
|
|
|
@@ -91,6 +107,12 @@ pub struct BinCount {
|
|
|
pub ratio_start: f64,
|
|
|
/// Ratio of max end pileup (representative reads) to total representative reads.
|
|
|
pub ratio_end: f64,
|
|
|
+ /// sum(low_qualities) / sum(depths)
|
|
|
+ pub ratio_lowq: f64,
|
|
|
+ /// max(fb_invs) / depths at that same position
|
|
|
+ pub ratio_fb: f64,
|
|
|
+ /// pooled sum(fwd) / sum(fwd+rev), bin-level
|
|
|
+ pub ratio_sens: f64,
|
|
|
/// Optional vector of outlier types for this bin.
|
|
|
pub outlier: Option<Vec<BinOutlier>>,
|
|
|
/// Per-base depth over the bin (comma-separated in TSV).
|
|
|
@@ -99,87 +121,91 @@ pub struct BinCount {
|
|
|
pub low_qualities: Vec<u32>,
|
|
|
/// Per-base fb_inv counts over the bin (comma-separated in TSV; optional column).
|
|
|
pub fb_invs: Vec<u32>,
|
|
|
+ /// fwd / (fwd + rev) per position, NaN if 0 total
|
|
|
+ pub depths_sens_ratio: Vec<f64>,
|
|
|
}
|
|
|
|
|
|
-impl From<&Bin> for BinCount {
|
|
|
- /// Converts a `Bin` reference to a `BinCount`.
|
|
|
- ///
|
|
|
- /// # Parameters
|
|
|
- /// - `bin: &Bin`: A reference to the `Bin` object to convert.
|
|
|
- ///
|
|
|
- /// # Returns
|
|
|
- /// A new `BinCount` instance populated with data from the `Bin`.
|
|
|
- fn from(bin: &Bin) -> Self {
|
|
|
- let n_reads = bin.n_reads();
|
|
|
- let (n_sa, n_start, n_end) = bin.count_reads_sa_start_end();
|
|
|
- let (ratio_sa, ratio_start, ratio_end) = if n_reads == 0 {
|
|
|
- (f64::NAN, f64::NAN, f64::NAN)
|
|
|
- } else {
|
|
|
- let n_reads_float = n_reads as f64;
|
|
|
- (
|
|
|
- n_sa as f64 / n_reads_float,
|
|
|
- n_start as f64 / n_reads_float,
|
|
|
- n_end as f64 / n_reads_float,
|
|
|
- )
|
|
|
- };
|
|
|
-
|
|
|
- Self {
|
|
|
- contig: bin.contig.clone(),
|
|
|
- start: bin.start,
|
|
|
- end: bin.end,
|
|
|
- n_reads: n_reads as u32,
|
|
|
- coverage: bin.mean_coverage_from_depths(),
|
|
|
- ratio_sa,
|
|
|
- ratio_start,
|
|
|
- ratio_end,
|
|
|
- outlier: None,
|
|
|
- depths: bin.depths.clone(),
|
|
|
- low_qualities: bin.low_qualities.clone(),
|
|
|
- fb_invs: bin.fb_invs.clone(),
|
|
|
- }
|
|
|
- }
|
|
|
-}
|
|
|
-
|
|
|
-impl From<&BinStats> for BinCount {
|
|
|
- fn from(bin: &BinStats) -> Self {
|
|
|
- let n_reads_f = bin.n_reads as f64;
|
|
|
-
|
|
|
- let (ratio_sa, ratio_start, ratio_end) = if bin.n_reads == 0 {
|
|
|
- (f64::NAN, f64::NAN, f64::NAN)
|
|
|
- } else {
|
|
|
- (
|
|
|
- bin.n_sa as f64 / n_reads_f,
|
|
|
- bin.max_start_count as f64 / n_reads_f,
|
|
|
- bin.max_end_count as f64 / n_reads_f,
|
|
|
- )
|
|
|
- };
|
|
|
-
|
|
|
- Self {
|
|
|
- contig: bin.contig.clone(),
|
|
|
- start: bin.start,
|
|
|
- end: bin.end,
|
|
|
- n_reads: bin.n_reads,
|
|
|
- coverage: bin.mean_coverage_from_depths(),
|
|
|
- ratio_sa,
|
|
|
- ratio_start,
|
|
|
- ratio_end,
|
|
|
- outlier: None,
|
|
|
- depths: bin.depths.clone(),
|
|
|
- low_qualities: bin.low_qualities.clone(),
|
|
|
- fb_invs: bin.fb_invs.clone(),
|
|
|
- }
|
|
|
- }
|
|
|
-}
|
|
|
+// impl From<&Bin> for BinCount {
|
|
|
+// /// Converts a `Bin` reference to a `BinCount`.
|
|
|
+// ///
|
|
|
+// /// # Parameters
|
|
|
+// /// - `bin: &Bin`: A reference to the `Bin` object to convert.
|
|
|
+// ///
|
|
|
+// /// # Returns
|
|
|
+// /// A new `BinCount` instance populated with data from the `Bin`.
|
|
|
+// fn from(bin: &Bin) -> Self {
|
|
|
+// let n_reads = bin.n_reads();
|
|
|
+// let (n_sa, n_start, n_end) = bin.count_reads_sa_start_end();
|
|
|
+// let (ratio_sa, ratio_start, ratio_end) = if n_reads == 0 {
|
|
|
+// (f64::NAN, f64::NAN, f64::NAN)
|
|
|
+// } else {
|
|
|
+// let n_reads_float = n_reads as f64;
|
|
|
+// (
|
|
|
+// n_sa as f64 / n_reads_float,
|
|
|
+// n_start as f64 / n_reads_float,
|
|
|
+// n_end as f64 / n_reads_float,
|
|
|
+// )
|
|
|
+// };
|
|
|
+//
|
|
|
+// Self {
|
|
|
+// contig: bin.contig.clone(),
|
|
|
+// start: bin.start,
|
|
|
+// end: bin.end,
|
|
|
+// n_reads: n_reads as u32,
|
|
|
+// coverage: bin.mean_coverage_from_depths(),
|
|
|
+// ratio_sa,
|
|
|
+// ratio_start,
|
|
|
+// ratio_end,
|
|
|
+// outlier: None,
|
|
|
+// depths: bin.depths.clone(),
|
|
|
+// low_qualities: bin.low_qualities.clone(),
|
|
|
+// fb_invs: bin.fb_invs.clone(),
|
|
|
+// }
|
|
|
+// }
|
|
|
+// }
|
|
|
+
|
|
|
+// impl From<&BinStats> for BinCount {
|
|
|
+// fn from(bin: &BinStats) -> Self {
|
|
|
+// let n_reads_f = bin.n_reads as f64;
|
|
|
+//
|
|
|
+// let (ratio_sa, ratio_start, ratio_end) = if bin.n_reads == 0 {
|
|
|
+// (f64::NAN, f64::NAN, f64::NAN)
|
|
|
+// } else {
|
|
|
+// (
|
|
|
+// bin.n_sa as f64 / n_reads_f,
|
|
|
+// bin.max_start_count as f64 / n_reads_f,
|
|
|
+// bin.max_end_count as f64 / n_reads_f,
|
|
|
+// )
|
|
|
+// };
|
|
|
+//
|
|
|
+// Self {
|
|
|
+// contig: bin.contig.clone(),
|
|
|
+// start: bin.start,
|
|
|
+// end: bin.end,
|
|
|
+// n_reads: bin.n_reads,
|
|
|
+// coverage: bin.mean_coverage_from_depths(),
|
|
|
+// ratio_sa,
|
|
|
+// ratio_start,
|
|
|
+// ratio_end,
|
|
|
+// outlier: None,
|
|
|
+// depths: bin.depths.clone(),
|
|
|
+// low_qualities: bin.low_qualities.clone(),
|
|
|
+// fb_invs: bin.fb_invs.clone(),
|
|
|
+// }
|
|
|
+// }
|
|
|
+// }
|
|
|
|
|
|
impl BinCount {
|
|
|
/// Converts the `BinCount` instance to a TSV (tab-separated) row.
|
|
|
///
|
|
|
/// ## Columns
|
|
|
- /// Always writes **12 columns** (including `fb_invs` as the last column).
|
|
|
- /// Downstream parsers may accept 11 columns for legacy files without `fb_invs`.
|
|
|
+ /// Always writes **16 columns**, including `fb_invs` (12), `depths_sens_ratio` (13), and the
|
|
|
+ /// pooled `ratio_lowq`/`ratio_fb`/`ratio_sens` scalars (14–16). Older parsers reading these
|
|
|
+ /// files should tolerate 11 or 12 columns for legacy data written before those fields existed;
|
|
|
+ /// this method itself never omits them.
|
|
|
pub fn to_tsv_row(&self) -> String {
|
|
|
format!(
|
|
|
- "{}\t{}\t{}\t{}\t{:.6}\t{:.6}\t{:.6}\t{:.6}\t{}\t{}\t{}\t{}",
|
|
|
+ "{}\t{}\t{}\t{}\t{:.6}\t{:.6}\t{:.6}\t{:.6}\t{}\t{}\t{}\t{}\t{}\t{:.6}\t{:.6}\t{:.6}",
|
|
|
self.contig,
|
|
|
self.start,
|
|
|
self.end,
|
|
|
@@ -211,27 +237,39 @@ impl BinCount {
|
|
|
.map(|e| e.to_string())
|
|
|
.collect::<Vec<_>>()
|
|
|
.join(","),
|
|
|
+ self.depths_sens_ratio
|
|
|
+ .iter()
|
|
|
+ .map(|v| format!("{v:.6}"))
|
|
|
+ .collect::<Vec<_>>()
|
|
|
+ .join(","),
|
|
|
+ self.ratio_lowq,
|
|
|
+ self.ratio_fb,
|
|
|
+ self.ratio_sens,
|
|
|
)
|
|
|
}
|
|
|
|
|
|
/// Parses a TSV row and creates a [`BinCount`].
|
|
|
///
|
|
|
- /// Accepts either:
|
|
|
- /// - **11 columns** (legacy; without `fb_invs`)
|
|
|
- /// - **12 columns** (extended; with `fb_invs`)
|
|
|
+ /// Accepts:
|
|
|
+ /// - **11 columns** (legacy; without `fb_invs`, `depths_sens_ratio`, or the pooled ratios)
|
|
|
+ /// - **12 columns** (adds `fb_invs`)
|
|
|
+ /// - **13 columns** (adds `depths_sens_ratio`)
|
|
|
+ /// - **16 columns** (adds `ratio_lowq`, `ratio_fb`, `ratio_sens`; current format)
|
|
|
+ ///
|
|
|
+ /// Column counts of 14 or 15 are rejected — fields 14-16 are only ever written together with
|
|
|
+ /// field 13, so a partial set means the row is corrupt rather than an older schema version.
|
|
|
///
|
|
|
/// ## Errors
|
|
|
/// Returns an error if:
|
|
|
- /// - Column count is not 11 or 12
|
|
|
+ /// - Column count is not one of 11, 12, 13, or 16
|
|
|
/// - Any numeric field fails to parse
|
|
|
/// - Outlier labels are not recognized
|
|
|
pub fn from_tsv_row(row: &str) -> anyhow::Result<Self> {
|
|
|
let fields: Vec<&str> = row.split('\t').collect();
|
|
|
- // wont block if fb_inv has been added or not
|
|
|
- if fields.len() != 11 && fields.len() != 12 {
|
|
|
+ // won't block if fb_inv and/or depths_sens_ratio have been added or not
|
|
|
+ if !matches!(fields.len(), 11 | 12 | 13 | 16) {
|
|
|
anyhow::bail!("Invalid number of fields in TSV row (got {})", fields.len());
|
|
|
}
|
|
|
-
|
|
|
let outlier = if !fields[8].is_empty() {
|
|
|
Some(
|
|
|
fields[8]
|
|
|
@@ -240,6 +278,9 @@ impl BinCount {
|
|
|
"SA" => Ok(BinOutlier::SA),
|
|
|
"Start" => Ok(BinOutlier::Start),
|
|
|
"End" => Ok(BinOutlier::End),
|
|
|
+ "LowQ" => Ok(BinOutlier::LowQ),
|
|
|
+ "Fb" => Ok(BinOutlier::Fb),
|
|
|
+ "Sens" => Ok(BinOutlier::Sens),
|
|
|
_ => Err(anyhow::anyhow!("Invalid outlier type: {}", s)),
|
|
|
})
|
|
|
.collect::<Result<Vec<BinOutlier>, _>>()?,
|
|
|
@@ -248,6 +289,16 @@ impl BinCount {
|
|
|
None
|
|
|
};
|
|
|
|
|
|
+ let (ratio_lowq, ratio_fb, ratio_sens) = if fields.len() == 16 {
|
|
|
+ (
|
|
|
+ parse_f64_allow_nan(fields[13])?,
|
|
|
+ parse_f64_allow_nan(fields[14])?,
|
|
|
+ parse_f64_allow_nan(fields[15])?,
|
|
|
+ )
|
|
|
+ } else {
|
|
|
+ (f64::NAN, f64::NAN, f64::NAN)
|
|
|
+ };
|
|
|
+
|
|
|
let depths = if fields[9].trim().is_empty() {
|
|
|
Vec::new()
|
|
|
} else {
|
|
|
@@ -260,7 +311,6 @@ impl BinCount {
|
|
|
})
|
|
|
.collect::<anyhow::Result<Vec<u32>>>()?
|
|
|
};
|
|
|
-
|
|
|
let low_qualities = if fields[10].trim().is_empty() {
|
|
|
Vec::new()
|
|
|
} else {
|
|
|
@@ -273,8 +323,7 @@ impl BinCount {
|
|
|
})
|
|
|
.collect::<anyhow::Result<Vec<u32>>>()?
|
|
|
};
|
|
|
-
|
|
|
- let fb_invs = if fields.len() == 12 {
|
|
|
+ let fb_invs = if fields.len() >= 12 {
|
|
|
if fields[11].trim().is_empty() {
|
|
|
Vec::new()
|
|
|
} else {
|
|
|
@@ -290,7 +339,22 @@ impl BinCount {
|
|
|
} else {
|
|
|
Vec::new()
|
|
|
};
|
|
|
-
|
|
|
+ let depths_sens_ratio = if fields.len() == 13 {
|
|
|
+ if fields[12].trim().is_empty() {
|
|
|
+ Vec::new()
|
|
|
+ } else {
|
|
|
+ fields[12]
|
|
|
+ .split(',')
|
|
|
+ .map(|s| {
|
|
|
+ s.trim()
|
|
|
+ .parse::<f64>() // NaN-safe: f64::from_str parses "NaN" natively
|
|
|
+ .with_context(|| format!("Failed to parse '{}' as f64", s))
|
|
|
+ })
|
|
|
+ .collect::<anyhow::Result<Vec<f64>>>()?
|
|
|
+ }
|
|
|
+ } else {
|
|
|
+ Vec::new()
|
|
|
+ };
|
|
|
Ok(BinCount {
|
|
|
contig: fields[0].to_string(),
|
|
|
start: fields[1].parse()?,
|
|
|
@@ -300,10 +364,14 @@ impl BinCount {
|
|
|
ratio_sa: parse_f64_allow_nan(fields[5])?,
|
|
|
ratio_start: parse_f64_allow_nan(fields[6])?,
|
|
|
ratio_end: parse_f64_allow_nan(fields[7])?,
|
|
|
+ ratio_lowq,
|
|
|
+ ratio_fb,
|
|
|
+ ratio_sens,
|
|
|
outlier,
|
|
|
depths,
|
|
|
low_qualities,
|
|
|
fb_invs,
|
|
|
+ depths_sens_ratio,
|
|
|
})
|
|
|
}
|
|
|
}
|
|
|
@@ -318,34 +386,6 @@ fn parse_f64_allow_nan(s: &str) -> anyhow::Result<f64> {
|
|
|
}
|
|
|
}
|
|
|
|
|
|
-/// Represents types of outliers that can be detected in a genomic bin.
|
|
|
-#[derive(Debug, Clone)]
|
|
|
-pub enum BinOutlier {
|
|
|
- /// Indicates an outlier in supplementary alignments.
|
|
|
- SA,
|
|
|
- /// Indicates an outlier in reads starting in this bin.
|
|
|
- Start,
|
|
|
- /// Indicates an outlier in reads ending in this bin.
|
|
|
- End,
|
|
|
-}
|
|
|
-
|
|
|
-impl fmt::Display for BinOutlier {
|
|
|
- /// Implements the `Display` trait for `BinOutlier`.
|
|
|
- ///
|
|
|
- /// # Parameters
|
|
|
- /// - `f: &mut fmt::Formatter<'_>`: The formatter to write the string representation.
|
|
|
- ///
|
|
|
- /// # Returns
|
|
|
- /// A `fmt::Result` indicating success or failure of the formatting operation.
|
|
|
- fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
|
|
|
- match self {
|
|
|
- BinOutlier::SA => write!(f, "SA"),
|
|
|
- BinOutlier::Start => write!(f, "Start"),
|
|
|
- BinOutlier::End => write!(f, "End"),
|
|
|
- }
|
|
|
- }
|
|
|
-}
|
|
|
-
|
|
|
pub struct CountsReader<R> {
|
|
|
reader: TabixReader<R>,
|
|
|
}
|
|
|
@@ -460,8 +500,8 @@ struct RepRec {
|
|
|
}
|
|
|
|
|
|
struct BinState {
|
|
|
- // alignment-based per-base arrays (bitwise identical to current)
|
|
|
- depths: Vec<u32>,
|
|
|
+ depths_fwd: Vec<u32>, // sens (forward strand)
|
|
|
+ depths_rev: Vec<u32>, // antisens (reverse strand)
|
|
|
lowq: Vec<u32>,
|
|
|
fb_invs: Vec<u32>,
|
|
|
|
|
|
@@ -472,7 +512,8 @@ struct BinState {
|
|
|
impl BinState {
|
|
|
fn new(len: usize) -> Self {
|
|
|
Self {
|
|
|
- depths: vec![0; len],
|
|
|
+ depths_fwd: vec![0; len],
|
|
|
+ depths_rev: vec![0; len],
|
|
|
lowq: vec![0; len],
|
|
|
fb_invs: vec![0; len],
|
|
|
reps: HashMap::new(),
|
|
|
@@ -587,6 +628,7 @@ pub fn scan_contig(
|
|
|
}
|
|
|
|
|
|
let is_lowq = rec.mapq() < min_mapq;
|
|
|
+ let is_reverse = rec.is_reverse();
|
|
|
let has_sa = matches!(rec.aux(b"SA"), Ok(rust_htslib::bam::record::Aux::String(_)));
|
|
|
|
|
|
// Representative fields (same ones used by old count_reads_sa_start_end on reads_store.values())
|
|
|
@@ -626,10 +668,14 @@ pub fn scan_contig(
|
|
|
let off0 = (ov_start - def.start) as usize;
|
|
|
let off1 = (ov_end - def.start) as usize;
|
|
|
for i in off0..=off1 {
|
|
|
- b.depths[i] += 1;
|
|
|
if is_lowq {
|
|
|
b.lowq[i] += 1;
|
|
|
}
|
|
|
+ if is_reverse {
|
|
|
+ b.depths_rev[i] += 1;
|
|
|
+ } else {
|
|
|
+ b.depths_fwd[i] += 1;
|
|
|
+ }
|
|
|
}
|
|
|
|
|
|
for pos in fb_positions.into_iter().flatten() {
|
|
|
@@ -655,16 +701,11 @@ pub fn scan_contig(
|
|
|
for (def_idx, def) in defs.into_iter().enumerate() {
|
|
|
let b = &mut bins[def_idx];
|
|
|
|
|
|
- // n_reads is unique-qname count (same as old reads_store.len())
|
|
|
let n_reads = b.reps.len() as u32;
|
|
|
-
|
|
|
- // SA count from representative records (same as old reads_store.values())
|
|
|
let n_sa = b.reps.values().filter(|r| r.has_sa).count() as u32;
|
|
|
|
|
|
- // max start/end pileups from representative records (same as old HashMap counting)
|
|
|
let mut start_counts: HashMap<u32, u32> = HashMap::new();
|
|
|
let mut end_counts: HashMap<u32, u32> = HashMap::new();
|
|
|
-
|
|
|
for r in b.reps.values() {
|
|
|
if (def.start..=def.end).contains(&r.aln_start) {
|
|
|
*start_counts.entry(r.aln_start).or_insert(0) += 1;
|
|
|
@@ -673,15 +714,31 @@ pub fn scan_contig(
|
|
|
*end_counts.entry(r.aln_end).or_insert(0) += 1;
|
|
|
}
|
|
|
}
|
|
|
-
|
|
|
let max_start = start_counts.values().copied().max().unwrap_or(0);
|
|
|
let max_end = end_counts.values().copied().max().unwrap_or(0);
|
|
|
|
|
|
- // coverage from alignment-based depths vector (same as old mean_coverage_from_depths)
|
|
|
- let coverage = if b.depths.is_empty() {
|
|
|
+ // total depths + per-position sens ratio, reconstructed from fwd/rev
|
|
|
+ let mut depths: Vec<u32> = Vec::with_capacity(b.depths_fwd.len());
|
|
|
+ let mut depths_sens_ratio: Vec<f64> = Vec::with_capacity(b.depths_fwd.len());
|
|
|
+ let mut sum_fwd: u64 = 0;
|
|
|
+ let mut sum_rev: u64 = 0;
|
|
|
+ for (&fwd, &rev) in b.depths_fwd.iter().zip(b.depths_rev.iter()) {
|
|
|
+ let total = fwd + rev;
|
|
|
+ depths.push(total);
|
|
|
+ depths_sens_ratio.push(if total == 0 {
|
|
|
+ f64::NAN
|
|
|
+ } else {
|
|
|
+ fwd as f64 / total as f64
|
|
|
+ });
|
|
|
+ sum_fwd += fwd as u64;
|
|
|
+ sum_rev += rev as u64;
|
|
|
+ }
|
|
|
+ let sum_depths = sum_fwd + sum_rev;
|
|
|
+
|
|
|
+ let coverage = if depths.is_empty() {
|
|
|
0.0
|
|
|
} else {
|
|
|
- b.depths.iter().sum::<u32>() as f64 / (b.depths.len() as f64)
|
|
|
+ sum_depths as f64 / (depths.len() as f64)
|
|
|
};
|
|
|
|
|
|
let n_reads_f = n_reads as f64;
|
|
|
@@ -695,6 +752,37 @@ pub fn scan_contig(
|
|
|
)
|
|
|
};
|
|
|
|
|
|
+ // NEW 1A: ratio_lowq = sum(low_qualities) / sum(depths)
|
|
|
+ let sum_lowq: u64 = b.lowq.iter().map(|&v| v as u64).sum();
|
|
|
+ let ratio_lowq = if sum_depths == 0 {
|
|
|
+ f64::NAN
|
|
|
+ } else {
|
|
|
+ sum_lowq as f64 / sum_depths as f64
|
|
|
+ };
|
|
|
+
|
|
|
+ // NEW 2C (depth-normalized): ratio_fb = max(fb_invs[i]) / depths[i]
|
|
|
+ let ratio_fb = b
|
|
|
+ .fb_invs
|
|
|
+ .iter()
|
|
|
+ .enumerate()
|
|
|
+ .max_by_key(|&(_, v)| *v)
|
|
|
+ .map(|(idx, &max_val)| {
|
|
|
+ let d = depths[idx];
|
|
|
+ if d == 0 {
|
|
|
+ f64::NAN
|
|
|
+ } else {
|
|
|
+ max_val as f64 / d as f64
|
|
|
+ }
|
|
|
+ })
|
|
|
+ .unwrap_or(f64::NAN);
|
|
|
+
|
|
|
+ // NEW 3A: ratio_sens = pooled sum(fwd) / sum(fwd+rev), bin-level
|
|
|
+ let ratio_sens = if sum_depths == 0 {
|
|
|
+ f64::NAN
|
|
|
+ } else {
|
|
|
+ sum_fwd as f64 / sum_depths as f64
|
|
|
+ };
|
|
|
+
|
|
|
out.push(BinCount {
|
|
|
contig: contig.to_string(),
|
|
|
start: def.start,
|
|
|
@@ -704,10 +792,14 @@ pub fn scan_contig(
|
|
|
ratio_sa,
|
|
|
ratio_start,
|
|
|
ratio_end,
|
|
|
+ ratio_lowq,
|
|
|
+ ratio_fb,
|
|
|
+ ratio_sens,
|
|
|
outlier: None,
|
|
|
- depths: std::mem::take(&mut b.depths),
|
|
|
+ depths,
|
|
|
low_qualities: std::mem::take(&mut b.lowq),
|
|
|
fb_invs: std::mem::take(&mut b.fb_invs),
|
|
|
+ depths_sens_ratio,
|
|
|
});
|
|
|
}
|
|
|
|
|
|
@@ -838,8 +930,8 @@ pub fn par_whole_scan(id: &str, time_point: &str, config: &Config) -> anyhow::Re
|
|
|
|
|
|
let nf = row.split('\t').count();
|
|
|
anyhow::ensure!(
|
|
|
- nf == 12,
|
|
|
- "BUG: to_tsv_row() produced {nf} fields instead of 12 for {contig}: {}:{}-{}\nrow={row}",
|
|
|
+ nf == 16,
|
|
|
+ "BUG: to_tsv_row() produced {nf} fields instead of 16 for {contig}: {}:{}-{}\nrow={row}",
|
|
|
bin.contig,
|
|
|
bin.start,
|
|
|
bin.end
|
|
|
@@ -924,8 +1016,8 @@ fn validate_count_file(
|
|
|
let fields = line.split_fields();
|
|
|
|
|
|
anyhow::ensure!(
|
|
|
- fields.len() == 12,
|
|
|
- "{path} line {n}: expected 12 fields, got {}",
|
|
|
+ fields.len() == 16,
|
|
|
+ "{path} line {n}: expected 16 fields, got {}",
|
|
|
fields.len()
|
|
|
);
|
|
|
anyhow::ensure!(
|
|
|
@@ -943,6 +1035,43 @@ fn validate_count_file(
|
|
|
Ok(())
|
|
|
}
|
|
|
|
|
|
+/// Represents types of outliers that can be detected in a genomic bin.
|
|
|
+#[derive(Debug, Clone)]
|
|
|
+pub enum BinOutlier {
|
|
|
+ /// Indicates an outlier in supplementary alignments.
|
|
|
+ SA,
|
|
|
+ /// Indicates an outlier in reads starting in this bin.
|
|
|
+ Start,
|
|
|
+ /// Indicates an outlier in reads ending in this bin.
|
|
|
+ End,
|
|
|
+ /// low_qualities ratio outlier
|
|
|
+ LowQ,
|
|
|
+ /// fb_invs ratio outlier
|
|
|
+ Fb,
|
|
|
+ /// Strand (sens/antisens) bias outlier
|
|
|
+ Sens,
|
|
|
+}
|
|
|
+
|
|
|
+impl fmt::Display for BinOutlier {
|
|
|
+ /// Implements the `Display` trait for `BinOutlier`.
|
|
|
+ ///
|
|
|
+ /// # Parameters
|
|
|
+ /// - `f: &mut fmt::Formatter<'_>`: The formatter to write the string representation.
|
|
|
+ ///
|
|
|
+ /// # Returns
|
|
|
+ /// A `fmt::Result` indicating success or failure of the formatting operation.
|
|
|
+ fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
|
|
|
+ match self {
|
|
|
+ BinOutlier::SA => write!(f, "SA"),
|
|
|
+ BinOutlier::Start => write!(f, "Start"),
|
|
|
+ BinOutlier::End => write!(f, "End"),
|
|
|
+ BinOutlier::LowQ => write!(f, "LowQ"),
|
|
|
+ BinOutlier::Fb => write!(f, "Fb"),
|
|
|
+ BinOutlier::Sens => write!(f, "Sens"),
|
|
|
+ }
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
/// Identifies and marks outliers in a slice of `BinCount` objects based on various ratio metrics.
|
|
|
///
|
|
|
/// # Parameters
|
|
|
@@ -986,11 +1115,28 @@ pub fn fill_outliers(bin_counts: &mut [BinCount]) {
|
|
|
let intermediate: Vec<_> = bin_counts
|
|
|
.par_iter()
|
|
|
.enumerate()
|
|
|
- .map(|(i, c)| (i, [c.ratio_sa, c.ratio_start, c.ratio_end]))
|
|
|
+ .map(|(i, c)| {
|
|
|
+ (
|
|
|
+ i,
|
|
|
+ [
|
|
|
+ c.ratio_sa,
|
|
|
+ c.ratio_start,
|
|
|
+ c.ratio_end,
|
|
|
+ c.ratio_lowq,
|
|
|
+ c.ratio_fb,
|
|
|
+ c.ratio_sens,
|
|
|
+ ],
|
|
|
+ )
|
|
|
+ })
|
|
|
.collect();
|
|
|
-
|
|
|
- let outlier_types = [BinOutlier::SA, BinOutlier::Start, BinOutlier::End];
|
|
|
-
|
|
|
+ let outlier_types = [
|
|
|
+ BinOutlier::SA,
|
|
|
+ BinOutlier::Start,
|
|
|
+ BinOutlier::End,
|
|
|
+ BinOutlier::LowQ,
|
|
|
+ BinOutlier::Fb,
|
|
|
+ BinOutlier::Sens,
|
|
|
+ ];
|
|
|
let outliers: Vec<(usize, BinOutlier)> = outlier_types
|
|
|
.iter()
|
|
|
.enumerate()
|
|
|
@@ -1002,13 +1148,11 @@ pub fn fill_outliers(bin_counts: &mut [BinCount]) {
|
|
|
(!ratio.is_nan()).then_some((*i, ratio))
|
|
|
})
|
|
|
.unzip();
|
|
|
-
|
|
|
filter_outliers_modified_z_score_with_indices(&ratios, indices)
|
|
|
.into_iter()
|
|
|
.map(move |i| (i, outlier_type.clone()))
|
|
|
})
|
|
|
.collect();
|
|
|
-
|
|
|
outliers.iter().for_each(|(i, outlier_type)| {
|
|
|
bin_counts[*i]
|
|
|
.outlier
|