Biostar Beta. Not for public use.
How to count SNPs, InDels
3
Entering edit mode
5.8 years ago
chchl7 • 30
@chchl717706

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 • 5.0k views
ADD COMMENTlink
0
Entering edit mode

in a BAM or in a VCF ?

ADD REPLYlink
4
Entering edit mode
5.8 years ago
Vivek ♦ 2.3k
@Vivek5462

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
1
Entering edit mode

Alex's solution is faster and better.

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
1
Entering edit mode

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
0
Entering edit mode

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

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
2
Entering edit mode
3.7 years ago
@Alex Reynolds20

Via BEDOPS convert2bed:

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

Login before adding your answer.

Similar Posts
Loading Similar Posts
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.3