Add an STRC regression test using public HG00733 HiFi data

tests/phaser_tests.rs

@@ -30,6 +30,33 @@ fn slice_matches_with_gaps(x: &[impl AsRef<[u8]>], y: &[impl AsRef<[u8]>]) -> bo
             .all(|z| x.iter().any(|w| w.as_ref().same_without_dels(&z.as_ref())))
 }
 
+fn strc_haplotype_class(name: &str) -> &'static str {
+    let identity = name
+        .rsplit('_')
+        .next()
+        .expect("STRC haplotype name should contain an identity label");
+    let has_numeric_suffix = |prefix: &str| {
+        identity.strip_prefix(prefix).is_some_and(|suffix| {
+            !suffix.is_empty() && suffix.chars().all(|digit| digit.is_ascii_digit())
+        })
+    };
+
+    if has_numeric_suffix("strchap") {
+        "strc"
+    } else if has_numeric_suffix("strcp1hap") {
+        "strcp1"
+    } else {
+        panic!("Unexpected STRC haplotype identity label in {name}");
+    }
+}
+
+#[test]
+fn strc_haplotype_class_is_closed_over_known_identities() {
+    assert_eq!(strc_haplotype_class("strc_strchap1"), "strc");
+    assert_eq!(strc_haplotype_class("strc_strcp1hap1"), "strcp1");
+    assert!(std::panic::catch_unwind(|| strc_haplotype_class("strc_otherhap1")).is_err());
+}
+
 #[test]
 fn smn1_ok() -> DResult {
     util::init_log(log::LevelFilter::Info);
@@ -104,6 +131,213 @@ fn smn1_ok() -> DResult {
     Ok(())
 }
 
+// Same-input oracle: Paraphase v3.5.0 c8016dff40d105501868719bdae0defd9937abfe
+// on the public HPRC HG00733 GRCh38 HiFi locus fixture (m54329U_2019).
+// Genome depth 31 comes from the same source BAM's bundled depth-probe union.
+#[test]
+fn strc_ok() -> DResult {
+    util::init_log(log::LevelFilter::Info);
+
+    let outdir = tempfile::TempDir::new()?;
+    let genome_bam = test_file("bams/HG00733.strc.bam");
+    let genome_path = if let Ok(x) = std::env::var("HG38") {
+        x.trim_end_matches(".mmi").to_string()
+    } else {
+        log::warn!("Skipping test: HG38 env is not set.");
+        return Ok(());
+    };
+    let gene_name = "strc";
+    let depth = depth::Result {
+        median: 31.0,
+        median_absolute_difference: 0.1,
+        sex: depth::Sex::Other,
+    };
+    let region_config = config::Region::try_load(None)?;
+    let settings = phaser::Settings::new(
+        "HG00733",
+        (genome_path, genome_bam),
+        outdir.path(),
+        gene_name,
+        &region_config,
+        /* genome depth= */ Some(depth),
+        /* sex = */ None,
+        String::from("38"),
+        None,
+        0.03,
+        false,
+    );
+
+    let gene_config = config::Gene::try_load(None)?;
+    let mut phaser = Phaser::new(
+        settings,
+        Some(gene_config),
+        None, // Option<SiteSelectionSettings>
+        None, // Option<RealignSettings>
+    )?;
+    let call = phaser.run()?;
+    assert_eq!(phaser.region_avg_depth[0], (71.0f32, 82.0f32));
+
+    let json_to_match = test_file("jsons/HG00733.strc.v3.5.0.json.xz");
+    let json_data = ParsedParaphaseOutputJSON::from_path(&json_to_match, Some(&region_config))?;
+    let expected_call = json_data
+        .object
+        .get(gene_name)
+        .and_then(serde_json::Value::as_object)
+        .ok_or("Missing STRC call in expected json")?;
+    let gene_data = json_data
+        .gene_data
+        .get(gene_name)
+        .ok_or("Missing STRC gene data in expected json")?;
+
+    assert_eq!(call.total_cn, Some(4));
+    assert_eq!(
+        call.region_specific_info
+            .get("gene_cn")
+            .and_then(serde_json::Value::as_i64),
+        Some(2)
+    );
+    assert_eq!(
+        call.region_specific_info
+            .get("intergenic_depth")
+            .and_then(serde_json::Value::as_f64),
+        Some(31.0)
+    );
+
+    let expected_final_haplotypes = serde_json::from_value::<BTreeMap<String, String>>(
+        expected_call
+            .get("final_haplotypes")
+            .ok_or("Missing STRC final_haplotypes in expected json")?
+            .clone(),
+    )?;
+    let final_classes_by_sequence = |haplotypes: &BTreeMap<String, String>| {
+        haplotypes
+            .iter()
+            .map(|(sequence, name)| (sequence.clone(), strc_haplotype_class(name)))
+            .collect::<BTreeMap<_, _>>()
+    };
+    assert_eq!(
+        final_classes_by_sequence(&call.final_haplotypes),
+        final_classes_by_sequence(&expected_final_haplotypes)
+    );
+
+    let expected_two_copy_haplotypes = serde_json::from_value::<Vec<String>>(
+        expected_call
+            .get("two_copy_haplotypes")
+            .ok_or("Missing STRC two_copy_haplotypes in expected json")?
+            .clone(),
+    )?;
+    let two_copy_classes = |haplotypes: &[String]| {
+        haplotypes
+            .iter()
+            .map(|name| strc_haplotype_class(name))
+            .sorted()
+            .collect::<Vec<_>>()
+    };
+    assert_eq!(
+        two_copy_classes(&call.two_copy_haplotypes),
+        two_copy_classes(&expected_two_copy_haplotypes)
+    );
+
+    let expected_sites = gene_data
+        .sites_for_phasing
+        .iter()
+        .map(std::string::ToString::to_string)
+        .sorted()
+        .collect::<Vec<_>>();
+    let found_sites = call
+        .sites_for_phasing
+        .iter()
+        .cloned()
+        .sorted()
+        .collect::<Vec<_>>();
+    assert_eq!(expected_sites, found_sites);
+
+    let marker_index = |site: &str| {
+        let indices = call
+            .sites_for_phasing
+            .iter()
+            .enumerate()
+            .filter_map(|(index, found_site)| (found_site == site).then_some(index))
+            .collect::<Vec<_>>();
+        assert_eq!(
+            indices.len(),
+            1,
+            "Expected exactly one {site} site for STRC identity markers"
+        );
+        indices[0]
+    };
+    let strc_snv_index = marker_index("43602487_C_G");
+    let strc_deletion_index = marker_index("43602630_del_314");
+    for (sequence, name) in &call.final_haplotypes {
+        let marker = |index: usize| {
+            sequence
+                .as_bytes()
+                .get(index)
+                .copied()
+                .map(char::from)
+                .unwrap_or_else(|| panic!("Missing STRC identity marker in haplotype {sequence}"))
+        };
+        let found_marker_pair = (marker(strc_snv_index), marker(strc_deletion_index));
+        let expected_marker_pair = match strc_haplotype_class(name) {
+            "strc" => ('1', '1'),
+            "strcp1" => ('2', '3'),
+            identity => unreachable!("Unhandled STRC haplotype identity {identity}"),
+        };
+        assert_eq!(
+            found_marker_pair, expected_marker_pair,
+            "Unexpected identity markers for {name} on haplotype {sequence}"
+        );
+    }
+
+    let expected_haplotypes = gene_data
+        .assembled_haps
+        .as_ref()
+        .ok_or("Missing STRC assembled_haplotypes in expected json")?
+        .iter()
+        .map(std::string::ToString::to_string)
+        .sorted()
+        .collect::<Vec<_>>();
+    let found_haplotypes = call
+        .assembled_haplotypes
+        .iter()
+        .cloned()
+        .sorted()
+        .collect::<Vec<_>>();
+    assert_eq!(expected_haplotypes, found_haplotypes);
+
+    let expected_haplotype_details =
+        serde_json::from_value::<BTreeMap<String, paraphase::phaser::HapInfoForJson>>(
+            expected_call
+                .get("haplotype_details")
+                .ok_or("Missing STRC haplotype_details in expected json")?
+                .clone(),
+        )?;
+    let variants_by_sequence =
+        |haplotypes: &BTreeMap<String, String>,
+         details: &BTreeMap<String, paraphase::phaser::HapInfoForJson>| {
+            haplotypes
+                .iter()
+                .map(|(sequence, name)| {
+                    let variants = details
+                        .get(name)
+                        .expect("Missing STRC haplotype details")
+                        .variants
+                        .iter()
+                        .cloned()
+                        .sorted()
+                        .collect::<Vec<_>>();
+                    (sequence.clone(), variants)
+                })
+                .collect::<BTreeMap<_, _>>()
+        };
+    assert_eq!(
+        variants_by_sequence(&call.final_haplotypes, &call.haplotype_details),
+        variants_by_sequence(&expected_final_haplotypes, &expected_haplotype_details)
+    );
+
+    Ok(())
+}
+
 #[test]
 fn amy1_ok() -> DResult {
     util::init_log(log::LevelFilter::Info);