How Can I Count Snps In My Final Vcf Files

Question: How Can I Count Snps In My Final Vcf Files

2.5 years ago by

Hi all,

I have a VCF file that containing 50 samples, i want to count the number of SNPs. My organism is non-model, So it does not have the chromosome.

Now, How can i count the number of SNPs for all 50 samples with this VCF?

Best Regard

Mostafa

snp • 3.6k views

ADD COMMENT • link •

modified 4 weeks ago by ghada.alqubati101 • 0 • written 2.5 years ago by mostafarafiepour • 110

run bcftools stats on vcf. It would summarize the VCF with most of the details you need.

ADD REPLY • link modified 2.5 years ago • written 2.5 years ago by cpad0112 ♦ 15k

try this

bcftools query -f '%POS\n' file.vcf.gz | wc -l

ADD REPLY • link written 4 weeks ago by ghada.alqubati101 • 0

2.5 years ago by

finswimmer ♦ 14k

Germany

finswimmer ♦ 14k wrote:

Hello mostafarafiepour,

you've started with one question. In the meantime there are three :)

1. How to read a vcf file

This is a very basic question. So you need starting some literature:

Wikipedia
the official specs (already posted by Pierre Lindenbaum )
Have a look into the header of your file. Lot's of explanations are there

If you doesn't understand any of the explanations, don't worry to ask.

2. How to count the variants in a vcf (your original question)

We already worked on this.

3. Is the resulted number of (2) correct?

Well, that's quite hard to say without knowing anything about your genome. How large is it? Is there a high diversity between individuals? As we just have the total number of different variants in all of your samples, it might be better to get a per sample count. The output of bcftool stats (as suggested by cpad0112 ) might be useful or have a look at this thread, especially the answers by Pierre and me.

fin swimmer

ADD COMMENT • link modified 2.5 years ago • written 2.5 years ago by finswimmer ♦ 14k

2.5 years ago by

finswimmer ♦ 14k

Germany

finswimmer ♦ 14k wrote:

Hello,

the total number you get by counting the lines in the vcf excluding the header lines.

$ grep -v "^#" input.vcf|wc -l

fin swimmer

ADD COMMENT • link written 2.5 years ago by finswimmer ♦ 14k

Hi swimmer,

Many thanks for your reply, I've done it and gave me a number(20546654). But I do not know how correct this number is?

Code i use:

grep -v "^#" Final_VCF_50Sample.vcf|wc -l

ADD REPLY • link modified 2.5 years ago by finswimmer ♦ 14k • written 2.5 years ago by mostafarafiepour • 110

Sorry swimmer,

Please describe the details in the photo for me.

enter image description here

ADD REPLY • link written 2.5 years ago by mostafarafiepour • 110

Please describe the details in the photo for me.

how is it related to your original question ? Are you sure you're using the correct terms ?

Each line of the VCF is a VARIANT. A Variant can be a SNP or an INDEL or etc...

The intersection of the Variant and the Samples' names is a GENOTYPE.

ADD REPLY • link written 2.5 years ago by Pierre Lindenbaum ♦ 134k

Yes, maybe I did not ask the exact question.

I want to know what the meaning of any of the terms in the picture is?

for example:

CHROM

POS

REF

ALT

QUAL

INFO

GT:AD:DP:GQ:PL

AND ......

ADD REPLY • link modified 2.5 years ago • written 2.5 years ago by mostafarafiepour • 110

I want to know what the meaning of any of the terms in the picture is?

https://samtools.github.io/hts-specs/VCFv4.3.pdf

ADD REPLY • link written 2.5 years ago by Pierre Lindenbaum ♦ 134k

5 weeks ago by

dario.galanti92 • 0

dario.galanti92 • 0 wrote:

Here is a quick way to count biallelic SNPs in vcf.gz files (use "cat" instead of "zcat" for uncompressed vcf files):

zcat input.vcf.gz | awk '{if ($4~/^[ACGT]$/ && $5~/^[ACGT]$/){c++}} END {print c}'

ADD COMMENT • link modified 4 weeks ago • written 5 weeks ago by dario.galanti92 • 0

I'm sorry but this is wrong.

> echo "chr1 123 . A C,<INDEL>,G" | awk '{if ($4~/[ATCG]/ && $5~/[ATCG]/) {c++}} END {print c}'
1

better use:

bcftools view --no-header -G -m 2 -M 2 --types snps input.vcf.gz | wc -l

ADD REPLY • link written 5 weeks ago by Pierre Lindenbaum ♦ 134k

Yes, that was indeed wrong, apologies. I fixed it now.

ADD REPLY • link modified 4 weeks ago • written 4 weeks ago by dario.galanti92 • 0

5 weeks ago by

4galaxy77 • 280

United Kingdom

4galaxy77 • 280 wrote:

If all your variants in the vcf are SNPS, then a very quick way is to first index and then index again with the -n flag.

bcftools index data.vcf 
bcftools index -n data.vcf

ADD COMMENT • link modified 5 weeks ago • written 5 weeks ago by 4galaxy77 • 280

4 weeks ago by

ghada.alqubati101 • 0

ghada.alqubati101 • 0 wrote:

try this

bcftools query -f '%POS\n' file.vcf.gz | wc -l

ADD COMMENT • link modified 4 weeks ago • written 4 weeks ago by ghada.alqubati101 • 0

why do you output %POS ?

 bcftools query -f 'x' input.vcf | wc -c

ADD REPLY • link written 4 weeks ago by Pierre Lindenbaum ♦ 134k

Please log in to add an answer.

Similar posts • Search »