AlistairNWard/vcfCTools — Issue #5: Files incorrectly parsed when one contains contigs not present in the other (with patch)
Author: kmsquire
Created: 2011-12-06 23:14:47+00:00
State: open

In intersectVcf (intersect.cpp), comments suggest that parsing will stop when the end of one file is reached and the variant structure for that file is empty. However, the code only stops when the end of both files are reached and the variant structures are both empty. This patch fixes that, and then also flushes the rest of the records in a vcf when that vcf contains contigs at the end of the vcf not present in the other vcf file. While this may be caused by misordered vcf files, it may also simply be that one or the the other vcf files does not include variants for certain contigs (e.g., MT or non-standard contigs). Patch below, also available as commit 081153a6d6f4763030ed49afc692cdfcfcb49811.

```diff

diff --git a/intersect.cpp b/intersect.cpp
index 908d01a..825d1dd 100644
--- a/intersect.cpp
+++ b/intersect.cpp
@@ -56,7 +56,7 @@ void intersect::intersectVcf(vcfHeader& header1, vcfHeader& header2, vcf& v1, va

// Parse and compare the two variant structures until the end of one of the files
// is reached and the variant structure for that file is empty.
- while ( !(var1.variantMap.size() == 0 && !v1.success) || !(var2.variantMap.size() == 0 && !v2.success) ) {
+ while ( !(var1.variantMap.size() == 0 && !v1.success) && !(var2.variantMap.size() == 0 && !v2.success) ) {

// If the two variant structures are built with the same reference sequence, compare
// the contents and parse through all varians for this reference sequence.
@@ -212,6 +212,26 @@ void intersect::intersectVcf(vcfHeader& header1, vcfHeader& header2, vcf& v1, va
}
}

+ // If one vcf file contains extra contigs at the end, flush these when finding the union or
+ // variants unique to that file
+ if (var1.variantMap.size() != 0 && var2.variantMap.size() == 0 &&
+ (flags.findUnion || (flags.findUnique && flags.writeFromFirst))) {
+ while(var1.variantMap.size() != 0) {
+ var1.clearReferenceSequence(header1, v1, flags, currentReferenceSequence, ofile, flags.writeFromFirst);
+ var1.buildVariantStructure(v1);
+ if (var1.variantMap.size() != 0) {var1.vmIter = var1.variantMap.begin();}
+ }
+ }
+ else if (var1.variantMap.size() == 0 && var2.variantMap.size() != 0 &&
+ (flags.findUnion || (flags.findUnique && !flags.writeFromFirst)) &&
+ !flags.annotate) {
+ while(var2.variantMap.size() != 0) {
+ var2.clearReferenceSequence(header2, v2, flags, currentReferenceSequence, ofile, !flags.writeFromFirst);
+ var2.buildVariantStructure(v2);
+ if (var2.variantMap.size() != 0) {var2.vmIter = var2.variantMap.begin();}
+ }
+ }
+
// If the variant structures are not empty, there was a problem and a warning is given.
if (var1.variantMap.size() != 0 || (var2.variantMap.size() != 0 && !flags.annotate) ) {
cerr << "WARNING: Not all records were flushed out of the variant structure." << endl;
```