bcftools roh fails for samples used to generate the reference genome
Nobody has claimed this yet.
Assessment
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Newbie friendliness
- 35/100
- Issue type
- Bug
- Clarity
- Needs clarification
- Activity status
- Stale
- Tech stack
- c
- Domain
- bioinformatics
Research direction
Start at the bcftools roh HMM likelihood calculations described in the issue and review how reference-derived genotypes affect the RR, RA, and AA terms. Reproduce the missing-ROH behavior with a reference individual, then determine and test a documented or implemented way to handle this genotype distribution.
Written by the indexing model from the issue text.
Description
Hi,
I have been using bcftools roh to estimate runs of homozygosity. I noticed that the individual used to generate the reference genome sequence had no roh calls, despite the fact that I have visualised a roh in this individual previously.
This individual will almost always have the genotypes 0/0 or 0/1 at a given site, as 1/1 would suggest an error in the genome assembly or variant calling.
The bcftools roh HMM calculates the likelihood of the hidden states (N = non-roh, H = roh) as follows:
P(Di|Xi=H)=(1–fi)P(Di|RR)+fiP(Di|AA)
P(Di|Xi=N)=(1–fi)2P(Di|RR)+2fi(1–fi)P(Di|RA)+fi2P(Di|AA)
My problem, I think, is that I have almost zero 'AA' genotypes (i.e. 1/1) when using the reference individual, so these probabilities are not meaningful.
Do you have any suggestions to get around this problem? One option would be to make a pseudo-vcf where alleles in this individual are flipped at a frequency of 0.4 (the default AF value).
Any advice would be very welcome,
Alex
- 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