Question: How to count SNPs, InDels
3
gravatar for chchl7
4.6 years ago by
chchl740
United States
chchl740 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 • 5.4k views
ADD COMMENTlink modified 4.6 years ago by Alex Reynolds29k • written 4.6 years ago by chchl740

in a BAM or in a VCF ?

ADD REPLYlink written 4.6 years ago by Pierre Lindenbaum124k
4
gravatar for Vivek
4.6 years ago by
Vivek2.3k
Denmark
Vivek2.3k 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 4.6 years ago • written 4.6 years ago by Vivek2.3k
1

Alex's solution is faster and better.

ADD REPLYlink written 2.7 years ago by SmallChess500

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 4.6 years ago by chchl740
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 4.3 years ago by nchuang200

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

ADD REPLYlink written 4.3 years ago by Pierre Lindenbaum124k

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

ADD REPLYlink written 4.3 years ago by nchuang200

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

ADD REPLYlink written 4.6 years ago by Pierre Lindenbaum124k
2
gravatar for Alex Reynolds
4.6 years ago by
Alex Reynolds29k
Seattle, WA USA
Alex Reynolds29k 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 2.5 years ago • written 4.6 years ago by Alex Reynolds29k
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: 1425 users visited in the last hour