check coverage when marking initial solid nodes

src/lib.rs

@@ -410,6 +410,23 @@ pub fn run_trio_analysis(settings: &TrioSettings) -> Result<(), Box<dyn Error>>
     );
     let trio_infos = trio::read_trio(&settings.markers)?;
 
+    let solid_cov_est = weighted_mean_solid_cov(&g, settings.solid_len);
+    if settings.suspect_homozygous_cov_coeff > 0. || settings.solid_homozygous_cov_coeff > 0. {
+        info!("Coverage estimate based on long nodes was {solid_cov_est}");
+        if solid_cov_est == 0. {
+            warn!("Looks like the graph didn't have coverage information, which we were hoping to use. \
+                    Consider providing it or changing --suspect-homozygous-cov-coeff and --solid-homozygous-cov-coeff");
+        }
+    }
+
+    let suspect_homozygous_cov = if settings.suspect_homozygous_cov_coeff < 0. {
+        None
+    } else {
+        Some(settings.suspect_homozygous_cov_coeff * solid_cov_est)
+    };
+
+    let solid_homozygous_cov = settings.solid_homozygous_cov_coeff * solid_cov_est;
+
     info!("Assigning initial parental groups to the nodes");
     let assignments = trio::assign_parental_groups(
         &g,
@@ -425,6 +442,7 @@ pub fn run_trio_analysis(settings: &TrioSettings) -> Result<(), Box<dyn Error>>
             issue_ratio: settings.issue_ratio.unwrap_or(settings.marker_ratio),
         },
         settings.solid_len,
+        solid_homozygous_cov,
     );
 
     let raw_cnts = trio_infos
@@ -440,23 +458,6 @@ pub fn run_trio_analysis(settings: &TrioSettings) -> Result<(), Box<dyn Error>>
         output_coloring(&g, &assignments, output, &hap_names)?;
     }
 
-    let solid_cov_est = weighted_mean_solid_cov(&g, settings.solid_len);
-    if settings.suspect_homozygous_cov_coeff > 0. || settings.solid_homozygous_cov_coeff > 0. {
-        info!("Coverage estimate based on long nodes was {solid_cov_est}");
-        if solid_cov_est == 0. {
-            warn!("Looks like the graph didn't have coverage information, which we were hoping to use. \
-                    Consider providing it or changing --suspect-homozygous-cov-coeff and --solid-homozygous-cov-coeff");
-        }
-    }
-
-    let suspect_homozygous_cov = if settings.suspect_homozygous_cov_coeff < 0. {
-        None
-    } else {
-        Some(settings.suspect_homozygous_cov_coeff * solid_cov_est)
-    };
-
-    let solid_homozygous_cov = settings.solid_homozygous_cov_coeff * solid_cov_est;
-
     info!("Marking homozygous nodes");
     let assigner = trio::HomozygousAssigner::new(
         &g,

src/trio.rs

@@ -220,6 +220,7 @@ pub fn assign_parental_groups(
     trio_infos: &[TrioInfo],
     settings: &GroupAssignmentSettings,
     solid_len: usize,
+    solid_cov: f64,
 ) -> AssignmentStorage {
     let mut assignments = AssignmentStorage::new();
 
@@ -230,13 +231,13 @@ pub fn assign_parental_groups(
             settings.issue_cnt, settings.issue_sparsity, settings.issue_ratio);
     assert!(settings.issue_ratio <= settings.assign_ratio);
 
-    let assign_node_f = |x: usize, y: usize, node_len: usize| {
+    let assign_node_f = |x: usize, y: usize, node_len: usize, node_cov: f64| {
         assert!(x >= y);
         let tot = x + y;
         tot >= settings.assign_cnt
             && node_len <= tot * settings.assign_sparsity
             && ((x as f64) > settings.assign_ratio * (y as f64) - 1e-6
-                || (node_len > solid_len && (x as f64) > settings.solid_ratio * (y as f64) - 1e-6))
+                || (node_len > solid_len && (x as f64) > settings.solid_ratio * (y as f64) - 1e-6 && node_cov < solid_cov - 1e-6))
     };
 
     let issue_node_f = |x: usize, y: usize, node_len: usize| {
@@ -251,6 +252,7 @@ pub fn assign_parental_groups(
     for trio_info in trio_infos {
         let node_id = g.name2id(&trio_info.node_name);
         let node_len = g.node_length(node_id);
+        let node_cov = g.node(node_id).coverage;
         debug!(
             "Looking at node {} (len={}), mat:pat={}",
             trio_info.node_name,
@@ -269,6 +271,7 @@ pub fn assign_parental_groups(
             max(trio_info.mat, trio_info.pat),
             min(trio_info.mat, trio_info.pat),
             node_len,
+            node_cov,
         ) {
             if trio_info.mat >= trio_info.pat {
                 debug!("Looks MATERNAL");