Bcftools how to add DP to FORMAT field (get per sample read depth for REF vs ALT alleles )

Question

[Linda]

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!

Answer 1

[Kevin Blighe]

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

• ClinicalGradeDNAseq#L300-L316

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

[Linda]

Brilliant thanks Kevin

[Linda]

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! :)

[Kevin Blighe]

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?

[Linda]

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

[Michal Nevo]

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"

[Kevin Blighe]

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).

[Michal Nevo]

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..