This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to count SNPs, InDels

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!

next-gen sequencing snp genome
• 14,635 views • link
updated by Ram Ram • written by chchl7 • chchl7

3 answers

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.

• 1 views • link
updated by Ram Ram • written by Vivek Vivek

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

• 1 views • link
written by chchl7 • chchl7

insertions:

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

deletions:

awk '! /\#/' variants.vcf | awk '{if(length($5) > 1) print}' | wc -l
• 0 views • link
updated by Ram Ram • written by nchuang nchuang

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

• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

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

• 0 views • link
updated by Ram Ram • written by nchuang nchuang

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

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Alex's solution is faster and better.

• 0 views • link
written by scchess scchess

Via BEDOPS convert2bed:

$ vcf2bed --snvs < foo.vcf | wc -l
$ vcf2bed --insertions < foo.vcf | wc -l
$ vcf2bed --deletions < foo.vcf | wc -l
• 1 views • link
updated by Ram Ram • written by Alex Reynolds Alex Reynolds

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

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

• 0 views • link
written by Charles Warden Charles Warden

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

• 0 views • link
updated by Ram Ram • written by chinhhq83 • chinhhq83

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.

• 0 views • link
written by Charles Warden Charles Warden

Log in to answer this question.

More posts like this