diff --git a/README.md b/README.md index 3d3d940..fee94a8 100644 --- a/README.md +++ b/README.md @@ -143,6 +143,8 @@ To include more details on intermediate results, run it again at higher verbosit Both input VCFs must be coordinate-sorted: each contig's records must be grouped together, and within a contig positions must not decrease. vcfdist exits with an error naming the offending record otherwise. +Every contig either input VCF carries must also be present in the reference FASTA, whether or not a BED file restricts which regions are evaluated. vcfdist exits with an error naming the contig and the callset otherwise. A query VCF carrying decoy or alt contigs (`chrEBV`, `HLA-*`) therefore needs a reference that includes them, or those records removed. + ## Wiki The [vcfdist wiki](https://github.com/TimD1/vcfdist/wiki) has helpful information on [command-line parameters](https://github.com/TimD1/vcfdist/wiki/02-Parameters-and-Usage), [output documentation](https://github.com/TimD1/vcfdist/wiki/09-Outputs), and [implementation](https://github.com/TimD1/vcfdist/wiki/01-Overview). diff --git a/src/bed.cpp b/src/bed.cpp index 1368e14..fb7a72d 100644 --- a/src/bed.cpp +++ b/src/bed.cpp @@ -481,14 +481,49 @@ static std::string ploidy_set_str(const std::set & ploidies) { } /** - * @brief Intersects reference FASTA, query VCF, truth VCF, and optional BED regions, retaining only common contigs. + * @brief Gives one callset an empty entry for every contig the other carries but it lacks. * - * @param[in] query_ptr A pointer to the query variantData. - * @param[in] truth_ptr A pointer to the truth variantData. + * Consumers index variants[hap][ctg] over the union of the two callsets' contig lists and + * dereference the result without a null check -- superclusterData's merge does, since a contig is + * superclustered whenever either callset has variants there. std::unordered_map::operator[] would + * insert a null shared_ptr for an absent contig, so the two lists are kept mutually inclusive here + * rather than guarded at each consumer. + * @param[in] ref_ptr A pointer to the reference fastaData, supplying the contig's length. + * @param[in] from_ptr The callset whose contigs are being copied across. + * @param[in,out] to_ptr The callset gaining an empty entry per contig it lacks. + */ +static void pair_missing_contigs(const std::shared_ptr & ref_ptr, + const std::shared_ptr & from_ptr, + const std::shared_ptr & to_ptr) { + for (const std::string & ctg : from_ptr->contigs) { + if (std::find(to_ptr->contigs.begin(), to_ptr->contigs.end(), ctg) + != to_ptr->contigs.end()) continue; + to_ptr->variants[HAP1][ctg] = std::shared_ptr(new ctgVariants(ctg)); + to_ptr->variants[HAP2][ctg] = std::shared_ptr(new ctgVariants(ctg)); + to_ptr->contigs.push_back(ctg); + // safe: intersect_contigs() has already rejected any VCF contig the FASTA lacks + to_ptr->lengths.push_back(ref_ptr->lengths.at(ctg)); + to_ptr->observed_ploidies.push_back({}); + } +} + +/** + * @brief Reconciles the contigs of the reference FASTA, query VCF, truth VCF, and optional BED. + * + * Every contig either VCF carries is retained, along with its reference sequence, so that variants + * on it remain reportable; contigs outside the BED are simply never evaluated, since + * parse_variants() has already discarded their variants as BED_OFFCTG. Contigs the BED names but a + * VCF lacks are injected empty, so that both callsets cover every evaluated region, and the two + * callsets are then paired so that each is indexable by every contig the other carries. + * + * @param[in,out] query_ptr A pointer to the query variantData, gaining any missing BED contigs. + * @param[in,out] truth_ptr A pointer to the truth variantData, gaining any missing BED contigs. * @param[in] ref_ptr A pointer to the reference fastaData. - * @throws WARNING if contigs in either VCF are not present in either the other VCF or BED file. - * @throws WARNING if corresponding contigs in the truth and query VCFs observed differing ploidies. - * @throws ERROR if a contig to be evaluated is not present in the reference FASTA. + * @throws WARNING if a BED contig is present in one VCF but not the other. + * @throws WARNING if corresponding evaluated contigs in the truth and query VCFs observed differing + * ploidies. + * @throws ERROR if a contig in either VCF is not present in the reference FASTA. + * @throws ERROR if a contig in the BED file is not present in the reference FASTA. */ void intersect_contigs( std::shared_ptr query_ptr, @@ -497,61 +532,32 @@ void intersect_contigs( if (g.verbosity >= 1) INFO(" "); if (g.verbosity >= 1) INFO(" Checking contigs:"); - if (g.bed_exists) { // use BED to determine contigs - - // remove all extraneous contigs in query VCF not in BED - std::vector::iterator itr = query_ptr->contigs.begin(); - while (itr != query_ptr->contigs.end()) { // query - if (std::find(g.bed.contigs.begin(), g.bed.contigs.end(), - *itr) == g.bed.contigs.end()) { - query_ptr->lengths.erase(query_ptr->lengths.begin() + - (itr - query_ptr->contigs.begin())); - query_ptr->observed_ploidies.erase(query_ptr->observed_ploidies.begin() + - (itr - query_ptr->contigs.begin())); - query_ptr->variants[HAP1].erase(*itr); - query_ptr->variants[HAP2].erase(*itr); - std::string dropped_ctg = *itr; // save name, erase() invalidates itr - itr = query_ptr->contigs.erase(itr); - if (g.verbosity >= 2) - WARN("Ignoring %s from QUERY VCF, not in BED file.", dropped_ctg.data()); - } else ++itr; - } - // remove all extraneous contigs in truth VCF not in BED - itr = truth_ptr->contigs.begin(); - while (itr != truth_ptr->contigs.end()) { // truth - if (std::find(g.bed.contigs.begin(), g.bed.contigs.end(), - *itr) == g.bed.contigs.end()) { - truth_ptr->lengths.erase(truth_ptr->lengths.begin() + - (itr - truth_ptr->contigs.begin())); - truth_ptr->observed_ploidies.erase(truth_ptr->observed_ploidies.begin() + - (itr - truth_ptr->contigs.begin())); - truth_ptr->variants[HAP1].erase(*itr); - truth_ptr->variants[HAP2].erase(*itr); - std::string dropped_ctg = *itr; // save name, erase() invalidates itr - itr = truth_ptr->contigs.erase(itr); - if (g.verbosity >= 2) - WARN("Ignoring %s from TRUTH VCF, not in BED file.", dropped_ctg.data()); - } else ++itr; - } - // remove all extraneous contigs in ref FASTA not in BED - auto itr2 = ref_ptr->fasta.begin(); - while (itr2 != ref_ptr->fasta.end()) { // fasta - if (std::find(g.bed.contigs.begin(), g.bed.contigs.end(), - itr2->first) == g.bed.contigs.end()) { - itr2 = ref_ptr->fasta.erase(itr2); - } else itr2++; + // Every contig in either callset is retained, so every one of them needs reference sequence: + // set_var_record() anchors an INS/DEL in ref->fasta when writing a variant back out. + const std::vector< std::pair > > callsets = + {{QUERY, query_ptr}, {TRUTH, truth_ptr}}; + for (const auto & [callset, vcf_ptr] : callsets) { + for (const std::string & ctg : vcf_ptr->contigs) { + if (ref_ptr->fasta.find(ctg) == ref_ptr->fasta.end()) + ERROR("Contig '%s' found in %s VCF but not reference FASTA.", + ctg.data(), callset_strs[callset].data()); } + } - // warn if list of truth and query contigs are not the same - for (std::string ctg : query_ptr->contigs) { - if (std::find(truth_ptr->contigs.begin(), - truth_ptr->contigs.end(), ctg) == truth_ptr->contigs.end()) + if (g.bed_exists) { // use BED to determine contigs + + // Warn if the truth and query contig lists differ, over the BED contigs alone: nothing is + // evaluated outside them, so a one-sided contig there yields neither FPs nor FNs to warn + // about. Runs before the injection below, which would otherwise mask every difference. + for (const std::string & ctg : g.bed.contigs) { + bool in_query = std::find(query_ptr->contigs.begin(), + query_ptr->contigs.end(), ctg) != query_ptr->contigs.end(); + bool in_truth = std::find(truth_ptr->contigs.begin(), + truth_ptr->contigs.end(), ctg) != truth_ptr->contigs.end(); + if (in_query && !in_truth) WARN("Contig '%s' found in query VCF but not truth VCF." " All query variants on '%s' will be false positives.", ctg.data(), ctg.data()); - } - for (std::string ctg : truth_ptr->contigs) { - if (std::find(query_ptr->contigs.begin(), - query_ptr->contigs.end(), ctg) == query_ptr->contigs.end()) + if (in_truth && !in_query) WARN("Contig '%s' found in truth VCF but not query VCF." " All truth variants on '%s' will be false negatives.", ctg.data(), ctg.data()); } @@ -586,59 +592,43 @@ void intersect_contigs( } else { // use truth VCF to determine contigs - // ensure fasta contains all contigs - for (std::string ctg : truth_ptr->contigs) { - if (ref_ptr->fasta.find(ctg) == ref_ptr->fasta.end()) - ERROR("Contig '%s' found in truth VCF but not reference FASTA. Please provide BED file.", ctg.data()); - } - - // ensure query/truth VCFs contain the same contigs (even if devoid of variants) - for (int i = 0; i < int(query_ptr->contigs.size()); i++) { - std::string ctg = query_ptr->contigs[i]; - if (std::find(truth_ptr->contigs.begin(), - truth_ptr->contigs.end(), ctg) == truth_ptr->contigs.end()) { + // warn if the truth and query contig lists differ; with no BED every contig is evaluated, + // so a one-sided contig here really does yield only false positives or only false negatives + for (const std::string & ctg : query_ptr->contigs) { + if (std::find(truth_ptr->contigs.begin(), + truth_ptr->contigs.end(), ctg) == truth_ptr->contigs.end()) WARN("Contig '%s' found in query VCF but not truth VCF." " All query variants on '%s' will be false positives.", ctg.data(), ctg.data()); - truth_ptr->variants[HAP1][ctg] = - std::shared_ptr(new ctgVariants(ctg)); - truth_ptr->variants[HAP2][ctg] = - std::shared_ptr(new ctgVariants(ctg)); - truth_ptr->contigs.push_back(ctg); - truth_ptr->lengths.push_back(ref_ptr->lengths.at(ctg)); - truth_ptr->observed_ploidies.push_back({}); - } } - for (int i = 0; i < int(truth_ptr->contigs.size()); i++) { - std::string ctg = truth_ptr->contigs[i]; - if (std::find(query_ptr->contigs.begin(), - query_ptr->contigs.end(), ctg) == query_ptr->contigs.end()) { + for (const std::string & ctg : truth_ptr->contigs) { + if (std::find(query_ptr->contigs.begin(), + query_ptr->contigs.end(), ctg) == query_ptr->contigs.end()) WARN("Contig '%s' found in truth VCF but not query VCF." " All truth variants on '%s' will be false negatives.", ctg.data(), ctg.data()); - query_ptr->variants[HAP1][ctg] = - std::shared_ptr(new ctgVariants(ctg)); - query_ptr->variants[HAP2][ctg] = - std::shared_ptr(new ctgVariants(ctg)); - query_ptr->contigs.push_back(ctg); - query_ptr->lengths.push_back(ref_ptr->lengths.at(ctg)); - query_ptr->observed_ploidies.push_back({}); - } - } - - // remove extra contigs from ref FASTA - auto itr = ref_ptr->fasta.begin(); - while (itr != ref_ptr->fasta.end()) { - if (std::find(truth_ptr->contigs.begin(), truth_ptr->contigs.end(), - itr->first) == truth_ptr->contigs.end()) { - itr = ref_ptr->fasta.erase(itr); - } else itr++; } } - // verify the observed ploidies match for all truth/query contigs + // Pair the two callsets' contig lists, whether or not a BED narrowed what is evaluated. The BED + // branch above only injects the contigs the BED names, so a one-sided contig outside the BED + // would otherwise survive with no counterpart for the merge to index. + pair_missing_contigs(ref_ptr, query_ptr, truth_ptr); + pair_missing_contigs(ref_ptr, truth_ptr, query_ptr); + + // verify the observed ploidies match for all evaluated truth/query contigs for (int i = 0; i < int(truth_ptr->contigs.size()); i++) { std::string ctg = truth_ptr->contigs[i]; - int query_ctg_idx = std::find(query_ptr->contigs.begin(), - query_ptr->contigs.end(), ctg) - query_ptr->contigs.begin(); + + // Ploidy is recorded before parse_variants() applies the BED filter, so a contig outside + // the BED carries ploidies no evaluation consults; disagreement there is not actionable. + if (g.bed_exists && std::find(g.bed.contigs.begin(), g.bed.contigs.end(), ctg) + == g.bed.contigs.end()) continue; + + // Unreachable: pair_missing_contigs() has given the query every contig the truth carries. + // Guarded anyway, since indexing on a failed find() would read past the end of the vector + // rather than report anything. + auto query_itr = std::find(query_ptr->contigs.begin(), query_ptr->contigs.end(), ctg); + if (query_itr == query_ptr->contigs.end()) continue; + int query_ctg_idx = query_itr - query_ptr->contigs.begin(); int truth_ctg_idx = i; const std::set & truth_ploidies = truth_ptr->observed_ploidies[truth_ctg_idx]; const std::set & query_ploidies = query_ptr->observed_ploidies[query_ctg_idx]; diff --git a/tests/integration/data/README.md b/tests/integration/data/README.md index 51096b8..6712099 100644 --- a/tests/integration/data/README.md +++ b/tests/integration/data/README.md @@ -111,6 +111,27 @@ skipped any contig the query does not call on, and the query-only direction pins no query call on `sc2` there is no phase block or phase to report there, so its truth record reports `PB=0` and `BS=.`. +### bed_absent_contig — a one-sided contig outside the BED (#246) + +Reuses the `one_sided_contig` VCF pair against `synthetic.bed` (`sc1` alone) instead of +`synthetic_2ctg.bed`, in both directions, so `sc2` is called by one callset and lies outside the +evaluated regions. `sc2` is retained rather than pruned, but `parse_variants()` has already +discarded its calls as off-contig, so both directions score SNP 2/2/0/0 — the same as running with +`sc2` removed from the input entirely. Each direction pins that `sc2` reaches `*summary.vcf`'s +header and carries no record there, which is what a later change writes retained calls into. + +### undeclared_contig — a contig one callset's header omits (#246) + +- `undeclared_contig_query.vcf` — declares `sc1` and `sc2`; SNP `sc1` 200 A>G and SNP `sc2` 50 T>C. +- `undeclared_contig_truth.vcf` — declares `sc1` only; SNP `sc1` 200 A>G. + +Run against `synthetic.bed`, so `sc2` is also outside the evaluated regions. Unlike the +`one_sided_contig` pair, whose files both declare `sc2`, the truth here never declares it at all — +the distinction `parse_variants()` draws, since it creates a `ctgVariants` per *header* contig but +appends to `contigs` per *record*. Superclustering merges `variants[hap][ctg]` across the union of +both contig lists and dereferences the result without a null check, so a contig left unpaired by +`intersect_contigs()` segfaults instead of being reported. Metrics match the declared-by-both case. + ### unsorted_position_query — a coordinate-unsorted VCF (#232) - `unsorted_position_query.vcf` — the `swallowed_snps` query's three SNPs with the last two swapped, diff --git a/tests/integration/data/undeclared_contig_query.vcf b/tests/integration/data/undeclared_contig_query.vcf new file mode 100644 index 0000000..a116f91 --- /dev/null +++ b/tests/integration/data/undeclared_contig_query.vcf @@ -0,0 +1,8 @@ +##fileformat=VCFv4.2 +##FILTER= +##FORMAT= +##contig= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE +sc1 200 . A G 50 PASS . GT 1/1 +sc2 50 . T C 50 PASS . GT 1/1 diff --git a/tests/integration/data/undeclared_contig_truth.vcf b/tests/integration/data/undeclared_contig_truth.vcf new file mode 100644 index 0000000..d5a1ae8 --- /dev/null +++ b/tests/integration/data/undeclared_contig_truth.vcf @@ -0,0 +1,6 @@ +##fileformat=VCFv4.2 +##FILTER= +##FORMAT= +##contig= +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE +sc1 200 . A G 50 PASS . GT 1/1 diff --git a/tests/integration/test-integration.yml b/tests/integration/test-integration.yml index ea064e8..e9f29b8 100644 --- a/tests/integration/test-integration.yml +++ b/tests/integration/test-integration.yml @@ -240,6 +240,88 @@ 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,0:.,.:.,.:.,.:50:0:0,0:0:0:.:.:.:.\t.:.:.:.:.:.:.:0:.:.:0:.:.:.:." +# --- contigs outside the BED file (synthetic 2-contig reference) --- +# The same VCF pair as the one-sided-contig scenarios above, run against synthetic.bed (sc1 only) +# instead of synthetic_2ctg.bed, so sc2 falls outside the evaluated regions in each direction. +# The contig is retained rather than pruned, so it reaches the output contig list; parse_variants() +# has already discarded its variants as BED_OFFCTG, so it contributes nothing to any metric. + +- name: test_bed-absent-contig_query-only-not-evaluated + # sc2 is called by the query alone and lies outside the BED. Contrast with + # test_one-sided-contig_query-only above, which runs the identical VCF pair against a BED that + # does cover sc2 and scores those same two calls as false positives. + command: bash -c '${VCFDIST_REPO_PATH}/src/vcfdist + "${VCFDIST_REPO_PATH}/tests/integration/data/one_sided_contig_both.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/one_sided_contig_sc1only.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.fasta" + -b "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.bed" + -p "bed_absent_contig_query_only_"' + files: + # only sc1 is evaluated, so the two sc2 calls are neither true nor false positives + - path: bed_absent_contig_query_only_precision-recall-summary.tsv + contains: + - "SNP\tNONE\t0\t2\t2\t0\t0" + - path: bed_absent_contig_query_only_query.tsv + contains: + - "sc1\t199\t0\tA\tG\t50.00\tSNP\tTP" + must_not_contain: + - "sc2" + # sc2 survives as far as the output contig list, but carries no record + - path: bed_absent_contig_query_only_summary.vcf + contains: + - "##contig=" + must_not_contain: + - "sc2\t50\t." + +- name: test_bed-absent-contig_truth-only-not-evaluated + # The mirror: sc2 is called by the truth alone and lies outside the BED, so nothing injects a + # query counterpart for it. Comparing the two callsets' observed ploidies on sc2 would then index + # the query's ploidy vector out of bounds, which this scenario is the integration-level guard for. + command: bash -c '${VCFDIST_REPO_PATH}/src/vcfdist + "${VCFDIST_REPO_PATH}/tests/integration/data/one_sided_contig_sc1only.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/one_sided_contig_both.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.fasta" + -b "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.bed" + -p "bed_absent_contig_truth_only_"' + stderr: + must_not_contain: + - "contig 'sc2' has ploid" + files: + - path: bed_absent_contig_truth_only_precision-recall-summary.tsv + contains: + - "SNP\tNONE\t0\t2\t2\t0\t0" + - path: bed_absent_contig_truth_only_truth.tsv + contains: + - "sc1\t199\t0\tA\tG\t50.00\tSNP\tTP" + must_not_contain: + - "sc2" + - path: bed_absent_contig_truth_only_summary.vcf + contains: + - "##contig=" + must_not_contain: + - "sc2\t50\t." + +- name: test_bed-absent-contig_undeclared-by-other-callset + # The one_sided_contig VCFs above both declare sc2 in their headers, so both callsets can be + # indexed by it. Here the truth's header omits sc2 entirely, which is what parse_variants() + # distinguishes: it creates a ctgVariants per header contig, and appends to `contigs` per record. + # Superclustering merges variants[hap][ctg] across the union of both contig lists and dereferences + # the result, so an unpaired sc2 segfaults rather than being reported. + command: bash -c '${VCFDIST_REPO_PATH}/src/vcfdist + "${VCFDIST_REPO_PATH}/tests/integration/data/undeclared_contig_query.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/undeclared_contig_truth.vcf" + "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.fasta" + -b "${VCFDIST_REPO_PATH}/tests/integration/data/synthetic.bed" + -p "undeclared_contig_"' + files: + # sc2 is outside the BED, so pairing it changes nothing that is counted + - path: undeclared_contig_precision-recall-summary.tsv + contains: + - "SNP\tNONE\t0\t2\t2\t0\t0" + - path: undeclared_contig_query.tsv + must_not_contain: + - "sc2" + # --- variants at reference position 0 (synthetic sc1 reference) --- - name: test_contig-start-snp-scored-as-tp diff --git a/tests/unit/src/test_bed.cpp b/tests/unit/src/test_bed.cpp index 50403e7..a57754f 100644 --- a/tests/unit/src/test_bed.cpp +++ b/tests/unit/src/test_bed.cpp @@ -1035,7 +1035,9 @@ TEST(CheckStrataContigs, AllZeroOverlapNamesLikelyCause) { /* intersect_contigs ******************************************************************************/ -TEST(IntersectContigs, BedDropsExtraneous) { +// A contig outside the BED survives so that its variants remain reportable. It is already empty: +// parse_variants() discards every variant on it as BED_OFFCTG, so keeping it costs an entry. +TEST(IntersectContigs, BedKeepsContigOutsideBed) { GlobalsGuard guard; g.bed_exists = true; g.bed = make_bed("chr1", {{0, 10}}); @@ -1048,16 +1050,154 @@ TEST(IntersectContigs, BedDropsExtraneous) { intersect_contigs(query, truth, ref); testing::internal::GetCapturedStderr(); - // chr2 is absent from the BED, so it is dropped from the query along with its parallel fields - EXPECT_EQ(std::vector({"chr1"}), query->contigs); - EXPECT_EQ(std::vector({12}), query->lengths); - EXPECT_EQ(std::vector< std::set >({{2}}), query->observed_ploidies); - EXPECT_EQ(size_t(0), query->variants[HAP1].count("chr2")); - EXPECT_EQ(size_t(0), query->variants[HAP2].count("chr2")); + // chr2 keeps its entry in contigs and in every field parallel to it + EXPECT_EQ(std::vector({"chr1", "chr2"}), query->contigs); + EXPECT_EQ(std::vector({12, 12}), query->lengths); + EXPECT_EQ(std::vector< std::set >({{2}, {2}}), query->observed_ploidies); + EXPECT_EQ(size_t(1), query->variants[HAP1].count("chr2")); + EXPECT_EQ(size_t(1), query->variants[HAP2].count("chr2")); + + // set_var_record() anchors an INS/DEL in the reference sequence, so it must survive too + EXPECT_EQ(size_t(2), ref->fasta.size()); + EXPECT_EQ("TTTTTTTTTTTT", ref->fasta.at("chr2")); +} + +// The BED contig list, not the truth VCF, decides which reference sequences are needed, and neither +// decides what may be freed: a contig absent from both is still reachable through a retained call. +TEST(IntersectContigs, NobedKeepsRefContigAbsentFromTruth) { + GlobalsGuard guard; + g.bed_exists = false; + std::shared_ptr query = make_variantData(QUERY, {"chr1"}, {12}, {{2}}); + std::shared_ptr truth = make_variantData(TRUTH, {"chr1"}, {12}, {{2}}); + std::shared_ptr ref = two_contig_ref(); + + testing::internal::CaptureStderr(); + intersect_contigs(query, truth, ref); + testing::internal::GetCapturedStderr(); + + EXPECT_EQ(size_t(2), ref->fasta.size()); + EXPECT_EQ("TTTTTTTTTTTT", ref->fasta.at("chr2")); +} + +// Retaining a contig means its reference sequence is now required, so a missing one is an error +// rather than a silent drop. Checked for both callsets, and whether or not a BED was supplied. +TEST(IntersectContigs, QueryContigMissingFromFastaErrors) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + EXPECT_EXIT({ + std::shared_ptr query = + make_variantData(QUERY, {"chr1", "chr2"}, {12, 12}, {{2}, {2}}); + std::shared_ptr truth = make_variantData(TRUTH, {"chr1"}, {12}, {{2}}); + std::shared_ptr ref = make_fasta("chr1", "ACGTACGTACGT"); + intersect_contigs(query, truth, ref); + }, testing::ExitedWithCode(1), + "Contig 'chr2' found in QUERY VCF but not reference FASTA"); +} + +TEST(IntersectContigs, BedTruthContigMissingFromFastaErrors) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + EXPECT_EXIT({ + std::shared_ptr query = make_variantData(QUERY, {"chr1"}, {12}, {{2}}); + std::shared_ptr truth = + make_variantData(TRUTH, {"chr1", "chr2"}, {12, 12}, {{2}, {2}}); + std::shared_ptr ref = make_fasta("chr1", "ACGTACGTACGT"); + intersect_contigs(query, truth, ref); + }, testing::ExitedWithCode(1), + "Contig 'chr2' found in TRUTH VCF but not reference FASTA"); +} - // the reference is pruned to the BED contigs too - EXPECT_EQ(size_t(1), ref->fasta.size()); - EXPECT_EQ(size_t(1), ref->fasta.count("chr1")); +// Nothing is evaluated on a contig outside the BED, so there are no false positives or false +// negatives there to warn about. +TEST(IntersectContigs, ContigOutsideBedDoesNotWarnOnMissingCounterpart) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + std::shared_ptr query = + make_variantData(QUERY, {"chr1", "chr2"}, {12, 12}, {{2}, {2}}); + std::shared_ptr truth = make_variantData(TRUTH, {"chr1"}, {12}, {{2}}); + std::shared_ptr ref = two_contig_ref(); + + testing::internal::CaptureStderr(); + intersect_contigs(query, truth, ref); + std::string out = testing::internal::GetCapturedStderr(); + + EXPECT_EQ(std::string::npos, out.find("Contig 'chr2'")) << out; +} + +// Ploidies are recorded before parse_variants() applies the BED filter, so a contig outside the BED +// carries observed ploidies that no evaluation ever consults. Disagreement there is not actionable. +TEST(IntersectContigs, ContigOutsideBedDoesNotWarnOnPloidy) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + std::shared_ptr query = + make_variantData(QUERY, {"chr1", "chr2"}, {12, 12}, {{2}, {1}}); + std::shared_ptr truth = + make_variantData(TRUTH, {"chr1", "chr2"}, {12, 12}, {{2}, {2}}); + std::shared_ptr ref = two_contig_ref(); + + testing::internal::CaptureStderr(); + intersect_contigs(query, truth, ref); + std::string out = testing::internal::GetCapturedStderr(); + + EXPECT_EQ(std::string::npos, out.find("contig 'chr2' has ploid")) << out; +} + +// A truth contig outside the BED still gets a query counterpart, because both callsets must remain +// indexable by every contig either one carries. Consumers index variants[hap][ctg] over the union +// of the two contig lists and dereference the result without a null check. +TEST(IntersectContigs, TruthOnlyContigOutsideBedIsPaired) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + std::shared_ptr query = make_variantData(QUERY, {"chr1"}, {12}, {{2}}); + std::shared_ptr truth = + make_variantData(TRUTH, {"chr1", "chr2"}, {12, 12}, {{2}, {1}}); + std::shared_ptr ref = two_contig_ref(); + + testing::internal::CaptureStderr(); + intersect_contigs(query, truth, ref); + std::string out = testing::internal::GetCapturedStderr(); + + // the query gains chr2, empty and having observed no ploidy of its own + EXPECT_EQ(std::vector({"chr1", "chr2"}), query->contigs); + EXPECT_EQ(std::vector({12, 12}), query->lengths); + EXPECT_EQ(std::vector< std::set >({{2}, {}}), query->observed_ploidies); + ASSERT_NE(nullptr, query->variants[HAP1]["chr2"]); + ASSERT_NE(nullptr, query->variants[HAP2]["chr2"]); + EXPECT_EQ(0, query->variants[HAP1]["chr2"]->n); + + // pairing it must not be mistaken for a ploidy disagreement, nor warned about as one-sided + EXPECT_EQ(std::string::npos, out.find("contig 'chr2' has ploid")) << out; + EXPECT_EQ(std::string::npos, out.find("Contig 'chr2'")) << out; +} + +// The mirror: a query contig outside the BED gains a truth counterpart. +TEST(IntersectContigs, QueryOnlyContigOutsideBedIsPaired) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 10}}); + std::shared_ptr query = + make_variantData(QUERY, {"chr1", "chr2"}, {12, 12}, {{2}, {1}}); + std::shared_ptr truth = make_variantData(TRUTH, {"chr1"}, {12}, {{2}}); + std::shared_ptr ref = two_contig_ref(); + + testing::internal::CaptureStderr(); + intersect_contigs(query, truth, ref); + std::string out = testing::internal::GetCapturedStderr(); + + EXPECT_EQ(std::vector({"chr1", "chr2"}), truth->contigs); + EXPECT_EQ(std::vector({12, 12}), truth->lengths); + EXPECT_EQ(std::vector< std::set >({{2}, {}}), truth->observed_ploidies); + ASSERT_NE(nullptr, truth->variants[HAP1]["chr2"]); + ASSERT_NE(nullptr, truth->variants[HAP2]["chr2"]); + EXPECT_EQ(0, truth->variants[HAP1]["chr2"]->n); + + EXPECT_EQ(std::string::npos, out.find("contig 'chr2' has ploid")) << out; + EXPECT_EQ(std::string::npos, out.find("Contig 'chr2'")) << out; } TEST(IntersectContigs, BedMissingInFastaErrors) { @@ -1166,7 +1306,7 @@ TEST(IntersectContigs, NobedTruthOnlyContigWarns) { EXPECT_EQ(0, query->variants[HAP2]["chr2"]->n); } -TEST(IntersectContigs, NobedFastaMissingErrors) { +TEST(IntersectContigs, NobedTruthContigMissingFromFastaErrors) { GlobalsGuard guard; g.bed_exists = false; EXPECT_EXIT({ @@ -1176,7 +1316,22 @@ TEST(IntersectContigs, NobedFastaMissingErrors) { std::shared_ptr ref = make_fasta("chr1", "ACGTACGTACGT"); intersect_contigs(query, truth, ref); }, testing::ExitedWithCode(1), - "Contig 'chr2' found in truth VCF but not reference FASTA"); + "Contig 'chr2' found in TRUTH VCF but not reference FASTA"); +} + +// Without a BED the missing sequence was previously reached as an unhandled std::out_of_range from +// lengths.at() while injecting the contig into the truth, rather than as a reported error. +TEST(IntersectContigs, NobedQueryContigMissingFromFastaErrors) { + GlobalsGuard guard; + g.bed_exists = false; + EXPECT_EXIT({ + std::shared_ptr query = + make_variantData(QUERY, {"chr1", "chr2"}, {12, 12}, {{2}, {2}}); + std::shared_ptr truth = make_variantData(TRUTH, {"chr1"}, {12}, {{2}}); + std::shared_ptr ref = make_fasta("chr1", "ACGTACGTACGT"); + intersect_contigs(query, truth, ref); + }, testing::ExitedWithCode(1), + "Contig 'chr2' found in QUERY VCF but not reference FASTA"); } TEST(IntersectContigs, PloidyMismatchWarns) { diff --git a/tests/unit/src/test_cluster.cpp b/tests/unit/src/test_cluster.cpp index 03e71f3..d2c0f58 100644 --- a/tests/unit/src/test_cluster.cpp +++ b/tests/unit/src/test_cluster.cpp @@ -10,6 +10,7 @@ #include "gtest/gtest.h" +#include "../../../src/bed.h" #include "../../../src/cluster.h" #include "../../../src/defs.h" #include "../../../src/globals.h" @@ -1698,6 +1699,40 @@ TEST(SuperclusterDataCtor, EndToEndSmoke) { EXPECT_EQ(std::vector({0}), tvars->superclusters); } +// The merge dereferences variants[hap][ctg] for every contig in the union of both contig lists, so +// a contig one callset never declared segfaults rather than being reported. intersect_contigs() is +// what rules that state out, by pairing every contig across the two callsets, and a BED covering +// only some of them must not be able to leave a contig unpaired. Runs the two together because +// neither function alone can show the invariant holds. +TEST(SuperclusterDataCtor, BedAbsentOneSidedContigIsPairedByIntersect) { + GlobalsGuard guard; + g.bed_exists = true; + g.bed = make_bed("chr1", {{0, 100}}); + + // chr2 is outside the BED and declared by the query alone, as a VCF header naming a contig the + // other callset's header omits entirely would leave it + std::shared_ptr qvd = + make_variantData(QUERY, {"chr1", "chr2"}, {1000, 500}, {{2}, {2}}); + std::shared_ptr tvd = make_variantData(TRUTH, {"chr1"}, {1000}, {{2}}); + std::shared_ptr ref = + make_fasta({{"chr1", std::string(1000, 'A')}, {"chr2", std::string(500, 'C')}}); + + testing::internal::CaptureStderr(); + intersect_contigs(qvd, tvd, ref); + testing::internal::GetCapturedStderr(); + + // the invariant the merge relies on: both callsets are indexable by every contig + ASSERT_NE(nullptr, tvd->variants[HAP1]["chr2"]); + ASSERT_NE(nullptr, tvd->variants[HAP2]["chr2"]); + + superclusterData sc_data(qvd, tvd, ref); + + EXPECT_EQ(std::vector({"chr1", "chr2"}), sc_data.contigs); + EXPECT_NE(nullptr, sc_data.superclusters["chr2"]->callset_vars[QUERY]); + EXPECT_NE(nullptr, sc_data.superclusters["chr2"]->callset_vars[TRUTH]); + EXPECT_EQ(0, sc_data.superclusters["chr2"]->callset_vars[TRUTH]->n); +} + /* wf_swg_cluster *********************************************************************************/ /**