Closed haochenz96 closed 1 year ago
It is possible that indel AD counts will be improved in this pull request https://github.com/samtools/bcftools/pull/1679
Thank you! I assume I would clone from here? https://github.com/jkbonfield/bcftools/tree/indel-tweak-jkb1
That's correct. The pull request will be integrated into develop, but may take a while.
Hey! Thanks again for the info. While waiting for the pull request, I came up with a temporary (ugly) solution. See here my code:
# @STEP1 mpileup and call
bcftools mpileup \
-Ou \
-R {params.PANEL_BED} \
-f {params.REF_GENOME} \
--annotate FORMAT/AD,FORMAT/DP,INFO/AD \
--max-depth 100000 \
--max-idepth 100000 \
--no-BAQ \
{input.SC_BAM} | \
bcftools call \
--keep-alts \
-C alleles \
-T {input.CANDIDATE_ALLELE} \
--multiallelic-caller \
-Oz \
-o {output.SC_MPILEUP_VCF} && \
tabix {output.SC_MPILEUP_VCF}
# @STEP2 get raw counts
# extract INFO/AD into a tab-delimited annotation file
SAMPLE_NAME=$(bcftools query -l {output.SC_MPILEUP_VCF})
bcftools query -f '%CHROM\t%POS\t%AD{{0}},%AD{{1}}\n' {output.SC_MPILEUP_VCF} | bgzip -c > $OUT_DIR/annot.txt.gz && \
tabix -s1 -b2 -e2 $OUT_DIR/annot.txt.gz
echo -e '##FORMAT=<ID=AD,Number=1,Type=Integer,Description="Total allelic depths (high-quality bases)">' >> $OUT_DIR/hdr.txt
bcftools annotate \
-s $SAMPLE_NAME \
-a $OUT_DIR/annot.txt.gz \
-h $OUT_DIR/hdr.txt \
q \
-Oz \
-o {output.SC_MPILEUP_VCF_raw} \
{output.SC_MPILEUP_VCF}
Essentially I let the bcftools call
proceed as it does, and then I extract the raw REF & ALT read counts (first two fields of AD
) for every mutated site of interest from the VCF output. I assume this is fine since every genomic site I input is biallelic. Do you think this will work?
Have you been able to apply your workaround successfully, @haochenz96? I see the same issue in my data. However, your workaround does not work for me. I am looking at simple 1-nt insertions or deletions. In all these cases, AC=0
and the INFO/AD
field does not show any depth for an alternate allele. My data is hybrid capture, but that should not make a difference.
Maybe your CC/AT
genotype is a special case?
Example:
bcftools mpileup -B -d 10000 -O u -f genome.fa -R regions.bed --annotate FORMAT/AD,FORMAT/DP,INFO/AD sample.bam | \
bcftools call -A -O v -C alleles -m -T als.tsv.gz > bcftools.vcf
1 46670206 . TC T 281.579 . DP=839;AD=588,0;MQ0F=0;AC=0;AN=2;DP4=312,276,0,0;MQ=60 GT:PL:DP:AD 0/0:0,255,255:588:588,0
Hi Ben- check out this issue #1711 regarding multiallelic site. Yours looks like a multiallelic site to me as the DP
in FORMAT
does not match that in INFO
.
Otherwise, feel free to upload a test case.
@gabeng Indels may be missed because of high depth, try increasing -L
-L, --max-idepth INT Maximum per-file depth for INDEL calling [250]
Unfortunately, the -L
option does not change the result.
I am attaching a small test case with three Indels that should be called.
OPS-3885.zip
Really? I am using the latest github version on the develop branch and this is what I am observing with your data:
$ bcftools mpileup -f hs37d5.fa test.bam -r 1:46670206 -L 1000 -a AD,INFO/AD |
bcftools call -mv |
bcftools query -f '%CHROM:%POS %REF %ALT %AD [ %GT] [ %AD]\n'
1:46670206 TCC TC 161,166 0/1 161,166
You are absolutely right, when simply calling any variants in the data, the Indels (or most of them) are correctly called.
However, this issue is about the (forced) genotyping at predefined coordinates using the options -C
and -T
.
I see it now, thanks for the test case. This should be fixed by https://github.com/samtools/bcftools/commit/815fd54e261ee400010082129ff860ff7023c5d0.
By the way, please open a separate issue next time and reference a related issue rather than append to an existing one.
Sorry, I honestly thought the original issues would be the same as mine because the command line options were similar. I did not notice the two AD values of the Indel in the original post. Thank you for the quick fix. I tried it out and have two questions:
I would have expected that the Indels are called exactly as in the input als.tsv, however, they are not:
# output of current develop branch
1 46670206 . TCC TC
1 172328767 . T TA
1 204502514 . TTCT TTCTGAAACAGGGTCT
# output of bcftools-1.16
1 46670206 . TC T
1 172328767 . T TA
1 204502514 . T TTCTGAAACAGGG
Should the output of bcftools mpileup/call
generally be normalized downstream?
Another thing I noticed (and this may be unrelated): when I use smaller target regions with the develop version that are only 25bp +/- around the Indels in the -R
option, 2 of 3 Indels are not called. The BED file and the VCF file are attached:
original-targets.zip
Is 25bp padding around Indels not sufficient?
Normalization is a separate step, yes, just stream through bcftools norm
. Regarding the difference, maybe the 1.16 version just inserted a missing site via the -i option that it could not match, not sure.
I was not able to reproduce the problem with the missing sites 1:172328767 and 1:204502514. If it's really not working for you with the new develop version, please open a new issue and include the full command.
Hi,
Given a list of known alleles, I was trying to genotype REF/ALT read counts in some samples' BAM files. Relevant to the issue below is that these are libraries generated by targeted amplicon-sequencing. I used the following command:
But for ALL of the MNP/Indels in the given allele list, the output VCF contains 0 ALT read in
FORMAT/AD
despite they show up in theINFO/AD
field, as shown for thechr12:25398284
locus below:Note that this applies to every indel/MNP so is very likely caused by some call setting not compatible with the sequencing technology. I feel like this could be a easy fix by disabling some base/mapping quality filters, since my goal is to simply get the REF/ALT counts for each known allele. I would appreciate any help!