correctly handle --variant-input

freebayes · May 21, 2015 · c15a283 · c15a283
1 parent 06ea37a
commit c15a283
Show file tree

Hide file tree

Showing 4 changed files with 8 additions and 5 deletions.
diff --git a/src/AlleleParser.cpp b/src/AlleleParser.cpp
@@ -2182,19 +2182,19 @@ void AlleleParser::getInputVariantsInRegion(string& seq, long start, long end) {
                     allelePos -= 1;
                     reflen = len + 2;
                     alleleSequence =
-                        uppercase(reference.getSubSequence(currentSequenceName, allelePos, 1))
+                        uppercase(reference.getSubSequence(currentVariant->sequenceName, allelePos, 1))
                         + alleleSequence
-                        + uppercase(reference.getSubSequence(currentSequenceName, allelePos+1+len, 1));
+                        + uppercase(reference.getSubSequence(currentVariant->sequenceName, allelePos+1+len, 1));
                     cigar = "1M" + convert(len) + "D" + "1M";
                 } else {
                     // we always include the flanking bases for these elsewhere, so here too in order to be consistent and trigger use
                     type = ALLELE_INSERTION;
                     // add previous base and post base to match format typically used for calling
                     allelePos -= 1;
                     alleleSequence =
-                        uppercase(reference.getSubSequence(currentSequenceName, allelePos, 1))
+                        uppercase(reference.getSubSequence(currentVariant->sequenceName, allelePos, 1))
                         + alleleSequence
-                        + uppercase(reference.getSubSequence(currentSequenceName, allelePos+1, 1));
+                        + uppercase(reference.getSubSequence(currentVariant->sequenceName, allelePos+1, 1));
                     len = variant.alt.size() - var.ref.size();
                     cigar = "1M" + convert(len) + "I" + "1M";
                     reflen = 2;

diff --git a/test/t/01_call_variants.t b/test/t/01_call_variants.t
@@ -5,7 +5,7 @@ BASH_TAP_ROOT=bash-tap
 
 PATH=../bin:$PATH # for freebayes
 
-plan tests 5
+plan tests 6
 
 is $(echo "$(comm -12 <(cat tiny/NA12878.chr22.tiny.giab.vcf | grep -v "^#" | cut -f 2 | sort) <(freebayes -f tiny/q.fa tiny/NA12878.chr22.tiny.bam | grep -v "^#" | cut -f 2 | sort) | wc -l) >= 13" | bc) 1 "variant calling recovers most of the GiAB variants in a test region"
 
@@ -67,3 +67,6 @@ is $(samtools view -u tiny/NA12878.chr22.tiny.bam | freebayes -f tiny/q.fa --std
     $(freebayes -f tiny/q.fa tiny/NA12878.chr22.tiny.bam | grep -v "^#" | wc -l) "reading from stdin or not makes no difference"
 
 is $(samtools view tiny/NA12878.chr22.tiny.bam | wc -l) $(freebayes -f tiny/q.fa tiny/NA12878.chr22.tiny.bam -d 2>&1 | grep ^alignment: | wc -l) "freebayes processes all alignments in input"
+
+# ensure targeting works even when there are no reads
+is $(freebayes -f tiny/q.fa -@ tiny/q.vcf.gz tiny/NA12878.chr22.tiny.bam | grep -v "^#" | wc -l) 19 "freebayes correctly handles variant input"
diff --git a/test/tiny/q.vcf.gz b/test/tiny/q.vcf.gz
diff --git a/test/tiny/q.vcf.gz.tbi b/test/tiny/q.vcf.gz.tbi