This is a test version of Biostars. For the public version, visit https://www.biostars.org.
heterozygosity filter through bcftools

Hi, I have to filter out variants sites with heterozygosity above 20% (0.2) and retain only sites with heterozygosity below 20%.

Can anyone kindly guide which bcftools command can be utilized?

Best regards,

heterozygosity bcftools
• 396 views • link
written by Bushra Bibi 8

With bcftools, it's always this: "filter" = bcftools view + appropriate include/exclude filters

• 0 views • link
written by Ram 4.5K

Thanks Ram!

Is it the correct one?

bcftools view \
    -i 'COUNT(GT="het")/(N_SAMPLES-COUNT(GT="mis")) <= 0.2' \
    -Oz \
    -o polymorphic.vcf.gz \
    raw_variants.vcf.gz
• 0 views • link
written by Bushra Bibi 8

The syntax looks fine. I cannot mentally parse your filter, you can judge its accuracy by trial and error if you establish proper test cases.

• 0 views • link
written by Ram 4.5K

1 answer

COUNT isn't the right function there - it aggregates over the values inside a tag, not over samples. To count samples matching a genotype condition you want N_PASS, and there's a built-in for the missing count as well:

bcftools view -i 'N_PASS(GT="het")/(N_SAMPLES-N_MISSING) <= 0.2' -Oz -o out.vcf.gz in.vcf.gz

Rather than trial and error on the whole file, bcftools query -f '%CHROM\t%POS[\t%GT]\n' on a dozen sites and counting the hets by hand is quicker for checking the denominator does what you meant.

• 0 views • link
updated by GenoMax 16K • | written by SeqBench 41

Thank you SeqBench.

My goal is to retain only those snps where heterozygosity is less than 20% and filter out snps not falling under threshold.

I do not want to remove samples but snps only.

Is N_PASS a filter for samples or SNPs?

I want some code that behaves as alternative of following code of tassel.

tassel-5-standalone/run_pipeline.pl -Xmx40g -vcf raw.vcf -FilterSiteBuilderPlugin -maxHeterozygous 0.2 -endPlugin -export "het_filtered" -exportType VCF

• 0 views • link
written by Bushra Bibi 8

SNPs. N_PASS counts samples at a site, but the whole expression is evaluated once per record, and bcftools view -i keeps or drops entire lines - it never touches the sample columns. Nobody gets removed, you just lose sites.

One thing I'd check rather than take from me: whether TASSEL's maxHeterozygous divides by all taxa or only those with a call. Mine uses non-missing as the denominator. If TASSEL uses all taxa, swap it for F_PASS(GT="het") <= 0.2, which divides by the full sample count.

• 0 views • link
written by SeqBench 41

I think right approach would be non-missing as a denominator. I will go with N_PASS

Thanks alot SeqBench

• 0 views • link
written by Bushra Bibi 8

Log in to answer this question.

More posts like this