refactor(bed)!: stop pruning contigs and reference sequence in intersect_contigs (#246)

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).

src/bed.cpp

@@ -481,14 +481,49 @@ static std::string ploidy_set_str(const std::set<int> & 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<fastaData> & ref_ptr,
+        const std::shared_ptr<variantData> & from_ptr,
+        const std::shared_ptr<variantData> & 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<ctgVariants>(new ctgVariants(ctg));
+        to_ptr->variants[HAP2][ctg] = std::shared_ptr<ctgVariants>(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<variantData> 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<std::string>::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<callset_t, std::shared_ptr<variantData> > > 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<ctgVariants>(new ctgVariants(ctg));
-                truth_ptr->variants[HAP2][ctg] = 
-                        std::shared_ptr<ctgVariants>(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<ctgVariants>(new ctgVariants(ctg));
-                query_ptr->variants[HAP2][ctg] = 
-                        std::shared_ptr<ctgVariants>(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<int> & truth_ploidies = truth_ptr->observed_ploidies[truth_ctg_idx];
         const std::set<int> & query_ploidies = query_ptr->observed_ploidies[query_ctg_idx];

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,

tests/integration/data/undeclared_contig_query.vcf

@@ -0,0 +1,8 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##contig=<ID=sc1,length=400>
+##contig=<ID=sc2,length=100>
+#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

tests/integration/data/undeclared_contig_truth.vcf

@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##contig=<ID=sc1,length=400>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE
+sc1	200	.	A	G	50	PASS	.	GT	1/1