bcftools filter using "." for missing values in integer or float fields
Nobody has claimed this yet.
Assessment
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Newbie friendliness
- 35/100
- Issue type
- Bug
- Clarity
- Mostly clear
- Activity status
- Stale
- Tech stack
- c
- Domain
- bioinformatics
Research direction
Start at the bcftools filter command and the expression handling for FMT/DP comparisons. Reproduce the reported filter on the merged VCF example, then verify that “.” matches missing Integer and Float values as it does for String fields, producing the expected genotype changes.
Written by the indexing model from the issue text.
Description
In bcftools 1.3 I cannot filter using "." to search for missing values in FORMAT fields which are defined as Integer or Float. According to the man page one can filter for missing values and the following examples are given:
"." to test missing values
DP=".", DP!=".", ALT="."
This would be useful, e.g. when merging VCF files. For instance, in VCF file [a] we excluded positions which were called as reference, keeping only high quality SNPs, and any positions where we had zero coverage, or where the confidence in the reference call was also low (as we might not want to assume a reference genotype). In VCF file [b] we do the same. Now for positions where [a] does not have a record but [b] does we might end up with a FORMAT field which looks like this (for 2 samples and 3 positions):
1 GT:AD:DP:RGQ ./.:.:.:. ./.:0:0:0
2 GT:AD:DP:RGQ ./.:.:.:. ./.:0:0:0
3 GT:AD:DP:RGQ ./.:.:.:. ./.:0:0:0
The first format field is ./.:.:.:. because it was missing from VCF [a], so we know it was a 0/0 genotype. The second format field is ./.:0:0:0. The additional information in the DP tag tells us that this position had zero coverage so it is good that the genotype is missing.
It would therefore be good to be able to filter these to set the genotype of ./.:.:.:. back into reference. We should be able to achieve this by the following filter command:
bcftools filter -S0 -e 'FMT/GT="." & FMT/DP = "."' merged.vcf
However bcftools will not look for missing values using the "." in Integer or Float fields. If I manually edit the header from:
FORMAT=<ID=DP,Number=1,Type=Integer,Description="")">
to:
FORMAT=<ID=DP,Number=1,Type=String,Description="")">
The filtering will work as expected and the above filter command returns:
1 GT:AD:DP:RGQ 0/0:.:.:. ./.:0:0:0
2 GT:AD:DP:RGQ 0/0:.:.:. ./.:0:0:0
3 GT:AD:DP:RGQ 0/0:.:.:. ./.:0:0:0
- Dominant language
- C
- Stars
- 891
- Forks
- 277
- Avg merge
- 2d 23h
- Merged PRs (30d)
- 2
Contributor guide
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
More from samtools/bcftools
-
needs-feedback
Difficulty 3/5 1-2 days Newbie friendliness 55/100
-
enhancement
Difficulty 5/5 Over a week Newbie friendliness 45/100
-
htslib-dependent P2: Desirable
Difficulty 3/5 1-2 days Newbie friendliness 48/100
-
Difficulty 5/5 Over a week Newbie friendliness 30/100
-
needs-feedback requires-test-case
Difficulty 3/5 1-2 days Newbie friendliness 55/100
All issues in samtools/bcftools
Similar issues
-
[adam] AdamNet network read doesn't cap to MAX_ADAM_PACKET_LEN, overflows client receive buffers Open
Difficulty 2/5 1-3 hours Newbie friendliness 88/100
FujiNetWIFI/fujinet-firmware#1649 · 2 comments ·
-
Difficulty 2/5 1-3 hours Newbie friendliness 68/100
HarbourMasters/Shipwright#7229 ·
-
Difficulty 2/5 1-3 hours Newbie friendliness 72/100
riscv-software-src/riscv-isa-sim#2435 · 1 comment ·
-
bug Self Built Image SNAPSHOT Supported Device target/ramips
Difficulty 2/5 1-3 hours Newbie friendliness 88/100
-
Difficulty 2/5 1-3 hours Newbie friendliness 76/100