This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Calculate Allele balance by sample

I have a VCF file which variants called and has several field like AD, DP, AF, AC, etc. I however also need to calculate the AB (allele balance for each sample) which the information available in this vcf file. I tried looking for a softwae online but I couldn't find any.

Here is a sample output of the vcf file:

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S1

X 5867701 . T C,G 231.77 . AC=1,1;AF=0.500,0.500;AN=2;DP=8;ExcessHet=3.0103;FS=0.000;MLEAC=1,1;MLEAF=0.500,0.500;MQ=45.48;QD=28.97;SOR=1.179 GT:AD:DP:GQ:PL 1/2:0,4,4:8:96:260,167,153,96,0,151

X 7753630 . A T,C 469.77 . AC=1,1;AF=0.500,0.500;AN=2;DP=13;ExcessHet=3.0103;FS=0.000;MLEAC=1,1;MLEAF=0.500,0.500;MQ=60.00;QD=30.18;SOR=2.833 GT:AD:DP:GQ:PL 1/2:0,10,3:13:88:498,88,96,416,0,407

Thanks in advance.

allele balance snp rna-seq
• 7,612 views • link
updated by ashotmarg2004 • ashotmarg2004 • written by kirannbishwa01 kirannbishwa01

2 answers

Not sure if this is the most elegant and efficient way but I tried with bcftools like this for the case of my tri-allelic sites, similar to the ones you have.

bcftools view -m3 -M3 someBCFfile.bcf | bcftools filter -e '(GT=\"het\" & FORMAT/AD[0:1]/FORMAT/AD[0:2]>4) || (GT=\"het\" & FORMAT/AD[0:1]/FORMAT/AD[0:2]<0.25)' -Ob -o outputBCFfile.bcf

The first part (view -m3 -M3) simply extracts the tri-allelic sites The second part (filter section) excludes heterozygous calls (GT=\"het\") that have allelic balance (AB) over 4 or less than 0.25 (this was just my choice, you should put your own values of course). In your example site: 5867701 will have AB of 4/4=1 while the second site: 7753630 will have AB of 10/3=3.33, i.e. both sites will be kept.

For the little bit tricky notation of FORMAT/AD[0:1]/FORMAT/AD[0:2]>4 check the expressions section of the bcftools's website, http://samtools.github.io/bcftools/bcftools.html#expressions

For bi-allelic site one should simply do FORMAT/AD[0:0]/FORMAT/AD[0:1]>4 instead to parse the correct values in the FORMAT/AD field.

• 0 views • link
written by ashotmarg2004 • ashotmarg2004

After I wrote my "manual" bcftools version, I've actually found this: https://gatk.broadinstitute.org/hc/en-us/articles/360040507131-FilterVcf-Picard- Seems that the --MIN_AB argument might just do what you are looking for, but I haven't used it so cannot comment on it.

• 0 views • link
written by ashotmarg2004 • ashotmarg2004

Log in to answer this question.

More posts like this