diff --git a/tests/data/bams/HG00733.strc.bam b/tests/data/bams/HG00733.strc.bam new file mode 100644 index 0000000..8cc4d04 Binary files /dev/null and b/tests/data/bams/HG00733.strc.bam differ diff --git a/tests/data/bams/HG00733.strc.bam.bai b/tests/data/bams/HG00733.strc.bam.bai new file mode 100644 index 0000000..92f1be9 Binary files /dev/null and b/tests/data/bams/HG00733.strc.bam.bai differ diff --git a/tests/data/jsons/HG00733.strc.v3.5.0.json.xz b/tests/data/jsons/HG00733.strc.v3.5.0.json.xz new file mode 100644 index 0000000..70cfc77 Binary files /dev/null and b/tests/data/jsons/HG00733.strc.v3.5.0.json.xz differ diff --git a/tests/phaser_tests.rs b/tests/phaser_tests.rs index 6532a91..7269ea6 100644 --- a/tests/phaser_tests.rs +++ b/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, + ®ion_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 + None, // Option + )?; + 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(®ion_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::>( + expected_call + .get("final_haplotypes") + .ok_or("Missing STRC final_haplotypes in expected json")? + .clone(), + )?; + let final_classes_by_sequence = |haplotypes: &BTreeMap| { + haplotypes + .iter() + .map(|(sequence, name)| (sequence.clone(), strc_haplotype_class(name))) + .collect::>() + }; + 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::>( + 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::>() + }; + 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::>(); + let found_sites = call + .sites_for_phasing + .iter() + .cloned() + .sorted() + .collect::>(); + 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::>(); + 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::>(); + let found_haplotypes = call + .assembled_haplotypes + .iter() + .cloned() + .sorted() + .collect::>(); + assert_eq!(expected_haplotypes, found_haplotypes); + + let expected_haplotype_details = + serde_json::from_value::>( + expected_call + .get("haplotype_details") + .ok_or("Missing STRC haplotype_details in expected json")? + .clone(), + )?; + let variants_by_sequence = + |haplotypes: &BTreeMap, + details: &BTreeMap| { + haplotypes + .iter() + .map(|(sequence, name)| { + let variants = details + .get(name) + .expect("Missing STRC haplotype details") + .variants + .iter() + .cloned() + .sorted() + .collect::>(); + (sequence.clone(), variants) + }) + .collect::>() + }; + 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);