From 750e344ed62ed79a8322a5786b580a4fafe7e258 Mon Sep 17 00:00:00 2001 From: Tim Dunn Date: Sat, 8 Aug 2026 00:24:51 -0400 Subject: [PATCH 1/2] refactor(phase): write the summary VCF with htslib instead of fprintf (#229) The summary VCF was emitted with fprintf, so every VCF semantic on the output side had to be re-implemented on strings. Build each record as a bcf1_t and write it with bcf_write instead, so htslib owns the encoding. write_summary_vcf() now opens the file with hts_open and writes a header built by summary_vcf_header(); print_var_info/print_var_empty/ print_var_sample become set_var_record(), which sets the fixed fields, and var_sample_fields()/empty_sample_fields(), which return one sample's FORMAT values for set_record_samples() to write. Per-allele values carry htslib's missing sentinels, and a sample of lower ploidy is padded with the end-of-vector marker rather than a shorter rendered list. The output changes in one visible way: BC is a Float, so htslib renders it compactly ("1" and "0.8", not "1.000000" and "0.800000"). Over the 72,585 records of the committed chr20 fixture that is the only field that differs; every other field is byte-identical and every BC value is unchanged at float32 precision. The header also declares PASS before fileDate, since bcf_hdr_init emits it first. Test assertions over the rendered floats are regenerated accordingly. --- src/phase.cpp | 81 +++---- src/variant.cpp | 314 ++++++++++++++++++++----- src/variant.h | 53 ++++- tests/integration/test-integration.yml | 32 +-- tests/unit/src/test_phase.cpp | 51 +++- 5 files changed, 390 insertions(+), 141 deletions(-) diff --git a/src/phase.cpp b/src/phase.cpp index 9280d4b..e438b84 100644 --- a/src/phase.cpp +++ b/src/phase.cpp @@ -24,50 +24,23 @@ * @note Contigs called by only one callset are included; a contig with no query variants has no * phase block to read, so its truth records are written with PB and BS defaulted * @throws ERROR if the output summary VCF file cannot be opened for writing + * @throws ERROR if the header or any record cannot be written * @throws ERROR if neither callset is selected next while variants remain */ void phaseblockData::write_summary_vcf(std::string out_vcf_fn) { // VCF header if (g.verbosity >= 1) INFO(" Writing summary VCF to '%s'", out_vcf_fn.data()); - FILE* out_vcf = fopen(out_vcf_fn.data(), "w"); + htsFile* out_vcf = hts_open(out_vcf_fn.data(), "w"); if (out_vcf == NULL) { ERROR("Failed to open summary VCF file '%s'", out_vcf_fn.data()); } - const std::chrono::time_point now{std::chrono::system_clock::now()}; - time_t tt = std::chrono::system_clock::to_time_t(now); - tm local_time = *localtime(&tt); - fprintf(out_vcf, "##fileformat=VCFv4.2\n"); - fprintf(out_vcf, "##fileDate=%04d%02d%02d\n", local_time.tm_year + 1900, - local_time.tm_mon + 1, local_time.tm_mday); - fprintf(out_vcf, "##CL=%s\n", g.cmd.data()); - for (size_t i = 0; i < this->contigs.size(); i++) { - fprintf(out_vcf, "##contig=\n", - this->contigs[i].data(), this->lengths[i]); + bcf_hdr_t* hdr = summary_vcf_header(this->contigs, this->lengths); + if (bcf_hdr_write(out_vcf, hdr) != 0) { + ERROR("Failed to write summary VCF header to '%s'", out_vcf_fn.data()); } - fprintf(out_vcf, "##FILTER=\n"); - // The per-haplotype fields carry one value per allele of the sample's GT, which is what VCF - // 4.4's Number=P declares. BCF_VL_P only reaches htslib in 1.23, so a consumer on any older - // bcftools or pysam would report a cardinality error; Number=. produces byte-identical records - // and merely gives up the declared cardinality, so the count and order are stated here instead. - const std::string per_allele = " One value per allele of this sample's GT, in GT allele order, " - "'.' for a reference allele."; - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT=\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT=\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT=\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT= 0, else '.').%s\">\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n", per_allele.data()); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT=\n"); - fprintf(out_vcf, "##FORMAT= 1/1 query, '-' if 1/1 truth -> 0/1 query, '.' otherwise)\">\n"); - fprintf(out_vcf, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tTRUTH\tQUERY\n"); + bcf1_t* rec = bcf_init(); + if (rec == NULL) ERROR("Failed to allocate summary VCF record"); // write variants for (std::string ctg : this->contigs) { @@ -133,7 +106,7 @@ void phaseblockData::write_summary_vcf(std::string out_vcf_fn) { phase_block++; /* if (next[QUERY] && ptrs[QUERY] < qvars->n) { */ - /* fprintf(out_vcf, "orig_gt: %s\tmatched_gt: %s\tcredit: %.2f|%.2f\tref_dist: %d|%d\tphase: %s\n", */ + /* fprintf(stderr, "orig_gt: %s\tmatched_gt: %s\tcredit: %.2f|%.2f\tref_dist: %d|%d\tphase: %s\n", */ /* gt_strs[vars[QUERY]->orig_gts[ptrs[QUERY]]].data(), */ /* gt_strs[vars[QUERY]->matched_gts[ptrs[QUERY]]].data(), */ /* vars[QUERY]->credit[HAP1][ptrs[QUERY]], */ @@ -144,39 +117,47 @@ void phaseblockData::write_summary_vcf(std::string out_vcf_fn) { /* ); */ /* } */ /* if (next[TRUTH] && ptrs[TRUTH] < tvars->n) { */ - /* fprintf(out_vcf, "orig_gt: %s\n", */ + /* fprintf(stderr, "orig_gt: %s\n", */ /* gt_strs[vars[TRUTH]->orig_gts[ptrs[TRUTH]]].data()); */ /* } */ + bcf_clear(rec); if (next[QUERY]) { // a positional tie between differing alleles is not a match: the two are written as // co-located records, the truth one on the next pass bool matched = next[TRUTH] && vars[QUERY]->refs[ptrs[QUERY]] == vars[TRUTH]->refs[ptrs[TRUTH]] && vars[QUERY]->alts[ptrs[QUERY]] == vars[TRUTH]->alts[ptrs[TRUTH]]; - vars[QUERY]->print_var_info(out_vcf, this->ref, ctg, ptrs[QUERY]); - if (matched) { - vars[TRUTH]->print_var_sample(out_vcf, ptrs[TRUTH], - sc_idx, phase_block, block_state == PHASE_SWAP, flip_error); - } else { - vars[TRUTH]->print_var_empty(out_vcf, sc_idx, phase_block); - } - vars[QUERY]->print_var_sample(out_vcf, ptrs[QUERY], - sc_idx, phase_block, block_state == PHASE_SWAP, flip_error, true); + vars[QUERY]->set_var_record(hdr, rec, this->ref, ctg, ptrs[QUERY]); + set_record_samples(hdr, rec, + matched ? vars[TRUTH]->var_sample_fields(ptrs[TRUTH], sc_idx, phase_block, + block_state == PHASE_SWAP, flip_error) : + empty_sample_fields(sc_idx, phase_block), + vars[QUERY]->var_sample_fields(ptrs[QUERY], sc_idx, phase_block, + block_state == PHASE_SWAP, flip_error, true)); ptrs[QUERY]++; if (matched) ptrs[TRUTH]++; } else if (next[TRUTH]) { - vars[TRUTH]->print_var_info(out_vcf, this->ref, ctg, ptrs[TRUTH]); - vars[TRUTH]->print_var_sample(out_vcf, ptrs[TRUTH], - sc_idx, phase_block, block_state == PHASE_SWAP, flip_error); - vars[QUERY]->print_var_empty(out_vcf, sc_idx, phase_block, true); + vars[TRUTH]->set_var_record(hdr, rec, this->ref, ctg, ptrs[TRUTH]); + set_record_samples(hdr, rec, + vars[TRUTH]->var_sample_fields(ptrs[TRUTH], sc_idx, phase_block, + block_state == PHASE_SWAP, flip_error), + empty_sample_fields(sc_idx, phase_block)); ptrs[TRUTH]++; } else { ERROR("No variants are selected next."); } + if (bcf_write(out_vcf, hdr, rec) != 0) { + ERROR("Failed to write summary VCF record at %s:%lld", ctg.data(), + static_cast(rec->pos+1)); + } } } - fclose(out_vcf); + bcf_destroy(rec); + bcf_hdr_destroy(hdr); + if (hts_close(out_vcf) != 0) { + ERROR("Failed to close summary VCF file '%s'", out_vcf_fn.data()); + } } diff --git a/src/variant.cpp b/src/variant.cpp index 7c53d46..d436019 100644 --- a/src/variant.cpp +++ b/src/variant.cpp @@ -3,6 +3,7 @@ * @brief Per-contig and per-callset variant containers with VCF parsing and output utilities. */ #include +#include #include #include #include @@ -334,100 +335,205 @@ void ctgVariants::set_var_matched_gt_on_hap(int var_idx, hap_t hap, bool set, /**************************************************************************************************/ /** - * @brief Writes fixed VCF fields (CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT) for one variant. - * @param[in] out_fp Open file pointer to output VCF + * @brief Builds the summary VCF header, declaring every FORMAT field and the TRUTH/QUERY samples. + * @param[in] contigs Contig names, in the order records are written + * @param[in] lengths Contig lengths, parallel to contigs + * @return Header owning its own memory, to be released by the caller with bcf_hdr_destroy() + * @throws ERROR The header cannot be allocated, a header line htslib rejects, a sample htslib + * rejects, or a header htslib cannot synchronize + */ +bcf_hdr_t* summary_vcf_header(const std::vector & contigs, + const std::vector & lengths) { + + // bcf_hdr_init() supplies the ##fileformat line + bcf_hdr_t* hdr = bcf_hdr_init("w"); + if (hdr == NULL) ERROR("Failed to allocate summary VCF header"); + + const std::chrono::time_point now{std::chrono::system_clock::now()}; + time_t tt = std::chrono::system_clock::to_time_t(now); + tm local_time = *localtime(&tt); + char file_date[32]; + snprintf(file_date, sizeof(file_date), "##fileDate=%04d%02d%02d", local_time.tm_year + 1900, + local_time.tm_mon + 1, local_time.tm_mday); + + // The per-haplotype fields carry one value per allele of the sample's GT, which is what VCF + // 4.4's Number=P declares. BCF_VL_P only reaches htslib in 1.23, so a consumer on any older + // bcftools or pysam would report a cardinality error; Number=. produces identical records and + // merely gives up the declared cardinality, so the count and order are stated here instead. + const std::string per_allele = " One value per allele of this sample's GT, in GT allele order, " + "'.' for a reference allele."; + std::vector lines = {file_date, "##CL=" + g.cmd}; + for (size_t i = 0; i < contigs.size(); i++) { + lines.push_back("##contig="); + } + // every record PASSes, and htslib rejects a filter its header does not declare; bcf_hdr_init() + // declares PASS itself, and appending an ID it already holds is a no-op rather than a duplicate + lines.push_back("##FILTER="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT= 0, else '.')." + per_allele + "\">"); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT="); + lines.push_back("##FORMAT= 1/1 query, '-' if 1/1 truth -> 0/1 query, '.' otherwise)\">"); + for (const std::string & line : lines) { + if (bcf_hdr_append(hdr, line.data()) != 0) + ERROR("Failed to add summary VCF header line '%s'", line.data()); + } + + // the samples are added in the order their columns are written + for (const char* sample : {"TRUTH", "QUERY"}) { + if (bcf_hdr_add_sample(hdr, sample) < 0) + ERROR("Failed to add sample '%s' to summary VCF header", sample); + } + if (bcf_hdr_sync(hdr) < 0) ERROR("Failed to synchronize summary VCF header"); + return hdr; +} + + +/** + * @brief Sets the fixed VCF fields (CHROM, POS, ID, REF, ALT, QUAL, FILTER) of one record. + * @param[in] hdr Summary VCF header, which must declare this contig + * @param[in,out] rec Cleared record to fill * @param[in] ref Reference FASTA data for retrieving flanking bases for indels * @param[in] ctg Contig name * @param[in] idx Variant index in this container + * @throws ERROR The contig is not declared in the header * @throws ERROR An INS/DEL sits at the contig start (0-based pos 0), leaving no preceding base to anchor * @throws ERROR The variant type is not TYPE_SUB, TYPE_INS, or TYPE_DEL + * @throws ERROR htslib rejects the record's FILTER or alleles */ -void ctgVariants::print_var_info(FILE* out_fp, std::shared_ptr ref, - const std::string & ctg, int idx) { - char ref_base; +void ctgVariants::set_var_record(const bcf_hdr_t* hdr, bcf1_t* rec, + std::shared_ptr ref, const std::string & ctg, int idx) const { + + rec->rid = bcf_hdr_name2id(hdr, ctg.data()); + if (rec->rid < 0) ERROR("Contig '%s' is not declared in the summary VCF header", ctg.data()); + bcf_float_set_missing(rec->qual); + if (bcf_add_filter(hdr, rec, bcf_hdr_id2int(hdr, BCF_DT_ID, "PASS")) < 0) + ERROR("Failed to set FILTER on summary VCF record at %s:%d", ctg.data(), this->poss[idx]); + + std::string ref_allele, alt_allele; switch (this->types[idx]) { case TYPE_SUB: - fprintf(out_fp, "%s\t%d\t.\t%s\t%s\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE", - ctg.data(), this->poss[idx]+1, this->refs[idx].data(), - this->alts[idx].data()); + rec->pos = this->poss[idx]; + ref_allele = this->refs[idx]; + alt_allele = this->alts[idx]; break; case TYPE_INS: - case TYPE_DEL: + case TYPE_DEL: { // INS/DEL are left-anchored on the preceding reference base; at contig start (0-based // pos 0) there is no preceding base, so guard against the out-of-bounds read of index -1 if (this->poss[idx] == 0) - ERROR("Cannot left-anchor INS/DEL at contig start (0-based pos 0) on '%s' in print_var_info", + ERROR("Cannot left-anchor INS/DEL at contig start (0-based pos 0) on '%s' in set_var_record", ctg.data()); - ref_base = ref->fasta.at(ctg)[this->poss[idx]-1]; - fprintf(out_fp, "%s\t%d\t.\t%s\t%s\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE", ctg.data(), - this->poss[idx], (ref_base + this->refs[idx]).data(), - (ref_base + this->alts[idx]).data()); + rec->pos = this->poss[idx] - 1; + char ref_base = ref->fasta.at(ctg)[this->poss[idx]-1]; + ref_allele = ref_base + this->refs[idx]; + alt_allele = ref_base + this->alts[idx]; break; + } default: - ERROR("print_var_info not implemented for type %d", static_cast(this->types[idx])); + ERROR("set_var_record not implemented for type %d", static_cast(this->types[idx])); } + + const char* alleles[2] = {ref_allele.data(), alt_allele.data()}; + if (bcf_update_alleles(hdr, rec, alleles, 2) < 0) + ERROR("Failed to set alleles on summary VCF record at %s:%d", ctg.data(), this->poss[idx]); } /** - * @brief Writes dot-separated empty sample fields for a variant with no call on this haplotype. - * @param[in] out_fp Open file pointer to output VCF - * @param[in] sc_idx Supercluster index for SC field - * @param[in] phase_block Phase block index for PB field - * @param[in] query If true, append newline (end of record); if false, tab (more samples follow) + * @brief Encodes the GT a sample reports for one variant, as the caller itself genotyped it. + * @param[in] orig_gt The caller's own genotype, never vcfdist's recovered matched_gt + * @param[in] ploidy Variant ploidy + * @return Phased allele indices: one for a haploid call, otherwise the diploid pair */ -void ctgVariants::print_var_empty(FILE* out_fp, int sc_idx, - int phase_block, bool query /* = false */) { - fprintf(out_fp, "\t.:.:.:.:.:.:.:%d:.:.:%d:.:.:.:.%s", sc_idx, phase_block, query ? "\n" : ""); +static std::vector genotype_alleles(gt_t orig_gt, ploidy_t ploidy) { + if (ploidy == PLOIDY_HAPLOID) return {bcf_gt_phased(1)}; + return {bcf_gt_phased(orig_gt == GT_ALT_REF || orig_gt == GT_ALT_ALT ? 1 : 0), + bcf_gt_phased(orig_gt == GT_REF_ALT || orig_gt == GT_ALT_ALT ? 1 : 0)}; } /** - * @brief Renders the GT a sample reports for one variant, as the caller itself genotyped it. - * @param[in] orig_gt The caller's own genotype, never vcfdist's recovered matched_gt - * @param[in] ploidy Variant ploidy - * @return "1" for a haploid call, otherwise the phased diploid pair + * @brief Returns the FORMAT values of a sample that made no call at a locus. + * @param[in] sc_idx Supercluster index for the SC field + * @param[in] phase_block Phase block index for the PB field + * @return Values reporting '.' for every field but SC and PB, which are locus-wide */ -static std::string display_gt(gt_t orig_gt, ploidy_t ploidy) { - return ploidy == PLOIDY_HAPLOID ? "1" : gt_strs[orig_gt]; +sample_fields empty_sample_fields(int sc_idx, int phase_block) { + sample_fields fields; + bcf_float_set_missing(fields.qq); + fields.sc = sc_idx; + fields.ps = bcf_int32_missing; + fields.pb = phase_block; + fields.bs = bcf_int32_missing; + fields.vp = bcf_int32_missing; + fields.fe = bcf_int32_missing; + return fields; } /** - * @brief Writes sample-specific FORMAT fields for one variant to output VCF. + * @brief Returns one sample's FORMAT values for a variant it called. * * One record is written per variant rather than per haplotype, so the per-haplotype fields (BD, BC, - * RD, QD, BK, SG) are comma-separated lists holding one value per haplotype: two for a diploid - * record, one for a haploid one. GT is rendered from orig_gt, the caller's own claim, so a - * haplotype carrying the reference allele has no evaluation data and every per-haplotype field - * reports "." for it. The evaluation lanes are keyed by matched_gt's haplotypes, which - * matched_gt_is_swapped() reports may be the reverse of orig_gt's. - * @param[in] out_fp Open file pointer to output VCF + * RD, QD, BK, SG) hold one value per haplotype: two for a diploid record, one for a haploid one. GT + * reports orig_gt, the caller's own claim, so a haplotype carrying the reference allele has no + * evaluation data and every per-haplotype field reports '.' for it. The evaluation lanes are keyed + * by matched_gt's haplotypes, which matched_gt_is_swapped() reports may be the reverse of orig_gt's. * @param[in] vi Variant index in this container - * @param[in] sc_idx Supercluster index for SC field - * @param[in] phase_block Phase block index for PB field + * @param[in] sc_idx Supercluster index for the SC field + * @param[in] phase_block Phase block index for the PB field * @param[in] phase_switch True if phase switched at this position * @param[in] phase_flip True if phase flipped (error) at this position - * @param[in] query If true, format as query sample; if false, as truth sample + * @param[in] query If true, report as the query sample; if false, as the truth sample + * @return This sample's FORMAT values, htslib-encoded */ -void ctgVariants::print_var_sample(FILE* out_fp, int vi, int sc_idx, int phase_block, - bool phase_switch, bool phase_flip, bool query /* = false */) { +sample_fields ctgVariants::var_sample_fields(int vi, int sc_idx, int phase_block, + bool phase_switch, bool phase_flip, bool query /* = false */) const { // ploidy is the count of genotype alleles, so it is also how many haplotypes to report on ploidy_t ploidy = this->ploidies[vi]; int haps = int(idx(ploidy)); - const std::string gt = display_gt(this->orig_gts[vi], ploidy); + + sample_fields fields; + fields.gt = genotype_alleles(this->orig_gts[vi], ploidy); + // QQ was printed with %d before this file wrote records through htslib, so it stays truncated + fields.qq = float(int(this->var_quals[vi])); + fields.sc = sc_idx; + fields.ps = this->phase_sets[vi]; + fields.pb = phase_block; + fields.bs = query ? (phase_switch ? 1 : 0) : bcf_int32_missing; + fields.vp = this->phases[vi] == PHASE_NONE ? + bcf_int32_missing : int32_t(idx(this->phases[vi])); + fields.fe = query ? (phase_flip ? 1 : 0) : bcf_int32_missing; + fields.ge = ac_strs[this->ac_errtype[vi]]; bool swap = this->matched_gt_is_swapped(vi); - std::string errtypes, credits, ref_eds, query_eds, match_types, sync_groups; + float missing_credit; + bcf_float_set_missing(missing_credit); + std::string errtypes, match_types; for (int hap_idx = 0; hap_idx < haps; hap_idx++) { const std::string sep = hap_idx ? "," : ""; hap_t hi = hap_t(hap_idx); // this haplotype carries the reference allele, so it was never evaluated if (!this->var_on_hap(vi, hi)) { - errtypes += sep + "."; credits += sep + "."; ref_eds += sep + "."; - query_eds += sep + "."; match_types += sep + "."; sync_groups += sep + "."; + errtypes += sep + "."; match_types += sep + "."; + fields.bc.push_back(missing_credit); + fields.rd.push_back(bcf_int32_missing); + fields.qd.push_back(bcf_int32_missing); + fields.sg.push_back(bcf_int32_missing); continue; } @@ -445,23 +551,111 @@ void ctgVariants::print_var_sample(FILE* out_fp, int vi, int sc_idx, int phase_b errtypes += sep + (query ? "FP" : "FN"); match_types += sep + "lm"; } - credits += sep + std::to_string(this->credit[hi_matched][vi]); - ref_eds += sep + (this->ref_ed[hi_matched][vi] == 0 ? "." : - std::to_string(this->ref_ed[hi_matched][vi])); - query_eds += sep + (this->ref_ed[hi_matched][vi] == 0 ? "." : - std::to_string(this->query_ed[hi_matched][vi])); - sync_groups += sep + std::to_string(int(this->sync_group[hi_matched][vi])); + fields.bc.push_back(this->credit[hi_matched][vi]); + fields.rd.push_back(this->ref_ed[hi_matched][vi] == 0 ? + bcf_int32_missing : this->ref_ed[hi_matched][vi]); + fields.qd.push_back(this->ref_ed[hi_matched][vi] == 0 ? + bcf_int32_missing : this->query_ed[hi_matched][vi]); + fields.sg.push_back(this->sync_group[hi_matched][vi]); } + fields.bd = errtypes; + fields.bk = match_types; + return fields; +} + - fprintf(out_fp, "\t%s:%s:%s:%s:%s:%s:%d:%d:%s:%d:%d:%s:%s:%s:%s%s", gt.data(), errtypes.data(), - credits.data(), ref_eds.data(), query_eds.data(), match_types.data(), - int(this->var_quals[vi]), sc_idx, sync_groups.data(), - this->phase_sets[vi], phase_block, - query ? (phase_switch ? "1" : "0") : "." , - phase_strs[this->phases[vi]].data(), - query ? (phase_flip ? "1" : "0") : "." , - ac_strs[this->ac_errtype[vi]].data(), - query ? "\n" : ""); +/** @brief Sets one value to the missing value of its type. */ +static void set_missing(int32_t & value) { value = bcf_int32_missing; } +static void set_missing(float & value) { bcf_float_set_missing(value); } + +/** @brief Sets one value to the end-of-vector marker of its type. */ +static void set_vector_end(int32_t & value) { value = bcf_int32_vector_end; } +static void set_vector_end(float & value) { bcf_float_set_vector_end(value); } + +/** + * @brief Concatenates two samples' per-allele values into one buffer of a shared value count. + * + * htslib stores the same number of values for every sample, so a haploid sample beside a diploid + * one is padded with the end-of-vector marker, which the writer prints as a shorter list. A sample + * with no values at all still occupies one slot, holding the missing value. + * @param[in] truth The truth sample's values + * @param[in] query The query sample's values + * @return The truth sample's padded values followed by the query sample's + */ +template +static std::vector pad_per_allele(const std::vector & truth, const std::vector & query) { + const std::vector* samples[2] = {&truth, &query}; + size_t n = std::max(std::max(truth.size(), query.size()), size_t(1)); + std::vector values(2 * n); + for (size_t si = 0; si < 2; si++) { + for (size_t i = 0; i < n; i++) { + T & value = values[si*n + i]; + if (i < samples[si]->size()) value = (*samples[si])[i]; + else if (i == 0) set_missing(value); + else set_vector_end(value); + } + } + return values; +} + +/** @brief Sets one integer FORMAT field of a record, holding one value per sample per allele. */ +static void update_format(const bcf_hdr_t* hdr, bcf1_t* rec, const char* key, + const std::vector & values) { + if (bcf_update_format_int32(hdr, rec, key, values.data(), int(values.size())) < 0) + ERROR("Failed to set FORMAT/%s on summary VCF record", key); +} + +/** @brief Sets one float FORMAT field of a record, holding one value per sample per allele. */ +static void update_format(const bcf_hdr_t* hdr, bcf1_t* rec, const char* key, + const std::vector & values) { + if (bcf_update_format_float(hdr, rec, key, values.data(), int(values.size())) < 0) + ERROR("Failed to set FORMAT/%s on summary VCF record", key); +} + +/** @brief Sets one string FORMAT field of a record, holding one string per sample. */ +static void update_format(const bcf_hdr_t* hdr, bcf1_t* rec, const char* key, + const std::string & truth, const std::string & query) { + const char* values[2] = {truth.data(), query.data()}; + if (bcf_update_format_string(hdr, rec, key, values, 2) < 0) + ERROR("Failed to set FORMAT/%s on summary VCF record", key); +} + + +/** + * @brief Sets every FORMAT field of one record from the two samples' values. + * + * Fields are added in FORMAT declaration order, which is the order htslib writes them in. + * @param[in] hdr Summary VCF header, which must declare every field + * @param[in,out] rec Record whose fixed fields are already set + * @param[in] truth The truth sample's values + * @param[in] query The query sample's values + * @throws ERROR htslib rejects any field's values + */ +void set_record_samples(const bcf_hdr_t* hdr, bcf1_t* rec, + const sample_fields & truth, const sample_fields & query) { + + // a sample with no call reports one missing allele, the '.' genotype + std::vector truth_gt = truth.gt, query_gt = query.gt; + if (truth_gt.empty()) truth_gt.push_back(bcf_gt_missing); + if (query_gt.empty()) query_gt.push_back(bcf_gt_missing); + const std::vector gt = pad_per_allele(truth_gt, query_gt); + if (bcf_update_genotypes(hdr, rec, gt.data(), int(gt.size())) < 0) + ERROR("Failed to set FORMAT/GT on summary VCF record"); + + update_format(hdr, rec, "BD", truth.bd, query.bd); + update_format(hdr, rec, "BC", pad_per_allele(truth.bc, query.bc)); + update_format(hdr, rec, "RD", pad_per_allele(truth.rd, query.rd)); + update_format(hdr, rec, "QD", pad_per_allele(truth.qd, query.qd)); + update_format(hdr, rec, "BK", truth.bk, query.bk); + update_format(hdr, rec, "QQ", std::vector{truth.qq, query.qq}); + update_format(hdr, rec, "SC", std::vector{truth.sc, query.sc}); + update_format(hdr, rec, "SG", pad_per_allele(truth.sg, query.sg)); + update_format(hdr, rec, "PS", std::vector{truth.ps, query.ps}); + update_format(hdr, rec, "PB", std::vector{truth.pb, query.pb}); + update_format(hdr, rec, "BS", std::vector{truth.bs, query.bs}); + update_format(hdr, rec, "VP", std::vector{truth.vp, query.vp}); + update_format(hdr, rec, "FE", std::vector{truth.fe, query.fe}); + update_format(hdr, rec, "GE", truth.ge, query.ge); } /**************************************************************************************************/ diff --git a/src/variant.h b/src/variant.h index cbb72bc..547f7a0 100644 --- a/src/variant.h +++ b/src/variant.h @@ -61,6 +61,33 @@ struct var_fields { EnumArray hap = {}; ///< per-haplotype results, indexed by HAP1 and HAP2 }; +/** + * @struct sample_fields + * @brief One sample's FORMAT values for a summary VCF record, in FORMAT declaration order. + * + * Values are stored htslib-encoded: a per-allele entry that was never evaluated holds + * bcf_int32_missing or a bcf_float_set_missing() float, and the two Number=. String fields hold the + * comma-joined list htslib stores as a single string. A sample that made no call at this locus + * leaves the per-allele vectors empty, which writes '.' for every one of its fields. + */ +struct sample_fields { + std::vector gt; ///< GT alleles, bcf_gt_phased()-encoded, in GT allele order + std::string bd = "."; ///< BD, per-allele benchmark decision (TP/FP/FN) + std::vector bc; ///< BC, per-allele benchmark credit + std::vector rd; ///< RD, per-allele reference edit distance + std::vector qd; ///< QD, per-allele query edit distance + std::string bk = "."; ///< BK, per-allele benchmark category ('gm', 'lm', or '.') + float qq = 0; ///< QQ, variant quality + int32_t sc = 0; ///< SC, supercluster index in contig + std::vector sg; ///< SG, per-allele sync group + int32_t ps = 0; ///< PS, input phase set + int32_t pb = 0; ///< PB, output phase block + int32_t bs = 0; ///< BS, block phase + int32_t vp = 0; ///< VP, variant phase + int32_t fe = 0; ///< FE, flip error + std::string ge = "."; ///< GE, allele count (genotype) error +}; + /** * @class ctgVariants * @brief Store all variant information for a single contig and callset. @@ -77,16 +104,13 @@ class ctgVariants { /** @brief Returns every field of one variant, for copying it into another container. */ var_fields get_var(int idx) const; - /** @brief Writes fixed VCF fields (CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT) for one variant. */ - void print_var_info(FILE* out_fp, std::shared_ptr ref, - const std::string & ctg, int idx); + /** @brief Sets the fixed VCF fields (CHROM, POS, ID, REF, ALT, QUAL, FILTER) of one record. */ + void set_var_record(const bcf_hdr_t* hdr, bcf1_t* rec, std::shared_ptr ref, + const std::string & ctg, int idx) const; - /** @brief Writes dot-separated empty sample fields for a variant with no call on this haplotype. */ - void print_var_empty(FILE* out_fp, int sc_idx, int phase_block, bool query = false); - - /** @brief Writes sample-specific FORMAT fields for one variant to output VCF. */ - void print_var_sample(FILE* out_fp, int vi, int sc_idx, int phase_block, - bool phase_switch, bool phase_flip, bool query = false); + /** @brief Returns one sample's FORMAT values for a variant it called. */ + sample_fields var_sample_fields(int vi, int sc_idx, int phase_block, + bool phase_switch, bool phase_flip, bool query = false) const; /** @brief Returns true if a variant is present on the specified haplotype. */ bool var_on_hap(int var_idx, hap_t hap, bool matched = false) const; @@ -170,6 +194,17 @@ class variantData { /** @brief Classifies a record's raw GT array into its parse-time genotype shape. */ gtparse_t classify_gt(const int32_t * gt, int ngt); +/** @brief Builds the summary VCF header, declaring every FORMAT field and the TRUTH/QUERY samples. */ +bcf_hdr_t* summary_vcf_header(const std::vector & contigs, + const std::vector & lengths); + +/** @brief Returns the FORMAT values of a sample that made no call at a locus. */ +sample_fields empty_sample_fields(int sc_idx, int phase_block); + +/** @brief Sets every FORMAT field of one record from the two samples' values. */ +void set_record_samples(const bcf_hdr_t* hdr, bcf1_t* rec, + const sample_fields & truth, const sample_fields & query); + /** @brief Parses variants from a VCF file into a variantData container, with filtering and validation. */ void parse_variants(const std::string & vcf_fn, std::shared_ptr variant_data, diff --git a/tests/integration/test-integration.yml b/tests/integration/test-integration.yml index 77026a9..94f6e29 100644 --- a/tests/integration/test-integration.yml +++ b/tests/integration/test-integration.yml @@ -29,8 +29,8 @@ - "sc1\t140\t.\tT\tC" - "sc1\t146\t.\tG\tA" - "sc1\t256\t.\tG\tA" - - ":TP,TP:1.000000,1.000000:" - - ":FN,FN:0.000000,0.000000:" + - ":TP,TP:1,1:" + - ":FN,FN:0,0:" must_not_contain: - ":FP" @@ -52,7 +52,7 @@ - "INDEL\tNONE\t0\t2\t2\t0\t0" - path: partial_match_del_ct07_summary.vcf contains: - - ":TP,TP:0.800000,0.800000:" + - ":TP,TP:0.8,0.8:" must_not_contain: - ":FN" - ":FP" @@ -74,8 +74,8 @@ - "INDEL\tNONE\t0\t0\t0\t2\t2" - path: partial_match_del_ct1_summary.vcf contains: - - ":FN,FN:0.000000,0.000000:" - - ":FP,FP:0.000000,0.000000:" + - ":FN,FN:0,0:" + - ":FP,FP:0,0:" must_not_contain: - ":TP" @@ -97,8 +97,8 @@ contains: - "sc1\t200\t.\tA\tG" - "sc1\t200\t.\tA\tC" - - ":FN,FN:0.000000,0.000000:" - - ":FP,FP:0.000000,0.000000:" + - ":FN,FN:0,0:" + - ":FP,FP:0,0:" must_not_contain: - ":TP" @@ -123,7 +123,7 @@ - "INDEL\tNONE\t0\t2\t2\t0\t0" - path: entangled_overlap_del_out_summary.vcf contains: - - ":TP,TP:0.800000,0.800000:" + - ":TP,TP:0.8,0.8:" - name: test_single-pass_adjacent-fn-no-overcredit # adjacent_fn_no_overcredit: regression guard for the consecutive-bypass excision fix. @@ -145,10 +145,10 @@ - path: adjacent_fn_no_overcredit_out_summary.vcf contains: - "sc1\t300\t.\tA\tC" - - ":TP,TP:1.000000,1.000000:1,1:0,0:" - - ":FN,FN:0.000000,0.000000:" + - ":TP,TP:1,1:1,1:0,0:" + - ":FN,FN:0,0:" must_not_contain: - - ":TP,TP:1.000000,1.000000:2," + - ":TP,TP:1,1:2," - ":FP" - name: test_summary-vcf_one-record-per-variant @@ -172,13 +172,13 @@ - path: record_shapes_out_summary.vcf contains: # hom SNP, hom CPX (as INS then DEL), and hom DEL: one record each, two values per field - - "sc1\t200\t.\tA\tG\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:TP,TP:1.000000,1.000000:1,1:0,0:gm,gm:50:0:3,3:0:0:.:.:.:.\t1|1:TP,TP:1.000000,1.000000:1,1:0,0:gm,gm:50:0:3,3:0:0:0:.:0:." + - "sc1\t200\t.\tA\tG\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:TP,TP:1,1:1,1:0,0:gm,gm:50:0:3,3:0:0:.:.:.:.\t1|1:TP,TP:1,1:1,1:0,0:gm,gm:50:0:3,3:0:0:0:.:0:." - "sc1\t209\t.\tG\tGTT\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:TP,TP:" - "sc1\t209\t.\tGCAAGA\tG\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:TP,TP:" - "sc1\t220\t.\tCAACT\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:TP,TP:" # het-alt: two co-located records, each carrying '.' for its reference allele - - "sc1\t250\t.\tA\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|0:TP,.:1.000000,.:1,.:0,.:gm,.:50:1:0,.:0:0:.:.:.:.\t1|0:TP,.:" - - "sc1\t250\t.\tA\tG\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t0|1:.,TP:.,1.000000:.,1:.,0:.,gm:50:1:.,0:0:0:.:.:.:.\t0|1:.,TP:" + - "sc1\t250\t.\tA\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|0:TP,.:1,.:1,.:0,.:gm,.:50:1:0,.:0:0:.:.:.:.\t1|0:TP,.:" + - "sc1\t250\t.\tA\tG\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t0|1:.,TP:.,1:.,1:.,0:.,gm:50:1:.,0:0:0:.:.:.:.\t0|1:.,TP:" must_not_contain: # the pre-#103 shape: a homozygous call written once per haplotype - "\t1|0:TP:" @@ -212,7 +212,7 @@ # truth-only direction below, where it used to be dropped - path: one_sided_contig_query_only_summary.vcf contains: - - "sc2\t50\t.\tT\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t.:.:.:.:.:.:.:0:.:.:0:.:.:.:.\t1|1:FP,FP:0.000000,0.000000:.,.:.,.:.,.:50:0:0,0:0:0:0:.:0:." + - "sc2\t50\t.\tT\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t.:.:.:.:.:.:.:0:.:.:0:.:.:.:.\t1|1:FP,FP:0,0:.,.:.,.:.,.:50:0:0,0:0:0:0:.:0:." - name: test_one-sided-contig_truth-only # The mirror of the case above: sc2 is called by the truth only. Superclustering used to @@ -238,7 +238,7 @@ # sample reports its own 1|1 call, one FN per allele. - path: one_sided_contig_truth_only_summary.vcf contains: - - "sc2\t50\t.\tT\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:FN,FN:0.000000,0.000000:.,.:.,.:.,.:50:0:0,0:0:0:.:.:.:.\t.:.:.:.:.:.:.:0:.:.:0:.:.:.:." + - "sc2\t50\t.\tT\tC\t.\tPASS\t.\tGT:BD:BC:RD:QD:BK:QQ:SC:SG:PS:PB:BS:VP:FE:GE\t1|1:FN,FN:0,0:.,.:.,.:.,.:50:0:0,0:0:0:.:.:.:.\t.:.:.:.:.:.:.:0:.:.:0:.:.:.:." # --- variants at reference position 0 (synthetic sc1 reference) --- diff --git a/tests/unit/src/test_phase.cpp b/tests/unit/src/test_phase.cpp index 2533a0d..98ebcd4 100644 --- a/tests/unit/src/test_phase.cpp +++ b/tests/unit/src/test_phase.cpp @@ -10,6 +10,8 @@ #include "gtest/gtest.h" +#include "htslib/vcf.h" + #include "../../../src/defs.h" #include "../../../src/globals.h" #include "../../../src/phase.h" @@ -1358,7 +1360,7 @@ TEST(WriteSummaryVcf, HomSnpIsOneRecordWithTwoValues) { std::vector sample = sole_query_sample(shape_vcf(dir, qvars), SPACING + 1); EXPECT_EQ("1|1", sample.at(FMT_GT)); EXPECT_EQ("TP,TP", sample.at(FMT_BD)); - EXPECT_EQ("1.000000,1.000000", sample.at(FMT_BC)); + EXPECT_EQ("1,1", sample.at(FMT_BC)); EXPECT_EQ("3,3", sample.at(FMT_RD)); EXPECT_EQ("0,0", sample.at(FMT_QD)); EXPECT_EQ("gm,gm", sample.at(FMT_BK)); @@ -1384,7 +1386,7 @@ TEST(WriteSummaryVcf, HomIndelIsOneRecordWithTwoValues) { std::vector sample = split(cols.at(QUERY_COL), ':'); EXPECT_EQ("1|1", sample.at(FMT_GT)); EXPECT_EQ("TP,TP", sample.at(FMT_BD)); - EXPECT_EQ("1.000000,1.000000", sample.at(FMT_BC)); + EXPECT_EQ("1,1", sample.at(FMT_BC)); } // A heterozygous call still carries one value per GT allele, but its reference allele was never @@ -1399,7 +1401,7 @@ TEST(WriteSummaryVcf, HetRecordDotsTheReferenceAllele) { std::vector sample = sole_query_sample(shape_vcf(dir, qvars), SPACING + 1); EXPECT_EQ("1|0", sample.at(FMT_GT)); EXPECT_EQ("TP,.", sample.at(FMT_BD)); - EXPECT_EQ("1.000000,.", sample.at(FMT_BC)); + EXPECT_EQ("1,.", sample.at(FMT_BC)); EXPECT_EQ("5,.", sample.at(FMT_RD)); EXPECT_EQ("2,.", sample.at(FMT_SG)); } @@ -1415,7 +1417,7 @@ TEST(WriteSummaryVcf, HaploidRecordCarriesOneValue) { std::vector sample = sole_query_sample(shape_vcf(dir, qvars), SPACING + 1); EXPECT_EQ("1", sample.at(FMT_GT)); EXPECT_EQ("TP", sample.at(FMT_BD)); - EXPECT_EQ("1.000000", sample.at(FMT_BC)); + EXPECT_EQ("1", sample.at(FMT_BC)); EXPECT_EQ("5", sample.at(FMT_RD)); EXPECT_EQ("2", sample.at(FMT_SG)); } @@ -1434,7 +1436,7 @@ TEST(WriteSummaryVcf, PerAlleleValuesFollowTheGenotypeSwap) { std::vector sample = sole_query_sample(shape_vcf(dir, qvars), SPACING + 1); EXPECT_EQ("1|0", sample.at(FMT_GT)); EXPECT_EQ("TP,.", sample.at(FMT_BD)); - EXPECT_EQ("1.000000,.", sample.at(FMT_BC)); + EXPECT_EQ("1,.", sample.at(FMT_BC)); EXPECT_EQ("2,.", sample.at(FMT_SG)); } @@ -1451,7 +1453,7 @@ TEST(WriteSummaryVcf, PerAlleleValuesAreInGenotypeAlleleOrder) { EXPECT_EQ("3,7", sample.at(FMT_RD)); EXPECT_EQ("0,2", sample.at(FMT_QD)); EXPECT_EQ("1,6", sample.at(FMT_SG)); - EXPECT_EQ("1.000000,0.750000", sample.at(FMT_BC)); + EXPECT_EQ("1,0.75", sample.at(FMT_BC)); } // A het-alt (1|2) source record is parsed into two entries with different ALTs, and the @@ -1487,6 +1489,43 @@ TEST(WriteSummaryVcf, HetAltStaysTwoColocatedRecords) { EXPECT_EQ("0|1", split(split(recs[1], '\t').at(QUERY_COL), ':').at(FMT_GT)); } +// The records are built as bcf1_t and written with bcf_write(), so htslib must be able to read +// back what it wrote: a field whose values disagree with its header declaration is rejected on the +// way in, which no assertion over the rendered text would catch. The two variants differ in ploidy +// so that the per-allele lists are padded to a shared length, the encoding most likely to break. +TEST(WriteSummaryVcf, RecordsReadBackThroughHtslib) { + GlobalsGuard guard; + TempDir dir; + pipeline_result result = run_pipeline(dir, make_ploidy_qvars({PLOIDY_DIPLOID, PLOIDY_HAPLOID})); + const std::string vcf_fn = dir.path("readback.vcf"); + { + StderrToFile redirect(dir.path("readback.log")); + result.data->write_summary_vcf(vcf_fn); + } + + htsFile* vcf = bcf_open(vcf_fn.data(), "r"); + ASSERT_NE(nullptr, vcf); + bcf_hdr_t* hdr = bcf_hdr_read(vcf); + ASSERT_NE(nullptr, hdr); + EXPECT_EQ(2, bcf_hdr_nsamples(hdr)); + EXPECT_STREQ("TRUTH", hdr->samples[0]); + EXPECT_STREQ("QUERY", hdr->samples[1]); + + bcf1_t* rec = bcf_init(); + int records = 0; + int read_ret = 0; + while ((read_ret = bcf_read(vcf, hdr, rec)) == 0) { + EXPECT_EQ(0, rec->errcode); + records++; + } + EXPECT_EQ(-1, read_ret); // end of file rather than a parse failure + EXPECT_EQ(2, records); + + bcf_destroy(rec); + bcf_hdr_destroy(hdr); + bcf_close(vcf); +} + // Number=P would declare the one-value-per-GT-allele cardinality these fields carry, but it is a // VCF 4.4 addition htslib only supports from 1.23, so Number=. is declared and the count and // order are stated in the description instead. From 1f1165ba9d9902e09f9d4645930a10099b2eb2f6 Mon Sep 17 00:00:00 2001 From: Tim Dunn Date: Sat, 8 Aug 2026 00:29:51 -0400 Subject: [PATCH 2/2] fix(tests): reach htslib through variant.h in test_phase.cpp The build copies htslib's headers into src/, so a bare "htslib/vcf.h" only resolves for a file in that directory; from tests/unit/src it built locally only because an unrelated Homebrew prefix happens to sit on the default include path. Include src/variant.h instead, as test_variant.cpp does. --- tests/unit/src/test_phase.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/unit/src/test_phase.cpp b/tests/unit/src/test_phase.cpp index 98ebcd4..fe1b8de 100644 --- a/tests/unit/src/test_phase.cpp +++ b/tests/unit/src/test_phase.cpp @@ -10,11 +10,12 @@ #include "gtest/gtest.h" -#include "htslib/vcf.h" - #include "../../../src/defs.h" #include "../../../src/globals.h" #include "../../../src/phase.h" +// htslib's headers are reached through variant.h: the build copies them into src/, so only a file +// in that directory resolves a bare "htslib/vcf.h" +#include "../../../src/variant.h" #include "test_helpers.h" namespace {