Question: How to count SNPs, InDels
3
gravatar for chchl7
3.9 years ago by
chchl730
United States
chchl730 wrote:

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!

sequencing snp next-gen genome • 4.7k views
ADD COMMENTlink modified 3.9 years ago by Alex Reynolds27k • written 3.9 years ago by chchl730

in a BAM or in a VCF ?

ADD REPLYlink written 3.9 years ago by Pierre Lindenbaum118k
4
gravatar for Vivek
3.9 years ago by
Vivek2.2k
Denmark
Vivek2.2k wrote:

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:

http://vcftools.sourceforge.net/htslib.html#stats

ADD COMMENTlink modified 3.9 years ago • written 3.9 years ago by Vivek2.2k
1

Alex's solution is faster and better.

ADD REPLYlink written 2.0 years ago by SmallChess480

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

ADD REPLYlink written 3.9 years ago by chchl730
1

insertions:

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

deletions:

awk '! /\#/' variants.vcf | awk '{if(length($5) > 1) print}' | wc -l
ADD REPLYlink written 3.7 years ago by nchuang190

wrong. if the variation is REF=A ALT=T,G or if the variation is symbolic, eg. ALT=<cnv100>

ADD REPLYlink written 3.7 years ago by Pierre Lindenbaum118k

the post above said ignoring multi-allelic. so i am assuming it's biallelic. 

ADD REPLYlink written 3.7 years ago by nchuang190

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

ADD REPLYlink written 3.9 years ago by Pierre Lindenbaum118k
2
gravatar for Alex Reynolds
3.9 years ago by
Alex Reynolds27k
Seattle, WA USA
Alex Reynolds27k wrote:

Via BEDOPS convert2bed:

$ vcf2bed --snvs < foo.vcf | wc -l
$ vcf2bed --insertions < foo.vcf | wc -l
$ vcf2bed --deletions < foo.vcf | wc -l
ADD COMMENTlink modified 22 months ago • written 3.9 years ago by Alex Reynolds27k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 800 users visited in the last hour