diff --git a/src/AlleleParser.cpp b/src/AlleleParser.cpp index 27af337f..1964282e 100644 --- a/src/AlleleParser.cpp +++ b/src/AlleleParser.cpp @@ -2182,9 +2182,9 @@ 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 @@ -2192,9 +2192,9 @@ void AlleleParser::getInputVariantsInRegion(string& seq, long start, long end) { // 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 index 962c6287..18c574a3 100644 --- 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 new file mode 100644 index 00000000..87d557e0 Binary files /dev/null and b/test/tiny/q.vcf.gz differ diff --git a/test/tiny/q.vcf.gz.tbi b/test/tiny/q.vcf.gz.tbi new file mode 100644 index 00000000..07810148 Binary files /dev/null and b/test/tiny/q.vcf.gz.tbi differ