AlistairNWard/vcfCTools — Issue #4: Not all variants filtered with multiple records of same variant type at same position (with patch)
Author: kmsquire
Created: 2011-12-06 23:06:19+00:00
State: open

When two vcf files have very different sizes, one buffer might be starved, leaving only one vcf record in its variant buffer. If a second variant of the same type at the same location is still not loaded, the variants at that position will be filtered incorrectly.

Solution was to split out updateVariantStructure() from buildVariantStructure(), and use it to keep the buffer full. Patch below, also available in commit acdbceb116181d90cd5d6acfad1df6971c0dbfc1.

```diff
diff --git a/intersect.cpp b/intersect.cpp
index 092f92e..908d01a 100644
--- a/intersect.cpp
+++ b/intersect.cpp
@@ -80,17 +80,11 @@ void intersect::intersectVcf(vcfHeader& header1, vcfHeader& header2, vcf& v1, va
// Clear the compared variants from the structure and add the next one from
// the file into the structure if it is from the same reference sequence.
var1.variantMap.erase(var1.vmIter);
- if (v1.variantRecord.referenceSequence == currentReferenceSequence && v1.success) {
- var1.addVariantToStructure(v1.position, v1.variantRecord);
- v1.success = v1.getRecord();
- }
+ v1.success = var1.updateVariantStructure(v1, currentReferenceSequence);
if (var1.variantMap.size() != 0) {var1.vmIter = var1.variantMap.begin();}

var2.variantMap.erase(var2.vmIter);
- if (v2.variantRecord.referenceSequence == currentReferenceSequence && v2.success) {
- var2.addVariantToStructure(v2.position, v2.variantRecord);
- v2.success = v2.getRecord();
- }
+ v2.success = var2.updateVariantStructure(v2, currentReferenceSequence);
if (var2.variantMap.size() != 0) {var2.vmIter = var2.variantMap.begin();}

// Variant from the first vcf file is at a larger coordinate than that in the
@@ -99,10 +93,7 @@ void intersect::intersectVcf(vcfHeader& header1, vcfHeader& header2, vcf& v1, va
} else if (var1.vmIter->first > var2.vmIter->first) {
if (flags.findCommon) {var2.filterUnique();}
var2.variantMap.erase(var2.vmIter);
- if (v2.variantRecord.referenceSequence == currentReferenceSequence && v2.success) {
- var2.addVariantToStructure(v2.position, v2.variantRecord);
- v2.success = v2.getRecord();
- }
+ v2.success = var2.updateVariantStructure(v2, currentReferenceSequence);

// Reset the iterator if there are still variants in the structure.
if (var2.variantMap.size() != 0) {
@@ -130,10 +121,7 @@ void intersect::intersectVcf(vcfHeader& header1, vcfHeader& header2, vcf& v1, va
} else if (var1.vmIter->first < var2.vmIter->first) {
if (flags.findCommon && !flags.annotate) {var1.filterUnique();}
var1.variantMap.erase(var1.vmIter);
- if (v1.variantRecord.referenceSequence == currentReferenceSequence && v1.success) {
- var1.addVariantToStructure(v1.position, v1.variantRecord);
- v1.success = v1.getRecord();
- }
+ v1.success = var1.updateVariantStructure(v1, currentReferenceSequence);

// Reset the iterator if there are still variants in the structure.
if (var1.variantMap.size() != 0) {
diff --git a/variant.cpp b/variant.cpp
index d5b677e..b63e1fd 100644
--- a/variant.cpp
+++ b/variant.cpp
@@ -72,35 +72,44 @@ void variant::determineVariantsToProcess(bool snps, bool mnps, bool indels, bool

// Build up a structure containing variants.
bool variant::buildVariantStructure(vcf& v) {
- unsigned int count = 0;

// When variants in the correct reference sequence are found, build the variant
// structure.
- count = 0;
string tempReferenceSequence = v.variantRecord.referenceSequence;

- // If this is not the first time that this reference sequence has been seen
- // when building, the structure then the records for this reference sequence
- // are not contiguous in the vcf file. For a number of tools (especially
- // the intersect tool), this will cause erroneous results to be generated.
- if (referenceSequenceInfo.count(tempReferenceSequence) != 0) {
- referenceSequenceInfo[tempReferenceSequence].contiguous = false;
- }
- while (v.success && count < recordsInMemory && v.variantRecord.referenceSequence == tempReferenceSequence) {
+ if (v.success && variantMap.size() < recordsInMemory) {

+ // If this is not the first time that this reference sequence has been seen
+ // when building, the structure then the records for this reference sequence
+ // are not contiguous in the vcf file. For a number of tools (especially
+ // the intersect tool), this will cause erroneous results to be generated.
+ if (referenceSequenceInfo.count(tempReferenceSequence) != 0) {
+ referenceSequenceInfo[tempReferenceSequence].contiguous = false;
+ }
+ else {
// Update the information about observed reference sequences.
- if (referenceSequenceInfo.count(tempReferenceSequence) == 0) {
- referenceSequenceInfo[tempReferenceSequence].numberRecords = 1;
+ referenceSequenceInfo[tempReferenceSequence].numberRecords = 0;
referenceSequenceInfo[tempReferenceSequence].usedInComparison = false;
referenceSequenceInfo[tempReferenceSequence].contiguous = true;
- } else {
- referenceSequenceInfo[tempReferenceSequence].numberRecords++;
}

+ return updateVariantStructure(v, tempReferenceSequence);
+ }
+
+ return v.success;
+}
+
+// Update the structure containing variants.
+bool variant::updateVariantStructure(vcf& v, string& currentReferenceSequence) {
+
+ while (v.success && variantMap.size() < recordsInMemory && v.variantRecord.referenceSequence == currentReferenceSequence) {
+
+ // Update the information about observed reference sequences.
+ referenceSequenceInfo[currentReferenceSequence].numberRecords++;
+
// Add the variant into the structure.
addVariantToStructure(v.position, v.variantRecord);
v.success = v.getRecord();
- count++;
}

return v.success;
diff --git a/variant.h b/variant.h
index da900b9..bfc9f37 100644
--- a/variant.h
+++ b/variant.h
@@ -124,6 +124,7 @@ class variant {
void annotateRecordVcf(variant&, int, unsigned int, string&, string&);
void buildOutputRecord(output&, vcfHeader&);
bool buildVariantStructure(vcf&);
+ bool updateVariantStructure(vcf&, string&);
void clearOriginalVariants(vcfHeader&, intFlags&, output&, bool);
void clearReferenceSequence(vcfHeader&, vcf&, intFlags, string, output&, bool);
void clearReferenceSequenceBed(vcfHeader&, vcf&, intFlags, string, output&);
```