Significant difference between IDV and AD when calling certain RNA-seq indels with mpileup
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
- Domain
- bioinformatics
Research direction
Start by reproducing the supplied bcftools mpileup command at 11:69587265 and compare its IDV, INFO/AD, FORMAT/AD, and DP values with the IGV view. Investigate how spliced RNA-seq reads carrying the insertion are treated, then document whether the difference is expected behavior or a defect and provide a focused regression case if a code change is needed.
Written by the indexing model from the issue text.
Description
I'm using mpileup to compute the mutation frequencies of a list of known variant sites in RNA-seq. For a small number of indels the AD numbers in the output are much lower than what I'm seeing in IGV. For instance, when I run the following command:
bcftools mpileup -r 11:69587265 -f <GRCm38 fasta> --annotate FORMAT/AD,FORMAT/DP,INFO/AD -F 0.001 --max-depth 10000 --max-idepth 10000 -Q20 -x -A --no-BAQ --tandem-qual 10000 <RNA-seq bam> | grep -v "^#"
Here is the output I'm getting:
11 69587265 . CA CAA 0 . INDEL;IDV=77;IMF=0.616;DP=124;AD=47,2;I16=24,23,1,1,1880,75200,80,3200,872,17236,40,800,1057,25169,46,1066;QS=0.908467,0.0915332;VDB=0.56;SGB=-0.453602;RPBZ=-0.571191;MQBZ=2.56432;MQSBZ=0;BQBZ=0;SCBZ=-1.87886;MQ0F=0 PL:DP:AD 0,106,55:49:47,2
And here is how this site looks like in IGV:
I am aware that some low quality reads are filtered during variant calling and the reported AD number will be lower than the raw IDV. However upon manual examination most of the reads seem to be high-quality and correctly aligned. Such a huge drop from 77 to 2 seems counter-intuitive to me and I wonder if this is the expected behavior of mpileup.
(One thing I notice is that most of the read pairs containing the insertion are spliced except 2, which happens to be the same as the AD reported)
- 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