diff --git a/.claude/settings.json b/.claude/settings.json deleted file mode 100644 index e6e2674b..00000000 --- a/.claude/settings.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "env": { - "ECC_DISABLED_HOOKS": "pre:bash:gateguard-fact-force,pre:edit-write:gateguard-fact-force" - } -} diff --git a/.gitignore b/.gitignore index 35fe5f24..e590c35d 100644 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,6 @@ testdata/reference/large/ .DS_Store validation/work/ validation/output/ + +# Local Claude Code settings (do not commit) +.claude/ diff --git a/CLAUDE.md b/CLAUDE.md index da221cb6..2b4b5218 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -23,7 +23,7 @@ Running log: `PORT_LOG.md`. ## Working rules -1. **Test before commit.** Every commit must leave the build green: `swift build && swift test` in `tools/conversion/` for Phase 0 work; `cmake --build build && ctest -V` for Phases 1+. +1. **Test before commit.** Every commit must leave the build green: `swift build && swift test` in `tools/conversion/` for Phase 0 work; `cmake --build build-macos && ctest -V` for Phases 1+. (Use `build-macos`, not `build` — a `build/` dir collides with the Bazel `BUILD` file on macOS's case-insensitive filesystem.) 2. **Never degrade scientific precision.** F1 thresholds are gates, not goals. If we slip below, we fix the root cause — we do not lower the bar. 3. **Never bypass an error.** No `--no-verify`, no swallowed exceptions, no commenting out of failing tests. Diagnose the root cause. 4. **Document every critical decision** in `PORT_LOG.md` with date, context, alternatives considered, and rationale. diff --git a/cmake/deps.cmake b/cmake/deps.cmake index 0fa86cbc..95c485f3 100644 --- a/cmake/deps.cmake +++ b/cmake/deps.cmake @@ -34,9 +34,18 @@ set_target_properties(htslib::htslib PROPERTIES IMPORTED_LOCATION "${HTSLIB_LIB}" INTERFACE_INCLUDE_DIRECTORIES "${HTSLIB_PREFIX}/include" ) +# Resolve libdeflate via Homebrew rather than a hardcoded /opt/homebrew path, +# so the build works under a non-default Homebrew prefix or keg-only layout. +execute_process( + COMMAND ${BREW_EXECUTABLE} --prefix libdeflate + OUTPUT_VARIABLE LIBDEFLATE_PREFIX + OUTPUT_STRIP_TRAILING_WHITESPACE +) +find_library(LIBDEFLATE_LIB NAMES libdeflate.a deflate + PATHS "${LIBDEFLATE_PREFIX}/lib" REQUIRED) target_link_libraries(htslib::htslib INTERFACE "-framework CoreFoundation" - /opt/homebrew/lib/libdeflate.a + "${LIBDEFLATE_LIB}" z bz2 lzma curl ) message(STATUS "htslib: ${HTSLIB_LIB}") @@ -83,10 +92,10 @@ FetchContent_Declare( URL https://github.com/mengyao/Complete-Striped-Smith-Waterman-Library/archive/v1.2.5.tar.gz URL_HASH SHA256=b294c0cb6f0f3d578db11b4112a88b20583b9d4190b0a9cf04d83bb6a8704d9a ) -FetchContent_GetProperties(libssw) -if(NOT libssw_POPULATED) - FetchContent_Populate(libssw) -endif() +# libssw ships no CMakeLists, so MakeAvailable just populates ${libssw_SOURCE_DIR} +# (no add_subdirectory) — and avoids the deprecated single-arg +# FetchContent_Populate(libssw) call. +FetchContent_MakeAvailable(libssw) # OVERLAY: replace the vendored sse2neon.h (Ratcliff/NVIDIA early version, # 8798 lines, missing fixes) with the modern DLTcollab fork (11744 lines, diff --git a/cmake/protos.cmake b/cmake/protos.cmake index 3496f5b5..0a4e43b8 100644 --- a/cmake/protos.cmake +++ b/cmake/protos.cmake @@ -46,7 +46,7 @@ endfunction() # 1. nucleus protos (self-contained, no TF imports) # --------------------------------------------------------------------------- set(NUCLEUS_PROTO_ROOT "${CMAKE_SOURCE_DIR}/third_party/nucleus/protos") -file(GLOB NUCLEUS_PROTOS "${NUCLEUS_PROTO_ROOT}/*.proto") +file(GLOB NUCLEUS_PROTOS CONFIGURE_DEPENDS "${NUCLEUS_PROTO_ROOT}/*.proto") set(NUCLEUS_PB_SRCS) foreach(_p ${NUCLEUS_PROTOS}) @@ -72,7 +72,7 @@ add_library(proto_tf_example ALIAS proto_nucleus) # 2. deepvariant protos # --------------------------------------------------------------------------- set(DV_PROTO_ROOT "${CMAKE_SOURCE_DIR}/deepvariant/protos") -file(GLOB DV_PROTOS "${DV_PROTO_ROOT}/*.proto") +file(GLOB DV_PROTOS CONFIGURE_DEPENDS "${DV_PROTO_ROOT}/*.proto") set(DV_PB_SRCS) foreach(_p ${DV_PROTOS}) diff --git a/deepvariant/CMakeLists.txt b/deepvariant/CMakeLists.txt index 4d35576a..ca8bd115 100644 --- a/deepvariant/CMakeLists.txt +++ b/deepvariant/CMakeLists.txt @@ -36,7 +36,7 @@ dv_library(dv_utils # --------------------------------------------------------------------------- # dv_channels — all pileup channel implementations # --------------------------------------------------------------------------- -file(GLOB DV_CHANNEL_SRCS "channels/*.cc") +file(GLOB DV_CHANNEL_SRCS CONFIGURE_DEPENDS "channels/*.cc") list(FILTER DV_CHANNEL_SRCS EXCLUDE REGEX "_test\\.cc$") dv_library(dv_channels SRCS ${DV_CHANNEL_SRCS} diff --git a/deepvariant/native/CMakeLists.txt b/deepvariant/native/CMakeLists.txt index df7f02a9..54704b46 100644 --- a/deepvariant/native/CMakeLists.txt +++ b/deepvariant/native/CMakeLists.txt @@ -82,8 +82,13 @@ target_include_directories(dv_metal_conv_kahan PUBLIC "${CMAKE_SOURCE_DIR}" "${ABSL_PREFIX}/include" ) +# NB: dv_metal_conv_kahan only needs ConvDesc from metal_conv_serial.h (a +# compile-time include available via target_include_directories), not any +# compiled MetalConvSerial symbol — so there is no link edge back to +# dv_metal_conv_serial. Adding one would create a static-library cycle +# (serial -> kahan is the only real edge, via MetalConvKahan::Create) that +# Apple ld reports as "ignoring duplicate libraries". target_link_libraries(dv_metal_conv_kahan PUBLIC - dv_metal_conv_serial # reuse ConvDesc absl::log "-framework Metal" "-framework Foundation" @@ -270,6 +275,24 @@ set(ME_COMPILE_OPTS "-include${CMAKE_SOURCE_DIR}/cmake/tf_stubs/tensorflow/core/platform/tf_compat.h" ) +# --------------------------------------------------------------------------- +# dv_haploid_regions — sex-chromosome haploid-calling helpers shared by +# make_examples (gVCF) and postprocess. The header is dependency-free and +# header-only; only LoadParRegions lives here because it uses zlib to read +# plaintext-or-gzipped PAR BEDs. +# --------------------------------------------------------------------------- +add_library(dv_haploid_regions STATIC + "${CMAKE_CURRENT_SOURCE_DIR}/haploid_regions.cc" +) +target_include_directories(dv_haploid_regions PUBLIC + "${CMAKE_SOURCE_DIR}" + "${ABSL_PREFIX}/include" +) +target_link_libraries(dv_haploid_regions PUBLIC + absl::strings + z # zlib (also pulled in transitively by htslib); for gzopen/gzgets +) + add_library(dv_make_examples_lib STATIC "${CMAKE_CURRENT_SOURCE_DIR}/make_examples_main.cc" "${CMAKE_CURRENT_SOURCE_DIR}/regions.cc" @@ -284,6 +307,7 @@ target_link_libraries(dv_make_examples_lib PUBLIC dv_tfrecord dv_small_model dv_direct_phasing # Phase 9 / Step 4 — per-region read phasing + dv_haploid_regions # haploid gVCF reference confidence (PAR BED loader) realigner proto_dv proto_nucleus @@ -294,6 +318,34 @@ target_link_libraries(dv_make_examples_lib PUBLIC nucleus_io ) +# --------------------------------------------------------------------------- +# dv_methylation_aware_phasing — faithful C++ port of upstream's +# methylation_aware_phasing.{cc,h} (Wilcoxon rank-sum 5mC haplotype voting). +# +# NOTE: this algorithm is compiled and unit-tested here (so it cannot bitrot +# against the native proto definitions), but it is NOT yet invoked by the +# native make_examples pipeline. The native phasing step +# (make_examples_main.cc) uses DirectPhasing only; there is no native +# --enable_methylation_aware_phasing flag, so a user cannot currently engage +# methylation-aware phasing. Wiring PerformMethylationAwarePhasing into the +# make_examples region loop (mirroring make_examples_core.py's methylated-ref- +# site partitioning + post-DirectPhasing call) is deferred follow-up work. +# This target deliberately is not linked into dv_make_examples_lib. +# --------------------------------------------------------------------------- +add_library(dv_methylation_aware_phasing STATIC + "${CMAKE_SOURCE_DIR}/deepvariant/methylation_aware_phasing.cc" +) +target_include_directories(dv_methylation_aware_phasing PUBLIC ${ME_INCLUDE_DIRS}) +target_compile_options(dv_methylation_aware_phasing PRIVATE ${ME_COMPILE_OPTS}) +target_link_libraries(dv_methylation_aware_phasing PUBLIC + proto_dv + proto_nucleus + absl::log + absl::strings + absl::flat_hash_map # raw_hash_set internals referenced by the .cc + absl::flat_hash_set +) + # --------------------------------------------------------------------------- # Phase 3 — postprocess_variants orchestration # --------------------------------------------------------------------------- @@ -305,6 +357,7 @@ target_include_directories(dv_postprocess_lib PUBLIC ${ME_INCLUDE_DIRS}) target_compile_options(dv_postprocess_lib PRIVATE ${ME_COMPILE_OPTS}) target_link_libraries(dv_postprocess_lib PUBLIC dv_tfrecord + dv_haploid_regions # haploid het-zeroing (PAR BED loader) proto_dv proto_nucleus absl::flags @@ -323,8 +376,7 @@ add_executable(debug_metal "${CMAKE_CURRENT_SOURCE_DIR}/debug_metal_main.cc") target_include_directories(debug_metal PRIVATE ${ME_INCLUDE_DIRS}) target_compile_options(debug_metal PRIVATE ${ME_COMPILE_OPTS}) target_link_libraries(debug_metal PRIVATE - dv_metal_inference - dv_weights + dv_metal_inference # transitively links dv_weights (PUBLIC); do not repeat it ) # --------------------------------------------------------------------------- diff --git a/deepvariant/native/call_variants_main.cc b/deepvariant/native/call_variants_main.cc index 6fe9fff8..283f7667 100644 --- a/deepvariant/native/call_variants_main.cc +++ b/deepvariant/native/call_variants_main.cc @@ -751,7 +751,10 @@ int RunCallVariants(int argc, char** argv) { } reader->Close(); - writer->Close(); + if (!writer->Close()) { + LOG(ERROR) << "Failed to flush/close output: " << outfile_path; + return 1; + } LOG(INFO) << "call_variants done: " << total_examples << " examples, " << total_batches << " batches → " << outfile_path; diff --git a/deepvariant/native/cli.cc b/deepvariant/native/cli.cc index 9c1c9838..10128244 100644 --- a/deepvariant/native/cli.cc +++ b/deepvariant/native/cli.cc @@ -3,9 +3,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -37,7 +39,8 @@ ABSL_FLAG(bool, include_alt_contigs, false, "FILTER parity at WG scale."); ABSL_FLAG(std::string, model_type, "WGS", - "Model type: WGS, WES, PACBIO, ONT, HYBRID_PACBIO_ILLUMINA"); + "Model type: WGS, WES, PACBIO, ONT (alias ONT_R104), " + "HYBRID_PACBIO_ILLUMINA"); ABSL_FLAG(std::string, output_vcf, "", "Output VCF path (run mode)."); ABSL_FLAG(std::string, output_gvcf, "", "Output gVCF path (run mode, optional)."); @@ -53,6 +56,9 @@ ABSL_FLAG(std::string, small_model_path, "", ABSL_DECLARE_FLAG(std::string, reads); ABSL_DECLARE_FLAG(std::string, ref); ABSL_DECLARE_FLAG(std::string, regions); +ABSL_DECLARE_FLAG(std::string, select_variant_types); +ABSL_DECLARE_FLAG(std::string, haploid_contigs); // postprocess_main.cc +ABSL_DECLARE_FLAG(std::string, par_regions_bed); // postprocess_main.cc ABSL_DECLARE_FLAG(int, num_shards); ABSL_DECLARE_FLAG(int, batch_size); ABSL_DECLARE_FLAG(std::string, inference_backend); @@ -80,6 +86,28 @@ inline void AppendAneSpeculateArgs(std::vector& cv_args, "--ane_speculate_confidence=", absl::GetFlag(FLAGS_ane_speculate_confidence))); } + +// Helper: forward the user's --select_variant_types (if any) to make_examples. +// make_examples_main.cc owns the flag and does the actual candidate filtering; +// here we just pass the value through to each per-mode make_examples worker. +inline void AppendSelectVariantTypes(std::vector& me_args) { + const std::string svt = absl::GetFlag(FLAGS_select_variant_types); + if (!svt.empty()) { + me_args.push_back(absl::StrCat("--select_variant_types=", svt)); + } +} + +// Helper: forward sex-chromosome haploid-calling flags to a worker's argv. +// Called for both the make_examples and postprocess args: postprocess_main.cc +// defines the flags and applies the het-zeroing correction, and +// make_examples_main.cc declares them and uses them for haploid-aware gVCF +// reference confidence. The arg vector is the per-mode worker's command line. +inline void AppendHaploidFlags(std::vector& args) { + const std::string hc = absl::GetFlag(FLAGS_haploid_contigs); + if (!hc.empty()) args.push_back(absl::StrCat("--haploid_contigs=", hc)); + const std::string par = absl::GetFlag(FLAGS_par_regions_bed); + if (!par.empty()) args.push_back(absl::StrCat("--par_regions_bed=", par)); +} } // namespace ABSL_DECLARE_FLAG(std::string, checkpoint); // Phase 9 / Step 1 — alt-aligned pileup mode (PacBio/ONT). Defined in @@ -362,6 +390,19 @@ int EffectiveBatchSize() { return AutoBatchSize(); } +// CanonicalModelType — canonicalize user-facing model_type aliases to the +// internal short form this file compares against. Upstream's canonical +// long-read ONT string is `ONT_R104` (also `ONT_R10` / `ONT_R9`), but every +// comparison in this file is against the short `"ONT"`. Without this, a user +// passing --model_type=ONT_R104 would be silently misrouted to the WGS/WES +// default branch. Upper-cases first to match the existing convention in this +// file (all comparisons here are against upper-case forms). +std::string CanonicalModelType(std::string mt) { + for (char& c : mt) c = static_cast(std::toupper(c)); + if (mt == "ONT_R104" || mt == "ONT_R10" || mt == "ONT_R9") return "ONT"; + return mt; +} + std::string ModelPath(const std::string& model_type) { if (!absl::GetFlag(FLAGS_model).empty()) { return absl::GetFlag(FLAGS_model); @@ -375,6 +416,69 @@ std::string ModelPath(const std::string& model_type) { return absl::StrCat(base, "/", type, ".mlpackage"); } +// MergeCvoFiles — concatenate small-model CVO shards (if any) followed by the +// big-model CVO into a single `merged_cvo_path`, by raw byte copy. TFRecord +// allows naive byte concatenation because each record is self-delimiting. +// +// `small_cvo_pattern` is either empty (no small model → only big_cvo is +// copied), a single file path, or a `name@N` shard spec written one-file-per +// make_examples worker thread; in the `@N` case we expand to the conventional +// `name-NNNNN-of-NNNNN` shard files and append each in order. `big_cvo_path` +// is the single file written by call_variants. +// +// This replaces an earlier unquoted `std::system("cat ... > ...")` merge in +// the somatic/pangenome paths, which broke on paths containing spaces or shell +// metacharacters. Returns false (with LOG(ERROR)) on any I/O failure. +// +// Critical subtlety preserved from the germline merge: we read each source +// into a buffer and only write when non-empty. `operator<<(streambuf*)` sets +// the output stream's failbit when the source is empty, which silently breaks +// ALL subsequent writes — and sharded small_cvo routinely has 0-record shards. +bool MergeCvoFiles(const std::string& small_cvo_pattern, + const std::string& big_cvo_path, + const std::string& merged_cvo_path) { + std::ofstream out(merged_cvo_path, std::ios::binary | std::ios::trunc); + if (!out) { + LOG(ERROR) << "Cannot open merged CVO: " << merged_cvo_path; + return false; + } + bool ok = true; + auto append_path = [&](const std::string& p) { + std::ifstream in(p, std::ios::binary); + if (!in) return; // missing shard (e.g. 0-record) → skip, not an error. + std::vector buf((std::istreambuf_iterator(in)), + std::istreambuf_iterator()); + if (!buf.empty()) out.write(buf.data(), buf.size()); + }; + // small_cvo first (matching the germline / `cat small big` order), then big. + if (!small_cvo_pattern.empty()) { + const auto at = small_cvo_pattern.find('@'); + if (at == std::string::npos) { + append_path(small_cvo_pattern); + } else { + const std::string prefix = small_cvo_pattern.substr(0, at); + int nshard = 0; + if (!absl::SimpleAtoi(small_cvo_pattern.substr(at + 1), &nshard) || + nshard <= 0) { + LOG(ERROR) << "Bad small_cvo shard spec: " << small_cvo_pattern; + return false; + } + for (int i = 0; i < nshard; ++i) { + append_path(absl::StrCat(prefix, "-", + absl::Dec(i, absl::kZeroPad5), + "-of-", absl::Dec(nshard, absl::kZeroPad5))); + } + } + } + append_path(big_cvo_path); + out.close(); + if (!out) { + LOG(ERROR) << "Failed writing merged CVO: " << merged_cvo_path; + ok = false; + } + return ok; +} + } // namespace // Forward decls. @@ -909,7 +1013,8 @@ int RunAll(int argc, char** argv) { return RunAllPangenome(argc, argv); } - const std::string model_type = absl::GetFlag(FLAGS_model_type); + const std::string model_type = + CanonicalModelType(absl::GetFlag(FLAGS_model_type)); const std::string reads_flag = absl::GetFlag(FLAGS_reads); const std::string ref_flag = absl::GetFlag(FLAGS_ref); const std::string output_vcf_flag = absl::GetFlag(FLAGS_output_vcf); @@ -918,6 +1023,20 @@ int RunAll(int argc, char** argv) { const std::string tmp_dir = absl::GetFlag(FLAGS_intermediate_results_dir); const int num_shards = EffectiveNumShards(); + // Ensure the intermediate results dir exists before any stage writes into it + // (the trio/somatic/pangenome dispatchers already create it; germline did + // not, so a non-existent --intermediate_results_dir failed mid-pipeline when + // make_examples tried to open its TFRecord output). + { + std::error_code ec; + std::filesystem::create_directories(tmp_dir, ec); + if (ec) { + LOG(ERROR) << "Failed to create intermediate dir " << tmp_dir << ": " + << ec.message(); + return 1; + } + } + if (reads_flag.empty() || ref_flag.empty() || output_vcf_flag.empty()) { LOG(ERROR) << "Usage: deepvariant run --reads= --ref= " "--output_vcf= [--model_type=WGS] [--regions=chr20]"; @@ -1044,6 +1163,8 @@ int RunAll(int argc, char** argv) { if (absl::GetFlag(FLAGS_use_direct_phasing)) { me_args.push_back("--use_direct_phasing=true"); } + AppendSelectVariantTypes(me_args); + AppendHaploidFlags(me_args); // haploid gVCF reference confidence auto argv_me = MakeArgv("deepvariant_make_examples", me_args); int n = static_cast(argv_me.size()) - 1; if (int rc = RunMakeExamples(n, argv_me.data()); rc != 0) { @@ -1084,42 +1205,9 @@ int RunAll(int argc, char** argv) { std::string postprocess_input = cvo_pattern; if (!small_model_path.empty()) { LOG(INFO) << "Stage 2.5: merge small_cvo + big_cvo → " << merged_cvo_path; - std::ofstream out(merged_cvo_path, std::ios::binary | std::ios::trunc); - if (!out) { - LOG(ERROR) << "Cannot open merged CVO: " << merged_cvo_path; + if (!MergeCvoFiles(small_cvo_path, cvo_pattern, merged_cvo_path)) { return 1; } - auto append_path = [&](const std::string& p) { - std::ifstream in(p, std::ios::binary); - if (!in) return; - // Read into buffer first — operator<<(streambuf*) sets failbit when - // the source is empty, which silently breaks ALL subsequent writes. - // Critical for sharded small_cvo where some shards have 0 records. - std::vector buf((std::istreambuf_iterator(in)), - std::istreambuf_iterator()); - if (!buf.empty()) out.write(buf.data(), buf.size()); - }; - // small_cvo: expand "@N" → per-shard files. - auto at = small_cvo_path.find('@'); - if (at == std::string::npos) { - append_path(small_cvo_path); - } else { - const std::string prefix = small_cvo_path.substr(0, at); - int nshard = 0; - if (!absl::SimpleAtoi(small_cvo_path.substr(at + 1), &nshard) || - nshard <= 0) { - LOG(ERROR) << "Bad small_cvo shard spec: " << small_cvo_path; - return 1; - } - for (int i = 0; i < nshard; ++i) { - append_path(absl::StrCat(prefix, "-", - absl::Dec(i, absl::kZeroPad5), - "-of-", absl::Dec(nshard, absl::kZeroPad5))); - } - } - // big_cvo: single file (call_variants writes once). - append_path(cvo_pattern); - out.close(); postprocess_input = merged_cvo_path; } @@ -1138,6 +1226,7 @@ int RunAll(int argc, char** argv) { } // Per-model postprocess flags (e.g. WES multiallelic_mode=min). for (const auto& f : PostprocessModelFlags(model_type)) pp_args.push_back(f); + AppendHaploidFlags(pp_args); auto argv_pp = MakeArgv("deepvariant_postprocess", pp_args); int n = static_cast(argv_pp.size()) - 1; if (int rc = RunPostprocessVariants(n, argv_pp.data()); rc != 0) { @@ -1158,7 +1247,8 @@ int RunAll(int argc, char** argv) { // ────────────────────────────────────────────────────────────────────── int RunAllTrio(int argc, char** argv) { absl::ParseCommandLine(argc, argv); - const std::string model_type = absl::GetFlag(FLAGS_model_type); + const std::string model_type = + CanonicalModelType(absl::GetFlag(FLAGS_model_type)); // Ensure intermediate_results_dir exists (may not be pre-created by caller). { std::system(absl::StrCat("mkdir -p '", absl::GetFlag(FLAGS_intermediate_results_dir), "'").c_str()); } @@ -1393,6 +1483,8 @@ int RunAllTrio(int argc, char** argv) { if (absl::GetFlag(FLAGS_use_direct_phasing)) { me_args.push_back("--use_direct_phasing=true"); } + AppendSelectVariantTypes(me_args); + AppendHaploidFlags(me_args); // haploid gVCF reference confidence auto argv_me = MakeArgv("deepvariant_make_examples", me_args); int n = static_cast(argv_me.size()) - 1; if (int rc = RunMakeExamples(n, argv_me.data()); rc != 0) { @@ -1479,6 +1571,7 @@ int RunAllTrio(int argc, char** argv) { if (!p.output_gvcf.empty()) { pp_args.push_back(absl::StrCat("--gvcf_outfile=", p.output_gvcf)); } + AppendHaploidFlags(pp_args); auto argv_pp = MakeArgv("deepvariant_postprocess", pp_args); int n = static_cast(argv_pp.size()) - 1; if (int rc = RunPostprocessVariants(n, argv_pp.data()); rc != 0) { @@ -1501,7 +1594,8 @@ int RunAllTrio(int argc, char** argv) { // ────────────────────────────────────────────────────────────────────── int RunAllSomatic(int argc, char** argv) { absl::ParseCommandLine(argc, argv); - const std::string model_type = absl::GetFlag(FLAGS_model_type); + const std::string model_type = + CanonicalModelType(absl::GetFlag(FLAGS_model_type)); { std::system(absl::StrCat("mkdir -p '", absl::GetFlag(FLAGS_intermediate_results_dir), "'").c_str()); } const std::string ref_flag = absl::GetFlag(FLAGS_ref); @@ -1633,6 +1727,8 @@ int RunAllSomatic(int argc, char** argv) { me_args.push_back(absl::StrCat("--population_vcfs=", pon)); } } + AppendSelectVariantTypes(me_args); + AppendHaploidFlags(me_args); // haploid gVCF reference confidence auto argv_me = MakeArgv("deepvariant_make_examples", me_args); int n = static_cast(argv_me.size()) - 1; if (int rc = RunMakeExamples(n, argv_me.data()); rc != 0) { @@ -1669,19 +1765,11 @@ int RunAllSomatic(int argc, char** argv) { // ── Stage 2.5: merge small_cvo into cvo (if SM was used). ── LOG(INFO) << "Somatic Stage 2.5: merge → " << merged_cvo_path; { - std::vector cmd = { - "/bin/sh", "-c", - absl::StrCat("cat ", cvo_path, " > ", merged_cvo_path)}; - if (!sm_path.empty()) { - // Pre-pend small_cvo records. - cmd[2] = absl::StrCat( - "cat ", - n_threads > 1 ? absl::StrCat(tmp_dir, "/small_cvo_tumor.tfrecord-*") - : small_cvo_pattern, - " ", cvo_path, " > ", merged_cvo_path); - } - int rc = std::system(cmd[2].c_str()); - if (rc != 0) { + // small_cvo_pattern is a `name@N` shard spec (n_threads>1) or a single + // file; MergeCvoFiles expands `@N` to the per-shard files. When sm_path is + // empty no small_cvo was written, so we merge only the big cvo. + const std::string small = sm_path.empty() ? "" : small_cvo_pattern; + if (!MergeCvoFiles(small, cvo_path, merged_cvo_path)) { LOG(ERROR) << "Somatic: merge step failed"; return 1; } @@ -1706,6 +1794,7 @@ int RunAllSomatic(int argc, char** argv) { pp_args.push_back(absl::StrCat("--pon_filtering=", user_pon)); } } + AppendHaploidFlags(pp_args); auto argv_pp = MakeArgv("deepvariant_postprocess", pp_args); int n = static_cast(argv_pp.size()) - 1; if (int rc = RunPostprocessVariants(n, argv_pp.data()); rc != 0) { @@ -1848,6 +1937,8 @@ int RunAllPangenome(int argc, char** argv) { // 25kb-partition reservoir sampling reduced to ~1, killing the candidate). // partition_size=1000 mirrors Docker's per-1kb reservoir granularity. me_args.push_back("--partition_size=1000"); + AppendSelectVariantTypes(me_args); + AppendHaploidFlags(me_args); // haploid gVCF reference confidence auto argv_me = MakeArgv("deepvariant_make_examples", me_args); int n = static_cast(argv_me.size()) - 1; if (int rc = RunMakeExamples(n, argv_me.data()); rc != 0) { @@ -1882,18 +1973,11 @@ int RunAllPangenome(int argc, char** argv) { // ── Stage 2.5: merge small_cvo into cvo. ── LOG(INFO) << "Pangenome Stage 2.5: merge → " << merged_cvo_path; { - std::vector cmd = { - "/bin/sh", "-c", - absl::StrCat("cat ", cvo_path, " > ", merged_cvo_path)}; - if (!sm_path.empty()) { - cmd[2] = absl::StrCat( - "cat ", - n_threads > 1 ? absl::StrCat(tmp_dir, "/small_cvo_reads.tfrecord-*") - : small_cvo_pattern, - " ", cvo_path, " > ", merged_cvo_path); - } - int rc = std::system(cmd[2].c_str()); - if (rc != 0) { + // small_cvo_pattern is a `name@N` shard spec (n_threads>1) or a single + // file; MergeCvoFiles expands `@N` to the per-shard files. When sm_path is + // empty no small_cvo was written, so we merge only the big cvo. + const std::string small = sm_path.empty() ? "" : small_cvo_pattern; + if (!MergeCvoFiles(small, cvo_path, merged_cvo_path)) { LOG(ERROR) << "Pangenome: merge step failed"; return 1; } @@ -1907,6 +1991,7 @@ int RunAllPangenome(int argc, char** argv) { absl::StrCat("--infile=", merged_cvo_path), absl::StrCat("--output_vcf_outfile=", out_vcf), }; + AppendHaploidFlags(pp_args); auto argv_pp = MakeArgv("deepvariant_postprocess", pp_args); int n = static_cast(argv_pp.size()) - 1; if (int rc = RunPostprocessVariants(n, argv_pp.data()); rc != 0) { diff --git a/deepvariant/native/coreml_inference.mm b/deepvariant/native/coreml_inference.mm index a706a099..9b3a0d9b 100644 --- a/deepvariant/native/coreml_inference.mm +++ b/deepvariant/native/coreml_inference.mm @@ -11,9 +11,59 @@ #include #include #include +#include namespace deepvariant { +namespace { + +// True when `arr`'s memory is densely packed C-contiguous (row-major) for its +// shape, i.e. the stride of each dimension equals the product of the extents of +// all dimensions inside it. Only then is a flat memcpy against a C-contiguous +// host buffer valid. MLMultiArray makes no layout guarantee — `initWithShape:` +// allocations and especially prediction-result arrays can carry padded strides. +bool IsPackedRowMajor(MLMultiArray* arr) { + NSArray* shape = arr.shape; + NSArray* strides = arr.strides; + NSInteger expected = 1; + for (NSInteger d = (NSInteger)shape.count - 1; d >= 0; --d) { + if (strides[d].integerValue != expected) return false; + expected *= shape[d].integerValue; + } + return true; +} + +// Copy `count` floats between a C-contiguous row-major host buffer and `arr`, +// honoring `arr.strides` (element units). When `to_array` is true the host +// buffer is the source (scatter into `arr`); otherwise `arr` is the source +// (gather into the host buffer). `count` must equal the product of the shape. +void StridedCopyFloat(MLMultiArray* arr, float* host, size_t count, + bool to_array) { + float* data = (float*)arr.dataPointer; + const NSUInteger nd = arr.shape.count; + std::vector dims(nd), strides(nd), idx(nd, 0); + for (NSUInteger d = 0; d < nd; ++d) { + dims[d] = arr.shape[d].integerValue; + strides[d] = arr.strides[d].integerValue; + } + for (size_t flat = 0; flat < count; ++flat) { + NSInteger off = 0; + for (NSUInteger d = 0; d < nd; ++d) off += idx[d] * strides[d]; + if (to_array) { + data[off] = host[flat]; + } else { + host[flat] = data[off]; + } + // Increment the multi-index, last dimension varying fastest (row-major). + for (NSInteger d = (NSInteger)nd - 1; d >= 0; --d) { + if (++idx[d] < dims[d]) break; + idx[d] = 0; + } + } +} + +} // namespace + struct CoreMLModel::Impl { MLModel* model = nil; NSString* input_name = @"x"; @@ -119,8 +169,14 @@ NSLog(@"MLMultiArray alloc failed: %@", error.localizedDescription); return false; } - std::memcpy(arr.dataPointer, images, - (size_t)N * (size_t)elemPerImage * sizeof(float)); + const size_t in_count = (size_t)N * (size_t)elemPerImage; + if (IsPackedRowMajor(arr)) { + std::memcpy(arr.dataPointer, images, in_count * sizeof(float)); + } else { + // Padded strides: scatter element-by-element so we write the right cells. + StridedCopyFloat(arr, const_cast(images), in_count, + /*to_array=*/true); + } MLDictionaryFeatureProvider* fp = [[MLDictionaryFeatureProvider alloc] @@ -145,8 +201,22 @@ return false; } // Output is FP32 (we requested it at conversion time) and shape (N, K). - const float* src = (const float*)out_arr.dataPointer; - std::memcpy(probs, src, (size_t)N * (size_t)num_classes * sizeof(float)); + const size_t out_count = (size_t)N * (size_t)num_classes; + // Require an exact logical-element match: a larger array means a different + // shape than (N, num_classes), and copying/gathering out_count elements + // against its real dims would silently scramble the per-genotype + // probabilities (the strided gather walks out_arr's own dims). + if ((size_t)out_arr.count != out_count) { + NSLog(@"Output array has %ld elements, expected %zu", + (long)out_arr.count, out_count); + return false; + } + if (IsPackedRowMajor(out_arr)) { + std::memcpy(probs, out_arr.dataPointer, out_count * sizeof(float)); + } else { + // Core ML may hand back a strided/padded array; gather honoring strides. + StridedCopyFloat(out_arr, probs, out_count, /*to_array=*/false); + } return true; } } diff --git a/deepvariant/native/dump_allele_counts_main.cc b/deepvariant/native/dump_allele_counts_main.cc index c14461a8..1d5b7bfb 100644 --- a/deepvariant/native/dump_allele_counts_main.cc +++ b/deepvariant/native/dump_allele_counts_main.cc @@ -50,8 +50,11 @@ int main(int argc, char** argv) { std::vector se = absl::StrSplit(parts[1], '-'); if (se.size() != 2) { std::fprintf(stderr, "bad region\n"); return 2; } int64_t s = 0, e = 0; - absl::SimpleAtoi(se[0], &s); - absl::SimpleAtoi(se[1], &e); + if (!absl::SimpleAtoi(se[0], &s) || !absl::SimpleAtoi(se[1], &e) || + s < 1 || e < s) { + std::fprintf(stderr, "bad region\n"); // 1-based, non-empty range + return 2; + } region.set_start(s - 1); // 1-based input → 0-based proto. region.set_end(e); } diff --git a/deepvariant/native/extract_pileup_at_pos_main.cc b/deepvariant/native/extract_pileup_at_pos_main.cc index 9d741728..4dc20755 100644 --- a/deepvariant/native/extract_pileup_at_pos_main.cc +++ b/deepvariant/native/extract_pileup_at_pos_main.cc @@ -1,7 +1,8 @@ // Phase 5.5c PASS-flip diagnostic: extract the pileup image at a // specific (chrom, pos, ref, alt) from an examples TFRecord and write -// it as a single (1, 100, 221, 7) NHWC FP32 .npy. Pixel encoding -// matches call_variants ((src - 128) / 128). +// it as a single (1, H, W, C) NHWC FP32 .npy, where H/W/C are read from +// the example's image/shape feature (WGS 100x221x7 when absent). Pixel +// encoding matches call_variants ((src - 128) / 128). // // Used to byte-compare our pileup image against Docker's at the same // site, isolating "inference drift" from "different pileup-image @@ -14,8 +15,10 @@ // start_1based is the conventional VCF coordinate (1-based); // internally we compare against variant.start() which is 0-based. +#include #include #include +#include #include #include #include @@ -65,10 +68,64 @@ std::string ExtractBytesListFirst(const uint8_t* buf, size_t len) { return {}; } -// Parse top-level tf.train.Example, return (image_encoded, variant_encoded). +// Decode a tf.train.Feature message holding an Int64List into its values. +// Handles both packed (proto3 default) and unpacked repeated-int64 encodings. +// Returns empty if the Feature is not an Int64List. +std::vector ParseInt64List(const uint8_t* buf, size_t len) { + std::vector out; + size_t i = 0; + while (i < len) { + uint64_t tag = ReadVarint(buf, len, i); + uint32_t field = static_cast(tag >> 3); + uint32_t wire = static_cast(tag & 7); + if (wire == 2) { + uint64_t seg_len = ReadVarint(buf, len, i); + if (i + seg_len > len) break; + if (field == 3) { // Feature.int64_list + const uint8_t* sub = buf + i; + size_t si = 0; + while (si < seg_len) { + uint64_t vtag = ReadVarint(sub, seg_len, si); + uint32_t vfield = static_cast(vtag >> 3); + uint32_t vwire = static_cast(vtag & 7); + if (vfield == 1 && vwire == 2) { // packed values + uint64_t plen = ReadVarint(sub, seg_len, si); + if (si + plen > seg_len) break; + size_t pend = si + plen; + while (si < pend) { + // Bound reads by pend (the packed-blob end), not seg_len, so a + // truncated trailing varint can't run into the rest of the message. + out.push_back(static_cast(ReadVarint(sub, pend, si))); + } + } else if (vfield == 1 && vwire == 0) { // single unpacked value + out.push_back(static_cast(ReadVarint(sub, seg_len, si))); + } else { + break; + } + } + return out; + } + i += seg_len; + } else if (wire == 0) { + ReadVarint(buf, len, i); + } else if (wire == 5) { + i += 4; + } else if (wire == 1) { + i += 8; + } else { + break; + } + } + return out; +} + +// Parse top-level tf.train.Example, return image bytes, variant bytes, shape. struct ExampleParts { std::string image_encoded; std::string variant_encoded; + std::vector image_shape; // [H, W, C] when present & well-formed. + bool image_shape_present = false; // true if the feature key was seen at all, + // independent of whether it decoded to 3. }; ExampleParts ParseExample(const std::string& payload) { @@ -114,6 +171,11 @@ ExampleParts ParseExample(const std::string& payload) { out.variant_encoded = ExtractBytesListFirst( reinterpret_cast(value_bytes.data()), value_bytes.size()); + } else if (key == "image/shape") { + out.image_shape_present = true; + out.image_shape = ParseInt64List( + reinterpret_cast(value_bytes.data()), + value_bytes.size()); } } } @@ -199,20 +261,28 @@ int main(int argc, char** argv) { const std::string tfr_path = argv[1]; const std::string out_path = argv[2]; const std::string chrom = argv[3]; - const int64_t start_0b = std::strtoll(argv[4], nullptr, 10) - 1; + // Checked parse of the 1-based start position. An unchecked strtoll() would + // silently turn bad input ("abc", "", "12x") into position 0, scanning the + // whole TFRecord and matching nothing — a confusing no-op. Reject anything + // that isn't a complete, in-range positive integer. + errno = 0; + char* end = nullptr; + const long long start_1based = std::strtoll(argv[4], &end, 10); + if (errno != 0 || end == argv[4] || *end != '\0' || start_1based < 1) { + std::fprintf(stderr, "error: must be a positive integer, " + "got \"%s\"\n", argv[4]); + return 2; + } + const int64_t start_0b = static_cast(start_1based) - 1; const std::string ref = argv[5]; const std::string alt = argv[6]; - constexpr int H = 100, W = 221, C = 7; - constexpr int64_t kElem = (int64_t)H * W * C; - auto reader = deepvariant::TFRecordReader::New(tfr_path); if (!reader) { std::fprintf(stderr, "cannot open %s\n", tfr_path.c_str()); return 1; } - std::vector img(kElem); int found = 0; long scanned = 0; while (reader->GetNext()) { @@ -220,6 +290,34 @@ int main(int argc, char** argv) { auto p = ParseExample(reader->record()); if (!VariantMatches(p.variant_encoded, chrom, start_0b, ref, alt)) continue; if (p.image_encoded.empty()) continue; + // Geometry comes from the matched example's image/shape so non-WGS models + // (WES/PacBio/ONT) work; fall back to WGS only when the feature is absent. + // A present-but-non-3-D or out-of-range shape is an error, not a WGS guess. + int H = 100, W = 221, C = 7; + if (p.image_shape_present) { + if (p.image_shape.size() != 3) { + std::fprintf(stderr, + "record %ld: image/shape has %zu values, expected 3 (H,W,C)\n", + scanned, p.image_shape.size()); + return 1; + } + constexpr int64_t kMaxDim = 100000; + const int64_t h = p.image_shape[0], w = p.image_shape[1], + c = p.image_shape[2]; + if (h <= 0 || w <= 0 || c <= 0 || + h > kMaxDim || w > kMaxDim || c > kMaxDim) { + std::fprintf(stderr, + "record %ld: image/shape %lldx%lldx%lld out of range (1..%lld)\n", + scanned, (long long)h, (long long)w, (long long)c, + (long long)kMaxDim); + return 1; + } + H = static_cast(h); + W = static_cast(w); + C = static_cast(c); + } + const int64_t kElem = static_cast(H) * W * C; + std::vector img(kElem); if ((int64_t)p.image_encoded.size() == kElem) { const uint8_t* src = reinterpret_cast(p.image_encoded.data()); constexpr float inv = 1.0f / 128.0f; @@ -230,8 +328,9 @@ int main(int argc, char** argv) { std::memcpy(img.data(), p.image_encoded.data(), (size_t)kElem * sizeof(float)); } else { - std::fprintf(stderr, "record %ld: bad image size %zu\n", scanned, - p.image_encoded.size()); + std::fprintf(stderr, "record %ld: bad image size %zu (expected %lld or " + "%lld for %dx%dx%d)\n", scanned, p.image_encoded.size(), + (long long)kElem, (long long)kElem * 4, H, W, C); continue; } if (!WriteNpyFp32(out_path, img.data(), 1, H, W, C)) { diff --git a/deepvariant/native/extract_pileup_npy_main.cc b/deepvariant/native/extract_pileup_npy_main.cc index 3cf3d771..8c7ab82b 100644 --- a/deepvariant/native/extract_pileup_npy_main.cc +++ b/deepvariant/native/extract_pileup_npy_main.cc @@ -1,6 +1,8 @@ // Profiling tool: extract the first N pileup images from a TFRecord (or // `name@N` shard spec) and write them as a NumPy `.npy` array of shape -// (N, 100, 221, 7) FP32 NHWC. Pixel encoding mirrors call_variants: +// (N, H, W, C) FP32 NHWC, where H/W/C are read from the example's +// image/shape feature (WGS 100x221x7 when absent). Pixel encoding mirrors +// call_variants: // uint8 src → (src - 128) / 128.0 → FP32 // or a passthrough when the input is already FP32. // @@ -68,7 +70,66 @@ std::string ExtractBytesListFirst(const uint8_t* buf, size_t len) { return {}; } -std::string ParseImageEncoded(const std::string& payload) { +// Decode a tf.train.Feature message holding an Int64List into its values. +// Handles both packed (proto3 default) and unpacked repeated-int64 encodings. +// Returns empty if the Feature is not an Int64List. +std::vector ParseInt64List(const uint8_t* buf, size_t len) { + std::vector out; + size_t i = 0; + while (i < len) { + uint64_t tag = ReadVarint(buf, len, i); + uint32_t field = static_cast(tag >> 3); + uint32_t wire = static_cast(tag & 7); + if (wire == 2) { + uint64_t seg_len = ReadVarint(buf, len, i); + if (i + seg_len > len) break; + if (field == 3) { // Feature.int64_list + const uint8_t* sub = buf + i; + size_t si = 0; + while (si < seg_len) { + uint64_t vtag = ReadVarint(sub, seg_len, si); + uint32_t vfield = static_cast(vtag >> 3); + uint32_t vwire = static_cast(vtag & 7); + if (vfield == 1 && vwire == 2) { // packed values + uint64_t plen = ReadVarint(sub, seg_len, si); + if (si + plen > seg_len) break; + size_t pend = si + plen; + while (si < pend) { + // Bound reads by pend (the packed-blob end), not seg_len, so a + // truncated trailing varint can't run into the rest of the message. + out.push_back(static_cast(ReadVarint(sub, pend, si))); + } + } else if (vfield == 1 && vwire == 0) { // single unpacked value + out.push_back(static_cast(ReadVarint(sub, seg_len, si))); + } else { + break; + } + } + return out; + } + i += seg_len; + } else if (wire == 0) { + ReadVarint(buf, len, i); + } else if (wire == 5) { + i += 4; + } else if (wire == 1) { + i += 8; + } else { + break; + } + } + return out; +} + +struct ParsedExample { + std::string image_encoded; + std::vector image_shape; // [H, W, C] when present & well-formed. + bool image_shape_present = false; // true if the feature key was seen at all, + // independent of whether it decoded to 3. +}; + +ParsedExample ParseExample(const std::string& payload) { + ParsedExample out; const uint8_t* buf = reinterpret_cast(payload.data()); size_t n = payload.size(); size_t i = 0; @@ -106,16 +167,21 @@ std::string ParseImageEncoded(const std::string& payload) { ei += elen; } if (key == "image/encoded" || key == "image") { - return ExtractBytesListFirst( + out.image_encoded = ExtractBytesListFirst( + reinterpret_cast(value_bytes.data()), + value_bytes.size()); + } else if (key == "image/shape") { + out.image_shape_present = true; + out.image_shape = ParseInt64List( reinterpret_cast(value_bytes.data()), value_bytes.size()); } } } - return {}; + return out; } -// Write a (N, 100, 221, 7) FP32 NHWC array to NumPy v1 .npy. +// Write a (N, H, W, C) FP32 NHWC array to NumPy v1 .npy. bool WriteNpyFp32NHWC(const std::string& path, int N, int H, int W, int C, const float* data) { std::ofstream f(path, std::ios::binary); @@ -161,24 +227,76 @@ int main(int argc, char** argv) { return 2; } - // Standard WGS DeepVariant pileup geometry. - constexpr int H = 100, W = 221, C = 7; - constexpr int64_t kElemPerImg = static_cast(H) * W * C; - auto reader = deepvariant::TFRecordReader::New(tfr_path); if (!reader) { std::fprintf(stderr, "cannot open %s\n", tfr_path.c_str()); return 1; } - std::vector all(static_cast(count) * kElemPerImg); + // Geometry is taken from the first record's image/shape feature so the tool + // adapts to any model type (WES/PacBio/ONT differ from WGS). The whole batch + // is packed into one (N, H, W, C) array, so later records must share the + // geometry — a mismatch is an error rather than a silent overwrite. + int H = 0, W = 0, C = 0; + int64_t kElemPerImg = 0; + std::vector all; int n_loaded = 0; for (int i = 0; i < count; ++i) { if (!reader->GetNext()) { std::fprintf(stderr, "EOF after %d records\n", i); break; } - const std::string img = ParseImageEncoded(reader->record()); + const ParsedExample ex = ParseExample(reader->record()); + if (i == 0) { + if (!ex.image_shape_present) { + H = 100; W = 221; C = 7; // WGS fallback only when the feature is absent. + std::fprintf(stderr, + "warning: record 0 has no image/shape; assuming WGS %dx%dx%d\n", + H, W, C); + } else if (ex.image_shape.size() != 3) { + // Present but not a 3-D [H,W,C] shape — don't silently guess WGS. + std::fprintf(stderr, + "record 0: image/shape has %zu values, expected 3 (H,W,C)\n", + ex.image_shape.size()); + return 1; + } else { + // Validate each int64 dim is a sane positive value before narrowing to + // int, so a corrupt shape can't wrap to garbage or request an absurd + // allocation. + constexpr int64_t kMaxDim = 100000; + const int64_t h = ex.image_shape[0], w = ex.image_shape[1], + c = ex.image_shape[2]; + if (h <= 0 || w <= 0 || c <= 0 || + h > kMaxDim || w > kMaxDim || c > kMaxDim) { + std::fprintf(stderr, + "record 0: image/shape %lldx%lldx%lld out of range (1..%lld)\n", + (long long)h, (long long)w, (long long)c, (long long)kMaxDim); + return 1; + } + H = static_cast(h); + W = static_cast(w); + C = static_cast(c); + } + kElemPerImg = static_cast(H) * W * C; + // count is already bounded to (0, 100000]; guard the batch allocation. + if (kElemPerImg <= 0 || kElemPerImg > INT64_MAX / count) { + std::fprintf(stderr, "record 0: batch too large (%lld elems x %d)\n", + (long long)kElemPerImg, count); + return 1; + } + all.resize(static_cast(count) * kElemPerImg); + } else if (ex.image_shape.size() == 3 && + (ex.image_shape[0] != H || ex.image_shape[1] != W || + ex.image_shape[2] != C)) { + std::fprintf(stderr, + "record %d: image/shape %lldx%lldx%lld differs from batch %dx%dx%d; " + "cannot pack heterogeneous geometry into one array\n", + i, static_cast(ex.image_shape[0]), + static_cast(ex.image_shape[1]), + static_cast(ex.image_shape[2]), H, W, C); + return 1; + } + const std::string& img = ex.image_encoded; float* dst = all.data() + static_cast(i) * kElemPerImg; if (static_cast(img.size()) == kElemPerImg) { // uint8 → (x - 128) / 128 — same path as call_variants. @@ -201,6 +319,13 @@ int main(int argc, char** argv) { ++n_loaded; } + if (n_loaded == 0) { + // Empty input / immediate EOF: don't emit a bogus (0,0,0,0) .npy. + std::fprintf(stderr, "no records read from %s; nothing written\n", + tfr_path.c_str()); + return 1; + } + if (!WriteNpyFp32NHWC(out_path, n_loaded, H, W, C, all.data())) { std::fprintf(stderr, "failed to write %s\n", out_path.c_str()); return 1; diff --git a/deepvariant/native/gvcf_emit.cc b/deepvariant/native/gvcf_emit.cc index ed178f06..0cc1c7e5 100644 --- a/deepvariant/native/gvcf_emit.cc +++ b/deepvariant/native/gvcf_emit.cc @@ -43,22 +43,31 @@ int Log10PtrueToPhred(double log10_p_ref, int max_gq) { return std::min(gq, max_gq); } +// log10 of an effectively-impossible probability, used to drive the +// heterozygous likelihood to zero on haploid contigs (matches upstream +// variant_caller.py:IMPOSSIBLE_PROBABILITY_LOG10 = 999.0). +constexpr double kImpossibleProbabilityLog10 = 999.0; + // Compute reference-confidence likelihoods + GQ for one site, given ref/total // read counts and per-base error rate. Returns log10 probs in [ref, het, alt]. +// When `is_haploid` (a haploid contig outside the PAR), the heterozygous +// likelihood is forced to zero, mirroring variant_caller.py's is_haploid +// branch in _calc_reference_confidence. void ReferenceConfidence(int n_ref, int n_total, double p_error, - double* log10_p_ref, double* log10_p_het, - double* log10_p_alt) { + bool is_haploid, double* log10_p_ref, + double* log10_p_het, double* log10_p_alt) { if (n_total <= 0) { - // No coverage: uniform. + // No coverage: uniform over the possible genotypes. Haploid drops het. *log10_p_ref = -1.0; - *log10_p_het = -1.0; + *log10_p_het = is_haploid ? -kImpossibleProbabilityLog10 : -1.0; *log10_p_alt = -1.0; } else { const int n_alts = n_total - n_ref; const double logp = std::log(p_error) / kLog10; const double log1p = std::log1p(-p_error) / kLog10; *log10_p_ref = n_ref * log1p + n_alts * logp; - *log10_p_het = -n_total * std::log10(2.0); + *log10_p_het = is_haploid ? -kImpossibleProbabilityLog10 + : -n_total * std::log10(2.0); *log10_p_alt = n_ref * logp + n_alts * log1p; } NormalizeLog10Probs(log10_p_ref, log10_p_het, log10_p_alt); @@ -98,10 +107,15 @@ std::vector MakeGvcfRows( const std::vector& summaries, const std::string& sample_name, - double p_error, int gq_resolution, int max_gq, bool include_med_dp) { + double p_error, int gq_resolution, int max_gq, bool include_med_dp, + const std::set* haploid_contigs, + const ParRegions* par_regions) { std::vector out; if (summaries.empty()) return out; + static const ParRegions kNoParRegions; + const ParRegions& par = par_regions ? *par_regions : kNoParRegions; + // 1. Compute per-site GQ + likelihoods. std::vector entries; entries.reserve(summaries.size()); @@ -115,9 +129,13 @@ std::vector MakeGvcfRows( // Skip non-canonical (N, IUPAC) — upstream does the same. continue; } + const bool is_haploid = + haploid_contigs && + IsHaploidPosition(e.ref_name, e.position, e.position + 1, + *haploid_contigs, par); ReferenceConfidence(s.ref_supporting_read_count(), s.total_read_count(), - p_error, &e.log10_probs[0], &e.log10_probs[1], - &e.log10_probs[2]); + p_error, is_haploid, &e.log10_probs[0], + &e.log10_probs[1], &e.log10_probs[2]); e.raw_gq = Log10PtrueToPhred(e.log10_probs[0], max_gq); e.quantized_gq = QuantizeGq(e.raw_gq, gq_resolution); e.gl_is_valid = @@ -129,6 +147,13 @@ std::vector MakeGvcfRows( // 2. Group consecutive entries with same (quantized_gq, gl_is_valid). Emit // one merged Variant row per group when gl_is_valid; emit one Variant per // site when not (uncalled `./.` rows). + // + // is_haploid is deliberately NOT part of the grouping key: upstream + // variant_caller.make_gvcfs groups only on (quantized_gq, has_valid_gl), so + // keying on ploidy here would emit more blocks than upstream and break gVCF + // parity. In practice a haploid site drops its het mass and renormalizes to a + // higher ref probability (higher GQ), so it usually lands in a different + // quantized-GQ bin than an adjacent diploid site and won't merge anyway. size_t i = 0; while (i < entries.size()) { const SiteEntry& first = entries[i]; @@ -157,7 +182,12 @@ std::vector MakeGvcfRows( dps.push_back(entries[k].n_total); } std::sort(dps.begin(), dps.end()); - const int med_dp = dps[dps.size() / 2]; + // True median: average the two middle values on even-length input, matching + // upstream int(statistics.median(...)). + const size_t n = dps.size(); + const int med_dp = (n % 2 == 1) + ? dps[n / 2] + : static_cast((dps[n / 2 - 1] + dps[n / 2]) / 2); if (first.gl_is_valid) { // Emit ONE merged Variant for [i, j). diff --git a/deepvariant/native/gvcf_emit.h b/deepvariant/native/gvcf_emit.h index faafef56..02958faf 100644 --- a/deepvariant/native/gvcf_emit.h +++ b/deepvariant/native/gvcf_emit.h @@ -11,9 +11,11 @@ #pragma once +#include #include #include +#include "deepvariant/native/haploid_regions.h" #include "deepvariant/protos/deepvariant.pb.h" #include "third_party/nucleus/protos/variants.pb.h" @@ -30,10 +32,13 @@ namespace deepvariant { // default). // max_gq: upper cap on GQ (typical: 50). // include_med_dp: emit MED_DP info field (default false). +// haploid_contigs: contigs to call as haploid (null/empty = all diploid). +// par_regions: PAR intervals exempted from haploid calling (null = none). // // Returns: coordinate-sorted Variant protos with `<*>` alt and // `END` info field. Each row may span multiple positions if their -// quantized GQs match. +// quantized GQs match. On a haploid contig outside the PAR, the +// reference-confidence het likelihood is forced to zero. std::vector MakeGvcfRows( const std::vector& summaries, @@ -41,6 +46,8 @@ std::vector MakeGvcfRows( double p_error = 1e-3, int gq_resolution = 5, // Matches upstream --gvcf_gq_binsize default. int max_gq = 50, - bool include_med_dp = false); + bool include_med_dp = false, + const std::set* haploid_contigs = nullptr, + const ParRegions* par_regions = nullptr); } // namespace deepvariant diff --git a/deepvariant/native/haploid_regions.cc b/deepvariant/native/haploid_regions.cc new file mode 100644 index 00000000..18f8e723 --- /dev/null +++ b/deepvariant/native/haploid_regions.cc @@ -0,0 +1,87 @@ +/* + * Copyright 2025 Google LLC. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +#include "deepvariant/native/haploid_regions.h" + +#include + +#include +#include +#include + +#include "absl/strings/match.h" +#include "absl/strings/numbers.h" +#include "absl/strings/str_cat.h" +#include "absl/strings/str_split.h" +#include "absl/strings/string_view.h" + +namespace deepvariant { + +bool LoadParRegions(const std::string& path, ParRegions* out, + std::string* err) { + // gzopen transparently reads both plaintext and gzipped/bgzipped input, so a + // .bed or .bed.gz both work (matching upstream's htslib BedReader). + gzFile f = gzopen(path.c_str(), "rb"); + if (f == nullptr) { + *err = absl::StrCat("cannot open --par_regions_bed: ", path); + return false; + } + // BED lines are short; 64 KiB is far more than any real interval line. + char buf[1 << 16]; + int lineno = 0; + bool ok = true; + while (gzgets(f, buf, sizeof(buf)) != nullptr) { + ++lineno; + absl::string_view line(buf); + while (!line.empty() && (line.back() == '\n' || line.back() == '\r')) { + line.remove_suffix(1); + } + if (line.empty() || line[0] == '#') continue; + if (absl::StartsWith(line, "track") || absl::StartsWith(line, "browser")) { + continue; + } + std::vector cols = + absl::StrSplit(line, absl::ByAnyChar("\t "), absl::SkipEmpty()); + int64_t start, end; + if (cols.size() < 3 || !absl::SimpleAtoi(cols[1], &start) || + !absl::SimpleAtoi(cols[2], &end) || start < 0 || end < start) { + *err = absl::StrCat("malformed BED line ", lineno, " in ", path, ": ", + line); + ok = false; + break; + } + out->by_contig[std::string(cols[0])].emplace_back(start, end); + } + gzclose(f); + return ok; +} + +} // namespace deepvariant diff --git a/deepvariant/native/haploid_regions.h b/deepvariant/native/haploid_regions.h new file mode 100644 index 00000000..aa5dd48f --- /dev/null +++ b/deepvariant/native/haploid_regions.h @@ -0,0 +1,127 @@ +/* + * Copyright 2025 Google LLC. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from this + * software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +// Sex-chromosome haploid-calling support shared by the native make_examples +// (gVCF reference confidence) and postprocess (genotype correction) stages. +// Mirrors upstream's --haploid_contigs / --par_regions_bed handling in +// variant_caller.py and postprocess_variants.py: a position on a haploid +// contig that does not overlap a pseudo-autosomal region is called haploid. + +#ifndef LEARNING_GENOMICS_DEEPVARIANT_NATIVE_HAPLOID_REGIONS_H_ +#define LEARNING_GENOMICS_DEEPVARIANT_NATIVE_HAPLOID_REGIONS_H_ + +#include +#include +#include +#include +#include +#include +#include + +#include "absl/strings/str_split.h" +#include "absl/strings/string_view.h" + +namespace deepvariant { + +// Pseudo-autosomal regions loaded from a BED file: contig -> 0-based +// half-open [start, end) intervals. A position on a haploid contig that +// overlaps one of these stays diploid. +struct ParRegions { + std::map>> by_contig; + + bool Overlaps(const std::string& chrom, int64_t start, int64_t end) const { + auto it = by_contig.find(chrom); + if (it == by_contig.end()) return false; + for (const auto& [s, e] : it->second) { + if (start < e && end > s) return true; + } + return false; + } +}; + +// Parse the whitespace/comma-separated --haploid_contigs flag into a set, +// mirroring upstream's split on ',' then on whitespace. +inline std::set ParseHaploidContigs(const std::string& flag) { + std::set out; + for (absl::string_view tok : + absl::StrSplit(flag, absl::ByAnyChar(", \t\r\n"), absl::SkipEmpty())) { + out.emplace(tok); + } + return out; +} + +// Load a PAR BED (chrom \t start \t end, 0-based half-open). Transparently +// reads plaintext or bgzipped/gzipped BED (matching upstream's htslib +// BedReader). Returns false on an unreadable file or a malformed line so a +// typo'd path fails fast rather than silently disabling the exemption. +// Defined in haploid_regions.cc (uses zlib, so it is kept out of this +// otherwise dependency-free, header-only API). +bool LoadParRegions(const std::string& path, ParRegions* out, std::string* err); + +// True when [start, end) on `chrom` should be called haploid: the contig is in +// `haploid_contigs` and the interval does not overlap a PAR region. Mirror of +// `reference_name in haploid_contigs and not par_regions.overlaps(...)`. +inline bool IsHaploidPosition(const std::string& chrom, int64_t start, + int64_t end, + const std::set& haploid_contigs, + const ParRegions& par_regions) { + if (haploid_contigs.empty()) return false; + if (!haploid_contigs.count(chrom)) return false; + return !par_regions.Overlaps(chrom, start, end); +} + +// Zero every heterozygous genotype's probability and renormalize, forcing a +// haploid call. Probabilities are in PL order F(j/k)=k*(k+1)/2+j (j<=k); het +// genotypes are those with j != k. Mirror of +// postprocess_variants.py:correct_nonautosome_probabilities. +inline void CorrectNonautosomeProbabilities(std::vector* like, + int n_alts) { + // Guard the caller's contract: `like` must hold one entry per diploid + // genotype for n_alts alts, i.e. (n_alts+1)(n_alts+2)/2. A mismatch would + // otherwise be a silent out-of-bounds write below. + if (n_alts < 0 || + like->size() != static_cast((n_alts + 1) * (n_alts + 2) / 2)) { + return; + } + for (int k = 0; k <= n_alts; ++k) { + for (int j = 0; j <= k; ++j) { + if (j != k) (*like)[k * (k + 1) / 2 + j] = 0.0; + } + } + double sum = 0.0; + for (double v : *like) sum += v; + if (sum <= 0.0) sum = 1.0; + for (double& v : *like) v /= sum; +} + +} // namespace deepvariant + +#endif // LEARNING_GENOMICS_DEEPVARIANT_NATIVE_HAPLOID_REGIONS_H_ diff --git a/deepvariant/native/haplotypes.cc b/deepvariant/native/haplotypes.cc index 827548b3..57e03b06 100644 --- a/deepvariant/native/haplotypes.cc +++ b/deepvariant/native/haplotypes.cc @@ -146,16 +146,12 @@ double Log10SumExp(const std::vector& xs) { return m + std::log10(s); } -// Subtract-max + log10(probs / sum) — mirror of +// Subtract log10sumexp so 10^x sums to 1 — mirror of // genomics_math.normalize_log10_probs. std::vector NormalizeLog10Probs(std::vector v) { if (v.empty()) return v; - double m = -std::numeric_limits::infinity(); - for (double x : v) m = std::max(m, x); - for (double& x : v) x -= m; // approximation: subtract max - // Exact: also normalise so 10^x sums to 1. The approximation upstream - // uses is the subtract-max version (line: scaled = [x - m for x in xs]); - // but it's not divided by sum. Both produce the same argmax. + const double lse = Log10SumExp(v); + for (double& x : v) x = std::min(x - lse, 0.0); return v; } diff --git a/deepvariant/native/make_examples_main.cc b/deepvariant/native/make_examples_main.cc index e897431d..a64c8b17 100644 --- a/deepvariant/native/make_examples_main.cc +++ b/deepvariant/native/make_examples_main.cc @@ -28,6 +28,7 @@ #include "deepvariant/direct_phasing.h" #include "deepvariant/make_examples_native.h" #include "deepvariant/native/gvcf_emit.h" +#include "deepvariant/native/haploid_regions.h" #include "deepvariant/native/numpy_mt19937.h" #include "deepvariant/native/realigner_native.h" #include "deepvariant/native/regions.h" @@ -151,6 +152,21 @@ ABSL_FLAG(int, partition_size, 1000, // stricter thresholds (20 / 14). ABSL_FLAG(int, min_mapping_quality, 5, "Min read mapping quality."); ABSL_FLAG(int, min_base_quality, 10, "Min base quality."); +// Mirrors make_examples_options.py's --select_variant_types: a +// whitespace-separated subset of {snps, indels, insertions, deletions, +// multi-allelics, all}. When set, only candidates whose variant matches one +// of the selectors (OR'd) are kept; empty means keep everything. +ABSL_FLAG(std::string, select_variant_types, "", + "Whitespace-separated variant types to keep when generating " + "examples: snps, indels, insertions, deletions, multi-allelics, " + "all. snps/indels/insertions/deletions select bi-allelic variants " + "of that type; multi-allelics selects any multi-allelic variant. " + "Empty (default) keeps all candidates."); +// Sex-chromosome haploid calling. The flags are DEFINED in postprocess_main.cc +// (same multi-call binary, so defining them twice would abort at startup); we +// only read them here to drive haploid gVCF reference confidence. +ABSL_DECLARE_FLAG(std::string, haploid_contigs); +ABSL_DECLARE_FLAG(std::string, par_regions_bed); // Small model first-pass. ABSL_FLAG(std::string, small_model, "", "Path to the small_model .mlpackage. Empty = no small model " @@ -1120,6 +1136,148 @@ bool IsSnpForIndices(const nucleus::genomics::v1::Variant& v, return true; } +// --------------------------------------------------------------------------- +// select_variant_types candidate filtering +// +// Faithful port of make_examples_core.filter_candidates + +// nucleus/util/variant_utils. The variant-type predicates exclude the gVCF +// '<*>' allele, the '' symbolic allele, and the '.' missing field +// from the alt set (variant_utils._non_excluded_alts) before classifying. +// --------------------------------------------------------------------------- + +// True for alts ignored by the type predicates (variant_utils default set). +bool IsExcludedAlt(const std::string& alt) { + return alt == "<*>" || alt == "" || alt == "."; +} + +// Alt alleles that count toward variant-type classification (non-excluded). +std::vector RelevantAlts( + const nucleus::genomics::v1::Variant& v) { + std::vector alts; + for (const auto& a : v.alternate_bases()) { + if (!IsExcludedAlt(a)) alts.push_back(&a); + } + return alts; +} + +// is_snp: REF is 1 bp and every non-excluded alt is 1 bp (>=1 such alt). +bool VarIsSnp(const nucleus::genomics::v1::Variant& v) { + const auto alts = RelevantAlts(v); + if (v.reference_bases().size() != 1 || alts.empty()) return false; + for (const auto* a : alts) { + if (a->size() != 1) return false; + } + return true; +} + +// is_indel: at least one non-excluded alt, and REF>1 or some alt>1 bp. +bool VarIsIndel(const nucleus::genomics::v1::Variant& v) { + const auto alts = RelevantAlts(v); + if (alts.empty()) return false; + if (v.reference_bases().size() > 1) return true; + for (const auto* a : alts) { + if (a->size() > 1) return true; + } + return false; +} + +bool VarIsBiallelic(const nucleus::genomics::v1::Variant& v) { + return RelevantAlts(v).size() == 1; +} + +bool VarIsMultiallelic(const nucleus::genomics::v1::Variant& v) { + return RelevantAlts(v).size() > 1; +} + +// has_insertion/has_deletion gate on is_indel but, like variant_utils, test +// the length condition over ALL alternate_bases (not just the non-excluded). +bool VarHasInsertion(const nucleus::genomics::v1::Variant& v) { + if (!VarIsIndel(v)) return false; + const size_t ref_len = v.reference_bases().size(); + for (const auto& a : v.alternate_bases()) { + if (ref_len < a.size()) return true; + } + return false; +} + +bool VarHasDeletion(const nucleus::genomics::v1::Variant& v) { + if (!VarIsIndel(v)) return false; + const size_t ref_len = v.reference_bases().size(); + for (const auto& a : v.alternate_bases()) { + if (ref_len > a.size()) return true; + } + return false; +} + +// Parsed --select_variant_types flag. `active` is false when no selector was +// requested, in which case filtering is a no-op (keep all candidates). +struct SelectedVariantTypes { + bool active = false; + bool snps = false; + bool indels = false; + bool insertions = false; + bool deletions = false; + bool multiallelics = false; + bool all = false; +}; + +// Parse the whitespace-separated flag. Mirrors make_examples_options.py: an +// unknown selector is a fatal command-line error (returns false + *err set). +bool ParseSelectVariantTypes(const std::string& flag, + SelectedVariantTypes* out, std::string* err) { + for (absl::string_view tok : + absl::StrSplit(flag, absl::ByAnyChar(" \t\r\n"), absl::SkipEmpty())) { + if (tok == "snps") { + out->snps = true; + } else if (tok == "indels") { + out->indels = true; + } else if (tok == "insertions") { + out->insertions = true; + } else if (tok == "deletions") { + out->deletions = true; + } else if (tok == "multi-allelics") { + out->multiallelics = true; + } else if (tok == "all") { + out->all = true; + } else { + *err = absl::StrCat( + "Select variant type '", tok, + "' not recognized. Allowed values are snps, indels, insertions, " + "deletions, multi-allelics, all"); + return false; + } + out->active = true; + } + return true; +} + +// Does the candidate's variant match any requested selector (OR'd)? +// snps/indels/insertions/deletions are bi-allelic-gated, matching the +// VARIANT_TYPE_SELECTORS table in make_examples_core. +bool CandidateSelected(const nucleus::genomics::v1::Variant& v, + const SelectedVariantTypes& sel) { + if (sel.all) return true; + if (sel.snps && VarIsSnp(v) && VarIsBiallelic(v)) return true; + if (sel.indels && VarIsIndel(v) && VarIsBiallelic(v)) return true; + if (sel.insertions && VarHasInsertion(v) && VarIsBiallelic(v)) return true; + if (sel.deletions && VarHasDeletion(v) && VarIsBiallelic(v)) return true; + if (sel.multiallelics && VarIsMultiallelic(v)) return true; + return false; +} + +// Drop candidates whose variant matches no requested selector, preserving +// order. No-op when no selector is active. +void FilterCandidatesBySelectedTypes(std::vector* candidates, + const SelectedVariantTypes& sel) { + if (!sel.active) return; + candidates->erase( + std::remove_if(candidates->begin(), candidates->end(), + [&](const DeepVariantCall& c) { + return !CandidateSelected(c.variant(), sel); + }), + candidates->end()); +} + // Phred = -10 * log10(p), truncated toward zero. Capped at 99. // // Truncation (not std::round) matches upstream's small_model @@ -1241,6 +1399,32 @@ int RunMakeExamples(int argc, char** argv) { LOG(ERROR) << "Required: --ref"; return 1; } + + // Validate --select_variant_types up front so a typo fails fast rather than + // silently keeping everything (the pre-fix behavior) or erroring mid-region. + SelectedVariantTypes selected_types; + { + std::string err; + if (!ParseSelectVariantTypes(absl::GetFlag(FLAGS_select_variant_types), + &selected_types, &err)) { + LOG(ERROR) << err; + return 1; + } + } + + // Sex-chromosome haploid calling config — drives haploid gVCF reference + // confidence on the --haploid_contigs (outside the PAR). + const std::set haploid_contigs = + ParseHaploidContigs(absl::GetFlag(FLAGS_haploid_contigs)); + ParRegions par_regions; + if (const std::string par_bed = absl::GetFlag(FLAGS_par_regions_bed); + !par_bed.empty()) { + std::string err; + if (!LoadParRegions(par_bed, &par_regions, &err)) { + LOG(ERROR) << err; + return 1; + } + } if (reads_path.empty() && !IsSomaticMode()) { LOG(ERROR) << "Required: --reads (or --reads_tumor for somatic mode)"; return 1; @@ -1740,6 +1924,20 @@ int RunMakeExamples(int argc, char** argv) { std::vector candidates = caller.CallsFromAlleleCounts(ac_map, C.name, C.role); + // DirectPhasing (is_phased/PS, below) must see the FULL candidate set so + // its SNP phasing graph matches upstream, which phases in + // candidates_in_region *before* filter_candidates. Capture the + // pre-prune set here (only when phasing is on) and phase that; the + // select_variant_types prune + small_model dispatch then run on the + // filtered set, and is_phased/PS is applied to the surviving + // big_candidates by position. The multi-sample path still filters + // before small_model so pruned types are not emitted via the + // small_model CVO. + std::vector phasing_candidates; + if (absl::GetFlag(FLAGS_use_direct_phasing)) { + phasing_candidates = candidates; // full set, pre-prune + } + FilterCandidatesBySelectedTypes(&candidates, selected_types); if (candidates.empty()) continue; C.total_candidates += candidates.size(); @@ -1844,7 +2042,10 @@ int RunMakeExamples(int argc, char** argv) { for (auto& r : reads_per_sample_v[s]) dp_read_ptrs.emplace_back(&r); ::learning::genomics::deepvariant::DirectPhasing dp( opts.direct_phasing_options()); - auto so = dp.PhaseReads(absl::MakeSpan(big_candidates), + // Phase the full pre-prune candidate set (not just big_candidates) + // so the phasing graph matches upstream; is_phased/PS is then applied + // to the surviving big_candidates by position below. + auto so = dp.PhaseReads(absl::MakeSpan(phasing_candidates), absl::MakeSpan(dp_read_ptrs)); if (so.ok()) { const auto phased = dp.GetPhasedVariants(); @@ -2075,7 +2276,11 @@ int RunMakeExamples(int argc, char** argv) { } LOG(INFO) << " read " << reads.size() << " reads from BAM"; - if (reads.empty()) continue; + // Upstream make_examples_core.py only aborts a zero-coverage region early + // when gVCF is disabled; with a gvcf_writer we must still fall through to + // emit per-position reference-confidence rows (the gVCF block below tolerates + // empty reads, and probe_candidates.empty() then `continue`s with no work). + if (reads.empty() && gvcf_writer == nullptr) continue; // RNA-seq: split reads on N (SKIP) CIGAR ops into per-exon sub-reads // before candidate discovery / realignment / pileup. Mirrors upstream @@ -2190,7 +2395,8 @@ int RunMakeExamples(int argc, char** argv) { absl::GetFlag(FLAGS_p_error), absl::GetFlag(FLAGS_gvcf_gq_binsize), /*max_gq=*/50, - absl::GetFlag(FLAGS_include_med_dp)); + absl::GetFlag(FLAGS_include_med_dp), + &haploid_contigs, &par_regions); for (const auto& v : gvcf_rows) { std::string serialized; v.SerializeToString(&serialized); @@ -2226,8 +2432,6 @@ int RunMakeExamples(int argc, char** argv) { caller.CallsFromAlleleCounter(counter); if (candidates.empty()) continue; - total_candidates += candidates.size(); - // Phase 5.5d/14 — DirectPhasing runs BEFORE small_model dispatch so the // 106-feature haplotype-expanded small_model (PacBio/ONT) sees the same // per-read phase that upstream's FeatureEncoder does. Upstream order @@ -2344,6 +2548,16 @@ int RunMakeExamples(int argc, char** argv) { } } + // select_variant_types pruning runs AFTER DirectPhasing so the phasing + // graph — and the per-read HP tags it produces, which feed the haplotype + // small_model's 106-feature vector — is built from the full candidate set. + // This matches upstream's order: filter_candidates runs in process() only + // after candidates_in_region has already phased. Pruning earlier would drop + // the SNP backbone DirectPhasing relies on and silently change phasing. + FilterCandidatesBySelectedTypes(&candidates, selected_types); + if (candidates.empty()) continue; + total_candidates += candidates.size(); + // Small-model first-pass dispatch. Mirror of upstream // `SmallModelVariantCaller.call_variants` + `make_small_model_examples. // get_set_of_allele_indices`: diff --git a/deepvariant/native/metal_inference.mm b/deepvariant/native/metal_inference.mm index 05a68a28..b6e99ade 100644 --- a/deepvariant/native/metal_inference.mm +++ b/deepvariant/native/metal_inference.mm @@ -1255,6 +1255,12 @@ FusedConv FoldConvBn(const DvwWeights& dvw, int conv_n, int bn_n) { } [cb commit]; [cb waitUntilCompleted]; + if (cb.status == MTLCommandBufferStatusError) { + LOG(ERROR) << "MetalInception::Predict(det): GPU command buffer failed: " + << (cb.error ? cb.error.localizedDescription.UTF8String + : "unknown"); + return false; + } // Phase 8 / Tier 6.0 — full-network det path: bypass MPSGraph entirely // by chaining all 11 Inception blocks + global avg pool, then copying @@ -1288,6 +1294,13 @@ FusedConv FoldConvBn(const DvwWeights& dvw, int conv_n, int bn_n) { } [cb_blk commit]; [cb_blk waitUntilCompleted]; + if (cb_blk.status == MTLCommandBufferStatusError) { + LOG(ERROR) << "MetalInception::Predict(det): GPU command buffer failed: " + << (cb_blk.error + ? cb_blk.error.localizedDescription.UTF8String + : "unknown"); + return false; + } // Copy (B, feature_dim) FP32 result to `output`. const size_t out_bytes = (size_t)batch_size * I.feature_dim * sizeof(float); std::memcpy(output, [I.gap_out_buf contents], out_bytes); diff --git a/deepvariant/native/microtest_neon_cigar_classify.cc b/deepvariant/native/microtest_neon_cigar_classify.cc index bcab89a5..d42b4d19 100644 --- a/deepvariant/native/microtest_neon_cigar_classify.cc +++ b/deepvariant/native/microtest_neon_cigar_classify.cc @@ -125,7 +125,7 @@ int main() { static const char alphabet[] = "ACGTNacgtX0"; std::mt19937 rng(0xFEEDFACEu); for (size_t i = 0; i < n; ++i) - read_buf[i] = alphabet[rng() % sizeof(alphabet) - 1]; + read_buf[i] = alphabet[rng() % (sizeof(alphabet) - 1)]; for (uint8_t q : test_quals) { qual_buf.assign(n, q); for (int leg = 0; leg <= 1; ++leg) { diff --git a/deepvariant/native/postprocess_main.cc b/deepvariant/native/postprocess_main.cc index 740e17b0..da3ca3ae 100644 --- a/deepvariant/native/postprocess_main.cc +++ b/deepvariant/native/postprocess_main.cc @@ -21,6 +21,8 @@ #include #include #include +#include +#include #include #include #include @@ -28,6 +30,7 @@ #include #include +#include "deepvariant/native/haploid_regions.h" #include "deepvariant/native/haplotypes.h" #include "deepvariant/native/tfrecord.h" #include "deepvariant/protos/deepvariant.pb.h" @@ -98,6 +101,17 @@ ABSL_FLAG(bool, process_somatic, false, // "min" (WES): take minimum probability across kept CVOs. ABSL_FLAG(std::string, multiallelic_mode, "product", "Multi-allelic CVO fusion: product (WGS default) or min (WES)."); +// Sex-chromosome haploid calling (mirror of upstream postprocess_variants.py +// --haploid_contigs / --par_regions_bed). On a haploid contig (e.g. chrX/chrY +// in an XY sample) outside the pseudo-autosomal regions, heterozygous +// genotypes are disallowed: their probabilities are zeroed and the vector +// renormalized, forcing a haploid (homozygous) call. +ABSL_FLAG(std::string, haploid_contigs, "", + "Comma/space-separated contigs to call as haploid (e.g. " + "\"chrX,chrY\" for GRCh38, \"X,Y\" for GRCh37). Empty = all diploid."); +ABSL_FLAG(std::string, par_regions_bed, "", + "BED of pseudo-autosomal regions exempted from haploid calling on " + "the --haploid_contigs (they stay diploid). Empty = no exemptions."); namespace deepvariant { @@ -109,6 +123,27 @@ namespace { constexpr int kMaxPhred = 99; +// --------------------------------------------------------------------------- +// Sex-chromosome haploid calling (--haploid_contigs / --par_regions_bed). +// ParRegions / ParseHaploidContigs / LoadParRegions / IsHaploidPosition are +// shared with the make_examples gVCF path via haploid_regions.h. +// --------------------------------------------------------------------------- + +// True when the variant should be called haploid: it sits on a --haploid_contig +// and does not overlap a PAR region. Mirror of postprocess_variants.py's +// `is_non_autosome(v) and not is_in_regions(v, par_regions)`. +bool IsHaploidVariant(const Variant& variant, + const std::set& haploid_contigs, + const ParRegions& par_regions) { + // Upstream is_in_regions -> RangeSet.variant_overlaps tests only the single + // 0-based point variant.start (not the variant span), so a multi-base + // variant whose start is outside every PAR interval is still corrected even + // if its tail reaches into one. Match that exactly: probe [start, start+1). + const int64_t start = variant.start(); + return IsHaploidPosition(variant.reference_name(), start, start + 1, + haploid_contigs, par_regions); +} + std::vector ExpandShards(const std::string& spec) { auto at = spec.find('@'); if (at == std::string::npos) return {spec}; @@ -545,7 +580,11 @@ int RunPostprocessVariants(int argc, char** argv) { ? contig_to_pos.at(std::get<0>(b)) : INT_MAX; if (pa != pb) return pa < pb; - return std::get<1>(a) < std::get<1>(b); + if (std::get<1>(a) != std::get<1>(b)) + return std::get<1>(a) < std::get<1>(b); + if (std::get<2>(a) != std::get<2>(b)) + return std::get<2>(a) < std::get<2>(b); + return std::get<3>(a) < std::get<3>(b); }); // ── Open VCF writer ─────────────────────────────────────────────────────── @@ -565,6 +604,19 @@ int RunPostprocessVariants(int argc, char** argv) { const double qual_filter = absl::GetFlag(FLAGS_qual_filter); const double homref_min_gq = absl::GetFlag(FLAGS_cnn_homref_call_min_gq); + // Sex-chromosome haploid calling config. + const std::set haploid_contigs = + ParseHaploidContigs(absl::GetFlag(FLAGS_haploid_contigs)); + ParRegions par_regions; + if (const std::string par_bed = absl::GetFlag(FLAGS_par_regions_bed); + !par_bed.empty()) { + std::string err; + if (!LoadParRegions(par_bed, &par_regions, &err)) { + LOG(ERROR) << err; + return 1; + } + } + int written = 0; int refcall = 0; int nocall = 0; @@ -750,6 +802,14 @@ int RunPostprocessVariants(int argc, char** argv) { } like = std::move(like_pruned); + // Haploid correction: on a haploid contig outside the PAR, disallow + // heterozygous genotypes before the genotype/QUAL/GQ/GL are derived from + // `like` (mirrors merge_predictions applying + // correct_nonautosome_probabilities to the merged probabilities). + if (IsHaploidVariant(variant, haploid_contigs, par_regions)) { + CorrectNonautosomeProbabilities(&like, n_alts); + } + // argmax genotype. int best = 0; for (int i = 1; i < n_gt; ++i) { @@ -1013,9 +1073,19 @@ int RunPostprocessVariants(int argc, char** argv) { for (const auto& v : variants_buffer) { std::string serialized; v.SerializeToString(&serialized); - w->WriteRecord(serialized); + if (!w->WriteRecord(serialized)) { + LOG(ERROR) << "Failed to write temp variant TFRecord: " + << tmp_var_tfrecord; + std::remove(tmp_var_tfrecord.c_str()); // don't leave a partial file + return 1; + } + } + if (!w->Close()) { + LOG(ERROR) << "Failed to flush temp variant TFRecord: " + << tmp_var_tfrecord; + std::remove(tmp_var_tfrecord.c_str()); // don't leave a partial file + return 1; } - w->Close(); } // Open the variant + non-variant readers and a second VcfWriter diff --git a/deepvariant/native/regions.cc b/deepvariant/native/regions.cc index 8c5d2d2f..df35ce20 100644 --- a/deepvariant/native/regions.cc +++ b/deepvariant/native/regions.cc @@ -9,6 +9,7 @@ #include "absl/log/log.h" #include "absl/strings/numbers.h" +#include "absl/strings/str_replace.h" #include "absl/strings/str_split.h" #include "third_party/nucleus/protos/range.pb.h" #include "third_party/nucleus/protos/reference.pb.h" @@ -39,16 +40,23 @@ bool ParseRegionString( auto dash = range_part.find('-'); int64_t start_1based, end_1based; + // Strip comma grouping from numeric substrings (e.g. "1,000,000") before + // parsing, matching upstream ranges.parse_literal which accepts [0-9,]+. if (dash == std::string::npos) { // Single position: "chr1:1000" → half-open [999, 1000) - if (!absl::SimpleAtoi(range_part, &start_1based)) { + const std::string pos_str = absl::StrReplaceAll(range_part, {{",", ""}}); + if (!absl::SimpleAtoi(pos_str, &start_1based)) { LOG(ERROR) << "Cannot parse position in region: " << s; return false; } end_1based = start_1based; } else { - if (!absl::SimpleAtoi(range_part.substr(0, dash), &start_1based) || - !absl::SimpleAtoi(range_part.substr(dash + 1), &end_1based)) { + const std::string start_str = + absl::StrReplaceAll(range_part.substr(0, dash), {{",", ""}}); + const std::string end_str = + absl::StrReplaceAll(range_part.substr(dash + 1), {{",", ""}}); + if (!absl::SimpleAtoi(start_str, &start_1based) || + !absl::SimpleAtoi(end_str, &end_1based)) { LOG(ERROR) << "Cannot parse range in region: " << s; return false; } diff --git a/deepvariant/native/small_model_features.cc b/deepvariant/native/small_model_features.cc index 7912469b..09148f09 100644 --- a/deepvariant/native/small_model_features.cc +++ b/deepvariant/native/small_model_features.cc @@ -19,6 +19,17 @@ using nucleus::genomics::v1::Variant; namespace { +// Count alternate_bases entries that are not symbolic/placeholder alleles. +// Mirrors variant_utils.is_multiallelic, which default-excludes "<*>", +// "", and "." before counting. +int NumNonSymbolicAlts(const Variant& v) { + int n = 0; + for (const auto& alt : v.alternate_bases()) { + if (alt != "<*>" && alt != "" && alt != ".") ++n; + } + return n; +} + // Pull the reads supporting the chosen alt allele indices into a single // flat vector of ReadSupport pointers. If `sample_filter` is non-empty, // only reads with matching `sample_name` are retained (mirrors upstream's @@ -197,7 +208,7 @@ std::vector EncodeSmallModelFeatures( del_len = std::max(del_len, d); } features.push_back(std::max(0, del_len)); // deletion_length - features.push_back(v.alternate_bases_size() > 1 ? 1 : 0); // is_multiallelic + features.push_back(NumNonSymbolicAlts(v) > 1 ? 1 : 0); // is_multiallelic features.push_back(alt_allele_indices.size() > 1 ? 1 : 0); // is_multiple_alt_alleles // ── VAF context (51 features, offsets -25..+25 inclusive) ───────────────── @@ -271,7 +282,7 @@ std::vector EncodeSmallModelFeaturesMultiSample( del_len = std::max(del_len, d); } features.push_back(std::max(0, del_len)); - features.push_back(v.alternate_bases_size() > 1 ? 1 : 0); + features.push_back(NumNonSymbolicAlts(v) > 1 ? 1 : 0); features.push_back(alt_allele_indices.size() > 1 ? 1 : 0); // ── VAF context (51) ───────────────────────────────────────────────────── @@ -388,7 +399,7 @@ std::vector EncodeSmallModelFeaturesHaplotype( } features.push_back(std::max(0, ins_len)); features.push_back(std::max(0, del_len)); - features.push_back(v.alternate_bases_size() > 1 ? 1 : 0); + features.push_back(NumNonSymbolicAlts(v) > 1 ? 1 : 0); features.push_back(static_cast(alt_allele_indices.size()) > 1 ? 1 : 0); { const auto& vaf_at_pos = candidate.allele_frequency_at_position(); diff --git a/deepvariant/native/tfrecord.cc b/deepvariant/native/tfrecord.cc index 8459b83b..28905780 100644 --- a/deepvariant/native/tfrecord.cc +++ b/deepvariant/native/tfrecord.cc @@ -4,6 +4,7 @@ #include "deepvariant/native/tfrecord.h" #include +#include #include #include #include @@ -103,38 +104,78 @@ bool TFRecordReader::GetNext() { while (true) { auto& s = impl_->stream; if (s.good()) { + const std::string& path = impl_->paths[impl_->current_index]; uint64_t length = 0; s.read(reinterpret_cast(&length), 8); - if (s.gcount() == 8) { - s.seekg(4, std::ios::cur); // skip length CRC (not verified) + const std::streamsize len_read = s.gcount(); + if (len_read == 8) { + // Read and verify the length CRC (masked CRC32C over the 8 length + // bytes) rather than seekg-skipping it. + uint32_t len_crc = 0; + s.read(reinterpret_cast(&len_crc), 4); + if (s.gcount() != 4) { + std::fprintf(stderr, + "tfrecord: truncated TFRecord (length CRC) in %s at " + "offset %lld\n", + path.c_str(), static_cast(offset_)); + return false; + } + const uint32_t expected_len_crc = + MaskedCrc32c(reinterpret_cast(&length), sizeof(length)); + if (len_crc != expected_len_crc) { + std::fprintf(stderr, + "tfrecord: length CRC mismatch in %s at offset %lld\n", + path.c_str(), static_cast(offset_)); + return false; + } record_.resize(length); s.read(record_.data(), static_cast(length)); if (static_cast(s.gcount()) != length) { - // BUG FIX (2026-05-10): the previous `return false` here would - // ABANDON all remaining shards in a multi-shard read whenever - // the LAST record of any shard was truncated. On a 14-shard - // WG run this caused 13/14 shards (~95 % of examples) to be - // silently dropped: call_variants only saw 69k of 954k - // examples → 947k PASS calls missing in the final VCF. - // - // Truncation cause: upstream's ExamplesGenerator destructor - // closes the writer without an explicit flush — the last - // partial-buffer write (1 record per shard, ≈10-150 KB out - // of 1 MiB buffer) is dropped on close. - // - // Fix: treat partial-payload same as EOF — fall through to - // shard-advance code. Loses the 1 truncated record per shard - // (unrecoverable since it was never written to disk) but - // preserves all following shards. WG impact: 14 lost records - // out of 954k = 0.0015 % vs 100 % loss before the fix. - // Fall through to shard-advance code below. + // A partial payload read (0 < gcount < length) is genuine + // truncation, not a clean record boundary: surface it as an + // error instead of silently advancing to the next shard. + std::fprintf(stderr, + "tfrecord: truncated TFRecord payload in %s at offset " + "%lld: read %lld of %llu bytes\n", + path.c_str(), static_cast(offset_), + static_cast(s.gcount()), + static_cast(length)); + return false; } else { - s.seekg(4, std::ios::cur); // skip payload CRC + // Read and verify the payload CRC (masked CRC32C over record_). + uint32_t data_crc = 0; + s.read(reinterpret_cast(&data_crc), 4); + if (s.gcount() != 4) { + std::fprintf(stderr, + "tfrecord: truncated TFRecord (payload CRC) in %s at " + "offset %lld\n", + path.c_str(), static_cast(offset_)); + return false; + } + const uint32_t expected_data_crc = + MaskedCrc32c(record_.data(), record_.size()); + if (data_crc != expected_data_crc) { + std::fprintf(stderr, + "tfrecord: payload CRC mismatch in %s at offset " + "%lld\n", + path.c_str(), static_cast(offset_)); + return false; + } offset_ += 8 + 4 + length + 4; return true; } + } else if (len_read != 0) { + // A partial length read at a record boundary is truncation. + std::fprintf(stderr, + "tfrecord: truncated TFRecord (length field) in %s at " + "offset %lld: read %lld of 8 bytes\n", + path.c_str(), static_cast(offset_), + static_cast(len_read)); + return false; } + // len_read == 0: clean EOF at a record boundary — fall through to the + // shard-advance code below. } // Current shard exhausted (or read failed at boundary). Try next shard. if (impl_->current_index + 1 >= impl_->paths.size()) return false; @@ -167,6 +208,8 @@ void TFRecordReader::Close() { namespace { constexpr size_t kBufBytes = 1 << 20; // 1 MiB write coalescing buffer +static_assert(kBufBytes % 4096 == 0, + "F_NOCACHE requires sector-aligned buffer"); } struct TFRecordWriter::Impl { @@ -279,7 +322,10 @@ bool TFRecordWriter::Close() { if (!impl_) return true; bool ok = impl_->FlushBuf(); if (impl_->fd >= 0) { - ::close(impl_->fd); + // Best-effort durability on macOS: flush the device's write cache so the + // record is on stable storage before we report success. + ::fcntl(impl_->fd, F_FULLFSYNC, 0); + if (::close(impl_->fd) != 0) ok = false; impl_->fd = -1; } return ok; diff --git a/patches/tfrecord_reader_macos.cc b/patches/tfrecord_reader_macos.cc index d6ea5823..65768eff 100644 --- a/patches/tfrecord_reader_macos.cc +++ b/patches/tfrecord_reader_macos.cc @@ -1,7 +1,7 @@ // POSIX replacement for third_party/nucleus/io/tfrecord_reader.cc. // TFRecord format: [uint64_le length][uint32_le masked_crc32c(len)] // [bytes payload][uint32_le masked_crc32c(payload)] -// Uncompressed only. CRC not verified (dev-speed path). +// Uncompressed only. Masked CRC32Cs are verified against the data. #include "third_party/nucleus/io/tfrecord_reader.h" @@ -10,15 +10,29 @@ #include #include #include +#include #include +#include "absl/crc/crc32c.h" +#include "absl/log/log.h" + namespace nucleus { namespace { +// Mask delta and helper must match those used by the writer +// (tfrecord_writer_macos.cc) so written CRCs verify here. +constexpr uint32_t kMaskDelta = 0xa282ead8UL; +uint32_t MaskedCrc32c(const char* data, size_t n) { + uint32_t crc = static_cast( + absl::ComputeCrc32c(std::string_view(data, n))); + return ((crc >> 15) | (crc << 17)) + kMaskDelta; +} + struct TFRRImpl { std::ifstream stream; - explicit TFRRImpl(const std::string& path) - : stream(path, std::ios::binary) {} + std::string path; + explicit TFRRImpl(const std::string& p) + : stream(p, std::ios::binary), path(p) {} }; std::mutex mu; std::unordered_map> impls; @@ -48,17 +62,57 @@ bool TFRecordReader::GetNext() { auto it = impls.find(this); if (it == impls.end()) return false; auto& s = it->second->stream; + const std::string& path = it->second->path; if (!s.good()) return false; + // Read the 8-byte little-endian length. A zero-byte read at a record + // boundary is a clean EOF; a partial read (1..7 bytes) is truncation. uint64_t length = 0; s.read(reinterpret_cast(&length), 8); - if (s.gcount() != 8) return false; - s.seekg(4, std::ios::cur); // skip length CRC + std::streamsize length_read = s.gcount(); + if (length_read == 0) return false; // clean EOF + if (length_read != 8) { + LOG(ERROR) << "Truncated TFRecord (short length field, " << length_read + << "/8 bytes) in " << path; + return false; + } + // Read and verify the masked CRC32C of the 8-byte length field. + uint32_t length_crc = 0; + s.read(reinterpret_cast(&length_crc), 4); + if (s.gcount() != 4) { + LOG(ERROR) << "Truncated TFRecord (short length CRC) in " << path; + return false; + } + uint32_t expected_length_crc = + MaskedCrc32c(reinterpret_cast(&length), 8); + if (length_crc != expected_length_crc) { + LOG(ERROR) << "Corrupt TFRecord (length CRC mismatch) in " << path; + return false; + } + + // Read the payload. A short payload read is truncation, never clean EOF + // (we are mid-record after consuming a valid length field). record_.resize(length); s.read(record_.data(), static_cast(length)); - if (static_cast(s.gcount()) != length) return false; - s.seekg(4, std::ios::cur); // skip payload CRC + if (static_cast(s.gcount()) != length) { + LOG(ERROR) << "Truncated TFRecord (short payload, " << s.gcount() << "/" + << length << " bytes) in " << path; + return false; + } + + // Read and verify the masked CRC32C of the payload. + uint32_t data_crc = 0; + s.read(reinterpret_cast(&data_crc), 4); + if (s.gcount() != 4) { + LOG(ERROR) << "Truncated TFRecord (short payload CRC) in " << path; + return false; + } + uint32_t expected_data_crc = MaskedCrc32c(record_.data(), length); + if (data_crc != expected_data_crc) { + LOG(ERROR) << "Corrupt TFRecord (payload CRC mismatch) in " << path; + return false; + } offset_ += 8 + 4 + length + 4; return true; diff --git a/patches/tfrecord_writer_macos.cc b/patches/tfrecord_writer_macos.cc index ddd5319f..913a51ff 100644 --- a/patches/tfrecord_writer_macos.cc +++ b/patches/tfrecord_writer_macos.cc @@ -77,9 +77,20 @@ bool TFRecordWriter::Close() { std::lock_guard lk(mu); auto it = impls.find(this); if (it == impls.end()) return true; - it->second->stream.flush(); - it->second->stream.close(); - return !it->second->stream.fail(); + auto& s = it->second->stream; + // Flush user-space buffers and surface any write error before closing. + // std::ofstream exposes no portable file descriptor, so we cannot + // ::fsync() the underlying file to force a kernel-level durability + // barrier; the strongest guarantee available is that all bytes were + // successfully handed to the OS (good() after flush) and that close() + // itself did not fail. + s.flush(); + const bool flush_ok = s.good(); + // Always close to release the underlying stream/fd, even if flush detected + // a write error — otherwise a failed flush leaks the open stream in `impls`. + s.close(); + const bool close_ok = !s.fail(); + return flush_ok && close_ok; } } // namespace nucleus diff --git a/release/build_glnexus.sh b/release/build_glnexus.sh index 346e8671..8262e863 100755 --- a/release/build_glnexus.sh +++ b/release/build_glnexus.sh @@ -66,7 +66,8 @@ cd "${WORK}" if [ ! -f "GLnexus-${VERSION}.tar.gz" ]; then echo "==> Downloading GLnexus v${VERSION} ..." - curl -sL "${URL}" -o "GLnexus-${VERSION}.tar.gz" + curl -fsL --retry 3 --connect-timeout 15 "${URL}" -o "GLnexus-${VERSION}.tar.gz.partial" + mv "GLnexus-${VERSION}.tar.gz.partial" "GLnexus-${VERSION}.tar.gz" fi if [ ! -d "GLnexus-${VERSION}" ]; then diff --git a/scripts/build-prereq-macos.sh b/scripts/build-prereq-macos.sh index 82304ff1..dd208acc 100755 --- a/scripts/build-prereq-macos.sh +++ b/scripts/build-prereq-macos.sh @@ -59,6 +59,10 @@ echo " pyenv: $(pyenv --version)" echo " brew: $(brew --version | head -1)" echo "==> ready. next:" -echo " cmake -S . -B build -G Ninja" -echo " cmake --build build --parallel" -echo " ctest --test-dir build --output-on-failure" +# NB: use 'build-macos', not 'build' — the repo has a Bazel 'BUILD' file at the +# root and macOS's default case-insensitive filesystem treats a 'build/' +# directory as the same name, clobbering it. README/docs/release already +# standardize on 'build-macos'. +echo " cmake -S . -B build-macos -G Ninja" +echo " cmake --build build-macos --parallel" +echo " ctest --test-dir build-macos --output-on-failure" diff --git a/tests/native/CMakeLists.txt b/tests/native/CMakeLists.txt index 55e0094f..f490fa6e 100644 --- a/tests/native/CMakeLists.txt +++ b/tests/native/CMakeLists.txt @@ -24,6 +24,31 @@ function(phase1_test name src lib) add_test(NAME ${name} COMMAND ${name}) endfunction() +# Like phase1_test, but the source lives outside tests/native/ (it is the +# upstream test, shared verbatim). Pass an absolute path to the .cc. +function(phase1_test_abs name abs_src lib) + add_executable(${name} "${abs_src}") + target_include_directories(${name} PRIVATE ${GATE_INCLUDE_DIRS}) + target_compile_options(${name} PRIVATE ${GATE_COMPILE_OPTS}) + target_link_libraries(${name} PRIVATE ${lib} gtest_main gmock) + add_test(NAME ${name} COMMAND ${name}) +endfunction() + +# Upstream methylation_aware_phasing unit test, run verbatim against the +# native build of the algorithm (dv_methylation_aware_phasing). +phase1_test_abs(methylation_aware_phasing_test + "${CMAKE_SOURCE_DIR}/deepvariant/methylation_aware_phasing_test.cc" + dv_methylation_aware_phasing) + +# Header-only haploid-calling helpers (deepvariant/native/haploid_regions.h): +# needs only absl + gtest, no dv library. +add_executable(haploid_regions_test + "${CMAKE_CURRENT_SOURCE_DIR}/haploid_regions_test.cc") +target_include_directories(haploid_regions_test PRIVATE ${GATE_INCLUDE_DIRS}) +target_link_libraries(haploid_regions_test PRIVATE + absl::strings gtest_main gmock) +add_test(NAME haploid_regions_test COMMAND haploid_regions_test) + phase1_test(nucleus_io_smoke nucleus_io_smoke_test.cc nucleus_io) phase1_test(realigner_smoke realigner_smoke_test.cc realigner) phase1_test(call_variants_smoke call_variants_smoke_test.cc dv_call_variants_lib) @@ -31,3 +56,19 @@ phase1_test(small_model_smoke small_model_smoke_test.cc dv_small_model) phase1_test(dv_weights_smoke dv_weights_smoke_test.cc dv_weights) phase1_test(metal_inference_smoke metal_inference_smoke_test.cc dv_metal_inference) phase1_test(bnns_finalize_smoke bnns_finalize_smoke_test.cc dv_bnns_finalize) + +# Self-contained determinism / fidelity microtests (defined as executables in +# deepvariant/native/CMakeLists.txt). Each returns non-zero on mismatch, so it +# is a real ctest gate. The fixture-dependent microtests (det_mixed5b, +# det_inception, bnns_stem — they read /tmp/dv_per_layer/*.npy) are left +# unregistered since they need pre-generated per-layer dumps. +foreach(t + microtest_numpy_rng + microtest_neon_base_color + microtest_neon_cigar_classify + microtest_conv_serial + microtest_conv_kahan) + if(TARGET ${t}) + add_test(NAME ${t} COMMAND ${t}) + endif() +endforeach() diff --git a/tests/native/haploid_regions_test.cc b/tests/native/haploid_regions_test.cc new file mode 100644 index 00000000..e7bfdc45 --- /dev/null +++ b/tests/native/haploid_regions_test.cc @@ -0,0 +1,92 @@ +// Unit tests for the sex-chromosome haploid-calling helpers shared by the +// make_examples gVCF path and postprocess (deepvariant/native/haploid_regions.h). + +#include "deepvariant/native/haploid_regions.h" + +#include +#include +#include + +#include "gtest/gtest.h" + +namespace deepvariant { +namespace { + +TEST(HaploidRegions, ParseHaploidContigsSplitsCommaAndWhitespace) { + EXPECT_TRUE(ParseHaploidContigs("").empty()); + EXPECT_EQ(ParseHaploidContigs("chrX,chrY"), + (std::set{"chrX", "chrY"})); + // Upstream splits on ',' then on whitespace, so both separators work. + EXPECT_EQ(ParseHaploidContigs("X Y"), (std::set{"X", "Y"})); + EXPECT_EQ(ParseHaploidContigs(" chrX , chrY "), + (std::set{"chrX", "chrY"})); +} + +TEST(HaploidRegions, ParRegionsOverlapIsHalfOpen) { + ParRegions par; + par.by_contig["chrX"] = {{10, 20}, {100, 200}}; + EXPECT_TRUE(par.Overlaps("chrX", 10, 11)); // at start + EXPECT_TRUE(par.Overlaps("chrX", 19, 25)); // straddles end + EXPECT_FALSE(par.Overlaps("chrX", 20, 25)); // [20,25) abuts [10,20) — no + EXPECT_FALSE(par.Overlaps("chrX", 0, 10)); // [0,10) abuts [10,20) — no + EXPECT_FALSE(par.Overlaps("chrY", 10, 11)); // wrong contig +} + +TEST(HaploidRegions, IsHaploidPositionRespectsContigSetAndPar) { + const std::set haploid = {"chrX", "chrY"}; + ParRegions par; + par.by_contig["chrX"] = {{1000, 2000}}; + + // On a haploid contig, outside PAR → haploid. + EXPECT_TRUE(IsHaploidPosition("chrX", 500, 501, haploid, par)); + // On a haploid contig, inside PAR → diploid. + EXPECT_FALSE(IsHaploidPosition("chrX", 1500, 1501, haploid, par)); + // Autosome → never haploid. + EXPECT_FALSE(IsHaploidPosition("chr20", 500, 501, haploid, par)); + // Empty haploid set → never haploid. + EXPECT_FALSE(IsHaploidPosition("chrX", 500, 501, {}, par)); +} + +TEST(HaploidRegions, CorrectNonautosomeBiallelicZerosHet) { + // PL order for 1 alt: [0/0, 0/1, 1/1]; 0/1 (index 1) is het. + std::vector like = {0.2, 0.5, 0.3}; + CorrectNonautosomeProbabilities(&like, /*n_alts=*/1); + ASSERT_EQ(like.size(), 3u); + EXPECT_DOUBLE_EQ(like[1], 0.0); + EXPECT_DOUBLE_EQ(like[0], 0.4); // 0.2 / (0.2 + 0.3) + EXPECT_DOUBLE_EQ(like[2], 0.6); // 0.3 / (0.2 + 0.3) +} + +TEST(HaploidRegions, CorrectNonautosomeTriallelicZerosAllHets) { + // PL order for 2 alts: 0/0,0/1,1/1,0/2,1/2,2/2. + // Het genotypes: 0/1(1), 0/2(3), 1/2(4). Hom: 0/0(0), 1/1(2), 2/2(5). + std::vector like = {0.1, 0.2, 0.1, 0.2, 0.3, 0.1}; + CorrectNonautosomeProbabilities(&like, /*n_alts=*/2); + EXPECT_DOUBLE_EQ(like[1], 0.0); + EXPECT_DOUBLE_EQ(like[3], 0.0); + EXPECT_DOUBLE_EQ(like[4], 0.0); + // Surviving hom mass {0.1, 0.1, 0.1} renormalizes to thirds. + EXPECT_DOUBLE_EQ(like[0], 1.0 / 3.0); + EXPECT_DOUBLE_EQ(like[2], 1.0 / 3.0); + EXPECT_DOUBLE_EQ(like[5], 1.0 / 3.0); +} + +TEST(HaploidRegions, CorrectNonautosomeAllZeroIsSafe) { + // Degenerate input (all het, everything zeroed) must not divide by zero. + std::vector like = {0.0, 1.0, 0.0}; + CorrectNonautosomeProbabilities(&like, /*n_alts=*/1); + for (double v : like) EXPECT_DOUBLE_EQ(v, 0.0); +} + +TEST(HaploidRegions, CorrectNonautosomeSizeMismatchIsNoOp) { + // A vector whose length doesn't match (n_alts+1)(n_alts+2)/2 is a contract + // violation; the function must leave it untouched rather than write OOB. + std::vector like = {0.2, 0.5}; // n_alts=1 expects 3 entries. + CorrectNonautosomeProbabilities(&like, /*n_alts=*/1); + ASSERT_EQ(like.size(), 2u); + EXPECT_DOUBLE_EQ(like[0], 0.2); + EXPECT_DOUBLE_EQ(like[1], 0.5); +} + +} // namespace +} // namespace deepvariant diff --git a/third_party/tf_example b/third_party/tf_example deleted file mode 120000 index 982d0421..00000000 --- a/third_party/tf_example +++ /dev/null @@ -1 +0,0 @@ -/Users/benjamin/deepvariant/tools/conversion/Protos/tensorflow \ No newline at end of file diff --git a/tools/conversion/README.md b/tools/conversion/README.md index 6b4f7046..d7493dcd 100644 --- a/tools/conversion/README.md +++ b/tools/conversion/README.md @@ -1,21 +1,21 @@ # Phase 0 — Inference framework bench (dev-time Python tooling) -> **Dev-time only — never shipped to users.** Python lives here because `coremltools`, `tf2onnx`, `tensorflow-metal`, and the SavedModel reader for MLX are all Python-only Apple/Google packages. Re-implementing them in Swift would re-introduce 7+ years of edge-case bug fixes (BatchNorm fused vs not, ConcatV2 axis handling, NHWC↔NCHW layout, FusedBatchNormV3 epsilon, etc.). Voie B accepts dev-time Python to keep this maturity, while the user-facing `deepvariant` binary stays 100 % C++/Obj-C++ with no embedded Python interpreter (verified by `otool -L` in Phase 5). See `CLAUDE.md` for the dev-time vs runtime split. +> **Dev-time only — never shipped to users.** Python lives here because `coremltools` and `mlx` are Python-only Apple packages, and coremltools' MIL builder carries years of edge-case maturity (BatchNorm fused vs not, ConcatV2 axis handling, NHWC↔NCHW layout, FusedBatchNormV3 epsilon, etc.) that we don't want to re-implement in Swift. The conversion pipeline is **TF-free**: weights are read straight from the SavedModel's TensorBundle by the pure-protobuf readers in this directory (`savedmodel_reader.py`, `tensor_bundle_reader.py`), with no TensorFlow, PyTorch, ONNX, or tensorflow-metal involved. The user-facing `deepvariant` binary stays 100 % C++/Obj-C++ with no embedded Python interpreter (verified by `otool -L` in Phase 5). See `CLAUDE.md` for the dev-time vs runtime split. -Three-way A/B/C bench (tensorflow-metal vs Core ML vs MLX) on the real WGS SavedModel: produces latency, throughput, GPU/ANE residency, and parity-vs-Linux measurements that feed the Phase 0 ADR. +Two-way A/B bench (Core ML vs MLX) on the real WGS SavedModel: produces latency, throughput, GPU/ANE residency, and parity-vs-Linux measurements that feed the Phase 0 ADR. (The earlier tensorflow-metal voie was dropped — tensorflow-metal 1.2.0 has been unmaintained since mid-2024.) ## Layout ```text tools/conversion/ ├── .python-version # pyenv pin: 3.11.x -├── requirements-coreml.txt # TF 2.16.2 + coremltools 7.2 -├── requirements-metal.txt # TF 2.16.2 + tensorflow-metal 1.2.0 -├── requirements-mlx.txt # MLX 0.21+ (TF only for SavedModel weight ingest) -├── setup_venvs.sh # creates venv-coreml / venv-metal / venv-mlx +├── requirements-coreml.txt # coremltools 9.0+ (TF-free MIL path) +├── requirements-mlx.txt # MLX 0.21+ + safetensors (TF-free) +├── setup_venvs.sh # creates venv-coreml / venv-mlx (TF-free) ├── fetch_savedmodel.sh # pulls gs://deepvariant/models/DeepVariant/1.10.0/ -├── convert_coreml.py # SavedModel -> .mlpackage (coremltools, with tf2onnx fallback) -├── convert_metal.py # passthrough (records TF/metal env metadata) +├── savedmodel_reader.py # pure-protobuf SavedModel reader (no TF) +├── tensor_bundle_reader.py # pure-protobuf TensorBundle weight reader (no TF) +├── convert_coreml.py # TensorBundle weights -> .mlpackage (coremltools MIL) ├── convert_mlx.py # SavedModel weights -> safetensors (MLX-friendly) ├── bench.py # latency + powermetrics GPU residency + softmax capture ├── parity_check.py # softmax max-abs / argmax disagreement vs reference @@ -25,7 +25,7 @@ tools/conversion/ ## Pipeline ```sh -# one-time setup (~30 min: pyenv install 3.11 + 3 venv pip installs) +# one-time setup (pyenv install 3.11 + two venv pip installs) ./setup_venvs.sh # pull the real WGS SavedModel (~700 MB) @@ -33,11 +33,7 @@ tools/conversion/ # convert each voie source venv-coreml/bin/activate -python convert_coreml.py --saved-model models/wgs --output models/wgs.mlpackage -deactivate - -source venv-metal/bin/activate -python convert_metal.py --saved-model models/wgs --output models/wgs.metal.json +python convert_coreml.py --bundle models/wgs/variables/variables --output models/wgs.mlpackage deactivate source venv-mlx/bin/activate @@ -51,30 +47,27 @@ python bench.py --backend coreml --model models/wgs.mlpackage \ --output ../../benchmarks/coreml_wgs.json \ --output-cv ../../benchmarks/coreml_wgs.cv.tfrecord deactivate -# (repeat for metal, mlx) +# (repeat for mlx) # parity vs Linux reference python parity_check.py \ --reference ../reference/output/wgs/call_variants_chr20.tfrecord \ --candidates ../../benchmarks/coreml_wgs.cv.tfrecord \ - ../../benchmarks/metal_wgs.cv.tfrecord \ ../../benchmarks/mlx_wgs.cv.tfrecord ``` -## Why three pinned venvs +## Why two pinned venvs -These three frameworks have **incompatible** Python/TF/numpy version requirements: +Core ML and MLX have incompatible Python dependency requirements, so each gets its own venv. Both are **TF-free** — `setup_venvs.sh` hard-fails if `tensorflow` is importable in either. -| venv | Python | TF | Other | -| --- | --- | --- | --- | -| venv-coreml | 3.11 | 2.16.2 | coremltools 7.2 (v1 found 9.0 + TF 2.20 hangs at 21 min / 3 GB RSS) | -| venv-metal | 3.11 | 2.16.2 | tensorflow-metal 1.2.0 (frozen at TF 2.16 since mid-2024) | -| venv-mlx | 3.11 | 2.16.2 (read-only for SavedModel ingest) | MLX 0.21+ | +| venv | Python | Key deps | +| --- | --- | --- | +| venv-coreml | 3.11 | coremltools 9.0+ (direct MIL path; no TF SavedModel conversion) | +| venv-mlx | 3.11 | MLX 0.21+, safetensors (Inception-v3 rebuilt in MLX at bench time) | -We pin `numpy < 2` everywhere because TF 2.16 was built against numpy 1. +Both venvs pin `numpy < 2` (see the requirements files). ## Stop conditions -- If `convert_coreml.py` hangs > 5 min on the WGS SavedModel, bail and try `convert_coreml.py --via-onnx` (fallback path: SavedModel → ONNX → Core ML via tf2onnx). - If the parity check shows argmax disagreement on **any** of the 1000 examples, the framework is rejected — no exceptions. - If `bench.py` reports `gpu_power=0` AND `ane_power=0` for a backend, that backend has a config bug, not a perf result. diff --git a/tools/conversion/convert_via_docker.sh b/tools/conversion/convert_via_docker.sh index 8277282b..357f581e 100755 --- a/tools/conversion/convert_via_docker.sh +++ b/tools/conversion/convert_via_docker.sh @@ -27,6 +27,7 @@ import json, sys d = json.load(open('${MODEL_DIR}/model.example_info.json')) print(','.join(str(x) for x in d['shape'])) ") + [[ -n "${SHAPE_JSON}" && "${SHAPE_JSON}" == *,*,* ]] || { echo "error: failed to parse shape from ${MODEL_DIR}/model.example_info.json" >&2; exit 1; } else SHAPE_JSON="100,221,7" fi diff --git a/tools/conversion/tensor_bundle_reader.py b/tools/conversion/tensor_bundle_reader.py index 0fe09358..8ff1dbc2 100644 --- a/tools/conversion/tensor_bundle_reader.py +++ b/tools/conversion/tensor_bundle_reader.py @@ -223,7 +223,7 @@ def _iterate_block(block: bytes): 3: (np.int32, 4), # DT_INT32 9: (np.int64, 8), # DT_INT64 10: (np.bool_, 1), # DT_BOOL - 14: (np.float32, 2), # DT_BFLOAT16 (decoded as fp32; bfloat16 isn't in numpy) + 14: (np.float32, 2), # DT_BFLOAT16 (2 bytes on disk; widened to fp32 on read since bfloat16 isn't in numpy) 17: (np.float16, 2), # DT_HALF } @@ -327,7 +327,7 @@ def read_tensor(self, name: str) -> np.ndarray: raise NotImplementedError( f"dtype {e.dtype} not supported (variable {name})" ) - np_dtype, _ = _NUMPY_DTYPE[e.dtype] + np_dtype, dtype_size = _NUMPY_DTYPE[e.dtype] path = self.shard_path(e.shard_id) with open(path, "rb") as f: f.seek(e.offset) @@ -337,7 +337,21 @@ def read_tensor(self, name: str) -> np.ndarray: f"truncated read of {name} from {path}: " f"want {e.size} bytes, got {len(raw)}" ) - arr = np.frombuffer(raw, dtype=np_dtype) + # Guard: the raw byte count must be a whole number of elements that + # matches the declared shape. Otherwise np.frombuffer would silently + # mis-decode (e.g. reinterpreting 2-byte values as 4-byte ones). + expected_bytes = int(np.prod(e.shape)) * dtype_size + if len(raw) != expected_bytes: + raise ValueError( + f"size mismatch for variable {name}: shape {e.shape} with " + f"{dtype_size}-byte elements implies {expected_bytes} bytes, " + f"but got {len(raw)} bytes" + ) + if e.dtype == 14: # DT_BFLOAT16: 2 raw bytes are the high 16 bits of fp32. + u16 = np.frombuffer(raw, dtype=np.uint16) + arr = (u16.astype(np.uint32) << 16).view(np.float32) + else: + arr = np.frombuffer(raw, dtype=np_dtype) if e.shape: arr = arr.reshape(e.shape) return arr diff --git a/tools/reference/capture_linux_x86.sh b/tools/reference/capture_linux_x86.sh index e7aea368..a761da92 100755 --- a/tools/reference/capture_linux_x86.sh +++ b/tools/reference/capture_linux_x86.sh @@ -67,6 +67,8 @@ if (( ${#EXAMPLES_GLOB[@]} > 0 && ${#CV_GLOB[@]} > 0 )); then cp "${CV_GLOB[0]}" "${OUT}/call_variants_chr20.tfrecord" fi +[[ -s "${OUT}/examples_chr20.tfrecord" ]] || { echo "error: examples_chr20.tfrecord missing — make_examples produced no output" >&2; exit 1; } + # Build a 1000-example slice for fast bench iteration. mkdir -p cache python3 - < Downloading ${out} (URL: ${url})" - curl -L --progress-bar -o "${out}" "${url}" + curl -fL --retry 3 --connect-timeout 15 --progress-bar -o "${out}.partial" "${url}" + mv "${out}.partial" "${out}" } # 1. Full GRCh38 reference FASTA from NCBI canonical no_alt_analysis_set. diff --git a/validation/download_giab_strats.sh b/validation/download_giab_strats.sh index 917c164b..33bc87b3 100755 --- a/validation/download_giab_strats.sh +++ b/validation/download_giab_strats.sh @@ -45,7 +45,8 @@ if [ ! -f "${ARCHIVE}" ]; then echo "==> Downloading GIAB stratifications v3.6 GRCh38 (~1.4 GB) ..." echo " URL: ${URL}" echo " Target: ${TARGET}/${ARCHIVE}" - curl -L --progress-bar -o "${ARCHIVE}" "${URL}" + curl -fL --retry 3 --connect-timeout 15 --progress-bar -o "${ARCHIVE}.partial" "${URL}" + mv "${ARCHIVE}.partial" "${ARCHIVE}" fi echo "==> Extracting ..." diff --git a/validation/run_hg002_wg_compare.sh b/validation/run_hg002_wg_compare.sh index bc870691..b2d5e678 100755 --- a/validation/run_hg002_wg_compare.sh +++ b/validation/run_hg002_wg_compare.sh @@ -49,6 +49,7 @@ echo "[$(date +%H:%M:%S)] Docker baseline PID=${DOCKER_PID}" echo "[$(date +%H:%M:%S)] Waiting for Docker WG (${DOCKER}) ..." until [ -f "${DOCKER}" ] && [ -s "${DOCKER}" ]; do + kill -0 "${DOCKER_PID}" 2>/dev/null || { echo "ERROR: Docker baseline (PID ${DOCKER_PID}) exited before producing ${DOCKER}; see /tmp/hg002_wg_docker.log" >&2; exit 1; } sleep 600 # 10 min poll — Docker takes ~22 h, no rush done echo "[$(date +%H:%M:%S)] Docker WG ready"