Bcftools how to add DP to FORMAT field (get per sample read depth for REF vs ALT alleles )
1
0
Entering edit mode
9 weeks ago
Linda ▴ 30

I'm trying to achieve what this post was looking for Add Dp Tag To Genotype Field Of Vcf File

Currently this is my command:

bcftools mpileup -Ou --max-depth 8000 --min-MQ 30 --min-BQ 30 -f reference.fasta sample1.sorted.bam | bcftools call --ploidy 1 -Ou -mv |  bcftools filter -s LowQual -e \'%QUAL<20\'  > sample1.flt.vcf

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  sample1.sorted.bam
Imtechella_halotolerans_length_3113269  4051    .   C   T   41.4148 PASS    DP=2;VDB=0.02;SGB=-0.453602;MQ0F=0;AC=1;AN=1;DP4=0,0,0,2;MQ=42  GT:PL   1:71,0
Imtechella_halotolerans_length_3113269  4081    .   C   T   45.4146 PASS    DP=2;VDB=0.02;SGB=-0.453602;MQ0F=0;AC=1;AN=1;DP4=0,0,0,2;MQ=42  GT:PL   1:75,0

As far as I can understand the DP specified in the INFO is depth of coverage across all samples. How can I get the depth per sample info in the format/genotype field?

Any help appreciated!

bcftools snp vcf calling • 569 views
ADD COMMENT
1
Entering edit mode
9 weeks ago

Hi, you can specify this via the mpileup command. Please take a look at this snippet from my historical pipeline:

That is:

bcftools mpileup \
  --redo-BAQ \
  --min-BQ 30 \
  --per-sample-mF \
  --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
  -f "${Ref_FASTA}" \
  Aligned_Sorted_PCRDuped_FiltMAPQ.bam |\
  bcftools call \
    --multiallelic-caller \
    --variants-only \
    -Ob > Aligned_Sorted_PCRDuped_FiltMAPQ.bcf ;

Kevin

ADD COMMENT
1
Entering edit mode

Brilliant thanks Kevin

ADD REPLY
0
Entering edit mode

Hi Kevin, sorry to dredge up an old question.

I used your annotate code snippet above and now my FORMAT col is:

GT:PL:DP:SP:ADF:ADR:AD

Is it possible to add GQ - genotype quality to this field? I tried adding FORMAT/GQ to the code above

--annotate FORMAT/GQ,FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR

but it's saying could not parse tag "FORMAT/GQ"

Thanks so much! :)

ADD REPLY
0
Entering edit mode

Hi Linda, I was not sure; so, had to search. It seems that GQ is added during the subsequent bcftools call command:

-f, --format-fields <list>      output format fields: GQ,GP (lowercase allowed) []

Found via this answer: How to get GQ after running BCFTools?

ADD REPLY
1
Entering edit mode

Hey thanks for your help, I found that link too. All good!

ADD REPLY
0
Entering edit mode

Hey,

I used your annotate code

Is it possible to add AO and RO?

I used this command:

bcftools mpileup \ --annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/AO,FORMAT/RO,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \ -O b -o $OUTPUT_BCF_F1 \

and I got: Could not parse tag "FORMAT/AO"

ADD REPLY
0
Entering edit mode

I don't believe AO and RO are in any way standard / common tags. Where have you seen them previously? - freebayes? As far as I am aware, they relate to a somatic variant context; so, there is no way for bcftools to calculate them from just a VCF (need the original BAMs).

ADD REPLY
0
Entering edit mode

I am trying to run VcfSampleCompare, and I get warnings as you can see in this link: https://github.com/hepcat72/vcfSampleCompare/issues/26

Robert is sayings that my vcf file is missing DP,AO,RO tags and I need to add the DP tag. After adding the DP tag, usings your script , I had another warning so I think maybe I need the AO RO too..

ADD REPLY

Login before adding your answer.

Traffic: 1730 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6