Question: SNP,INDEL counting per chromosomes in vcf
1
gravatar for Siavash Salek Ardestani
3.0 years ago by

Hi all, how can I counting SNPs and INDELs per each chromosome in a vcf file?

snp next-gen vcf • 2.4k views
ADD COMMENTlink modified 3.0 years ago by Alex Reynolds31k • written 3.0 years ago by Siavash Salek Ardestani20
2
gravatar for Alex Reynolds
3.0 years ago by
Alex Reynolds31k
Seattle, WA USA
Alex Reynolds31k wrote:

Using BEDOPS vcf2bed and standard Unix tools like wc to count lines:

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

Cf. http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html

The above does counts over all chromosomes. You need to do an extra step to count per chromosome.

To perform this count exercise quickly per-chromosome, per-variant class:

$ vcf2bed --snvs < variants.vcf > variants.snvs.bed
$ for chr in `bedextract --list-chr variants.snvs.bed`; do echo $chr; bedextract $chr variants.snvs.bed | wc -l; done

The bedextract tool is part of BEDOPS. If you install BEDOPS to get vcf2bed, you get bedextract, as well.

Repeat this procedure for --insertions and --deletions, replacing snvs accordingly, to count those classes.

ADD COMMENTlink modified 3.0 years ago • written 3.0 years ago by Alex Reynolds31k

it is correct. I humbly thank you,

ADD REPLYlink modified 3.0 years ago • written 3.0 years ago by Siavash Salek Ardestani20

I run these commands and I have number of SNP per chromosomes furthermore I have many 'Un' (I think this is unknown chromosomes SNP). how can I split or cut mt,x and unknown chromosomes from vcf?

ADD REPLYlink written 3.0 years ago by Siavash Salek Ardestani20

Use grep on the output of bedextract --list-chr. You can pass grep a file containing chromosome names you want to include (or, conversely, exclude) via grep -f.

$ for chr in `bedextract --list-chr variants.snvs.bed | grep -f wanted-chromosomes.txt -`; do ... ; done

Run man grep to see a full list of options.

ADD REPLYlink written 3.0 years ago by Alex Reynolds31k

many thanks but I want cut mt, unknown and X chromosomes SNPs from a vcf file.

I means, my prior vcf have many SNPs that these SNPs are in mitochondrial DNA, X chromosome and unknown chromosomes. I want to have a vcf file without mitochondrial , X chromosome and unknown chromosomes SNPs information.

ADD REPLYlink written 3.0 years ago by Siavash Salek Ardestani20

You could use grep in the same way. It may help to review some examples of how it works so you can understand Unix streams and how grep can be used here.

ADD REPLYlink written 3.0 years ago by Alex Reynolds31k

Hi, I also tried using this script but gave me only the number of snvs for each chromosome and did not count the insertions and deletions??

In fact, the outputs related to insertions and deletions are empty and there is no data in it?

ADD REPLYlink written 2.4 years ago by mostafarafiepour100

Are you sure you have insertions/deletions?

ADD REPLYlink written 2.4 years ago by _r_am32k

yes, Because when I use this command, i find 182456 deletions and 0 insertions.

insertions:

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

insertions:

awk '! /\#/' variants.vcf | awk '{if(length($4) > 1 ) print}' | wc -l
ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by mostafarafiepour100
0
gravatar for Pierre Lindenbaum
3.0 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum133k wrote:

using bioalcidaejdk: http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

$ java -jar dist/bioalcidaejdk.jar -e 'stream().map(V->V.getType().name()+" "+V.getContig()).collect(Collectors.groupingBy(Function.identity(), Collectors.counting())).forEach((K,V)->{println(K+" : "+V);});'  input.vcf

SNP 1 : 44
INDEL 1 : 6
ADD COMMENTlink written 3.0 years ago by Pierre Lindenbaum133k

thanks but Did not respond

ADD REPLYlink written 3.0 years ago by Siavash Salek Ardestani20

Hi, i have been trying to do your scripts. But counted the only number of SNPs and did not count the InDels. ? Please see the photo in the attachment.

enter image description here

ADD REPLYlink modified 2.4 years ago • written 2.4 years ago by mostafarafiepour100

Read the part that begins with [WARN]

ADD REPLYlink written 2.4 years ago by _r_am32k

the warning is meaningless here.

ADD REPLYlink written 2.4 years ago by Pierre Lindenbaum133k

Ah I see. It's probably my second guess then - no InDels to pick.

ADD REPLYlink written 2.4 years ago by _r_am32k
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: 2104 users visited in the last hour
_