diff --git a/Cargo.lock b/Cargo.lock index 49d33b3..9f65421 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -329,7 +329,7 @@ checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" [[package]] name = "recmap" -version = "0.3.4" +version = "0.3.5" dependencies = [ "clap", "csv", diff --git a/Cargo.toml b/Cargo.toml index d29322a..be68999 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "recmap" -version = "0.3.4" +version = "0.3.5" edition = "2021" license = "MIT" authors = ["Vince Buffalo "] diff --git a/src/numeric.rs b/src/numeric.rs index 4ec9b5c..d13f748 100644 --- a/src/numeric.rs +++ b/src/numeric.rs @@ -82,7 +82,7 @@ pub fn recomb_dist_matrix( positions_x: &[T], positions_y: &[T], haldane: bool, - min_rec: Option, + rec_floor: Option, ) -> Array2 { let mut dist_matrix = Array2::::zeros((positions_x.len(), positions_y.len())); for (i, &pos_x) in positions_x.iter().enumerate() { @@ -93,7 +93,7 @@ pub fn recomb_dist_matrix( } else { haldane_inverse(dist) }; - if let Some(min_rec_rate) = min_rec { + if let Some(min_rec_rate) = rec_floor { rf = rf.max(min_rec_rate); } dist_matrix[[i, j]] = rf; @@ -164,7 +164,7 @@ where let y2 = ToPrimitive::to_f64(&y[idx])?; let x0 = ToPrimitive::to_f64(&x0)?; - // Perform linear interpolation + // linear interpolation let y0 = y1 + (y2 - y1) * (x0 - x1) / (x2 - x1); NumCast::from(y0) diff --git a/src/recmap.rs b/src/recmap.rs index 95e69e4..753946c 100644 --- a/src/recmap.rs +++ b/src/recmap.rs @@ -176,14 +176,13 @@ impl RateMap { self.rates .iter() .zip(self.span().iter()) - .map(|(&x, &y)| x * y as RateFloat) + .map(|(&rate, &span)| rate * span as RateFloat) .collect() } - /// Calculate the cumulative map length at each position. + /// Calculate the cumulative map length in Morgans at each position. pub fn calc_cumulative_mass(&mut self) { let mass = self.mass(); - //println!("mass: {:?}", mass); let cumulative_sum: Vec<_> = mass .iter() .scan(0.0, |state, &x| { @@ -358,7 +357,8 @@ impl RecMap { new_rate_map.ends.push(end); new_rate_map.rates.push(rate); - // if there is a fourth column, let's use it for validation + // if there is a fourth column, we could use it for validation + // TODO if has_fourth_column { //new_rate_map.calc_cumulative_mass(); //assert_floats_eq(&new_rate_map.map_pos, map_positions.as_slice(), 0.01); @@ -392,6 +392,8 @@ impl RecMap { map: rec_map, metadata, }; + // generate the map positions from the marker positions + // and the per-basepair rates.. rec_map.generate_map_positions(); Ok(rec_map) } @@ -434,7 +436,9 @@ impl RecMap { .get(name) .ok_or(RecMapError::NoChrom(name.to_string()))?; let ends = &rate_map.ends; - let interp_result = interp1d(&ends[0..ends.len() - 1], &rate_map.map_pos, position); + let interp_result = interp1d(&ends[0..ends.len() - 1], + &rate_map.map_pos, + position); let interpolated_map_pos = interp_result.ok_or(RecMapError::LookupOutOfBounds(name.to_string(), position))?; Ok(interpolated_map_pos) @@ -480,11 +484,11 @@ impl RecMap { positions_x: &[Position], positions_y: &[Position], haldane: bool, - min_rec: Option, + rec_floor: Option, ) -> Result, RecMapError> { let x_pos = self.interpolate_map_positions(chrom, positions_x)?; let y_pos = self.interpolate_map_positions(chrom, positions_y)?; - Ok(recomb_dist_matrix(&x_pos, &y_pos, haldane, min_rec)) + Ok(recomb_dist_matrix(&x_pos, &y_pos, haldane, rec_floor)) } /// Write recombination map to HapMap-formatted file. @@ -665,7 +669,7 @@ mod tests { &chr1_map.map_pos, &to_morgans(vec![0.0, 5.5, 13., 34.65]), 1e-3, - ); + ); } #[test]