How to count SNPs, InDels

Question

[chchl7]

I have used Bowtie 2 to align my reads to reference genome and Samtools to call variances (SNPs and InDels). Does anyone know how to count how many SNPs and InDels I got? Thanks!

[Pierre Lindenbaum]

in a BAM or in a VCF ?

Answer 1

[Vivek]

You could write a quick awk line ignoring multi-allelic loci:

SNPs:

awk '! /\#/' variants.vcf | awk '{if(length($4) == 1 && length($5) == 1) print}' | wc -l

Indels:

awk '! /\#/' variants.vcf | awk '{if(length($4) > 1 || length($5) > 1) print}' | wc -l

but for something more comprehensive, use bcftools stats.

[scchess]

Alex's solution is faster and better.

[chchl7]

Thank you very much! The scripts are working well. How can I use awk to count Insertions and deletions separately? Thank you!

[nchuang]

insertions:

awk '! /\#/' variants.vcf | awk '{if(length($4) > 1 ) print}' | wc -l

deletions:

awk '! /\#/' variants.vcf | awk '{if(length($5) > 1) print}' | wc -l

[Pierre Lindenbaum]

Wrong. If the variation is REF=A ALT=T,G or if the variation is symbolic, e.g., ALT=<cnv100>

[nchuang]

The post above said ignoring multi-allelic. So I am assuming it's biallelic.

[Pierre Lindenbaum]

your awk script is wrong for multi-allelic sites: like `REF=A ALT=T,G`

Answer 2

[Alex Reynolds]

Via BEDOPS convert2bed:

$ vcf2bed --snvs < foo.vcf | wc -l
$ vcf2bed --insertions < foo.vcf | wc -l
$ vcf2bed --deletions < foo.vcf | wc -l

Answer 3

[Charles Warden]

You can also use the summary from bcftools stats to answer this question.

http://samtools.github.io/bcftools/bcftools.html#stats

[chinhhq83]

Hi Charles,

I crossed your post about deepvariant caller and think you can advise me on a bioinformatic pipeline for analyzing genetic data. We generated DNA sequencing by using TrueSight Cancer 94, Illumina, and were not sure how we should analyze them, from aligment to variant calling to annotate to variant classification, and how filter variant. I would appreciate if you can help me. Thanks for your time and consideration!

Chinh

[brunobsouzaa]

This might help you!

[Charles Warden]

Hi - I think you meant this as a comment for my other answer.

Also, I think you might be referring to this blog post:

http://cdwscience.blogspot.com/2019/05/precisionfda-and-custom-scripts-for.html

My sense is that it helps to be able to visually inspect variation calls, and I am not sure I would be comfortable with a result if you couldn't.

I have some sense that you might need to take some extra caution with targeted gene panels, which I think is a slightly different blog post (even though it is relevant to the idea that I thought the default set of DeepVariant calls could benefit from some additional filtering, kind of like the unfiltered GATK calls):

http://cdwscience.blogspot.com/2019/06/general-comments-on-low-frequency.html

However, to be honest, I don't have specific advice for the TrueSight Cancer 94 gene panel.