Question: SNP call using bcftools
2
gravatar for blur
3.9 years ago by
blur90
European Union
blur90 wrote:

Hi,

 

I am trying to do SNP call using bcftools, I tried running this in samtools:

samtools mpileup -u -f CBS138.fasta 11_304.sorted.bam | bcftools view -bvcg - > 304.raw.bcf

The samtools part works fine and I get an outfile that looks OK, but the bcf parameters are wrong

(-bvcg doesn't work).

I looked this up online and everywhere I look this is the command line I see - 

can anyone tell me where I went wrong?

 

Thanks,

snp samtools bcftools • 11k views
ADD COMMENTlink modified 3.9 years ago by Evgeniia Golovina980 • written 3.9 years ago by blur90

Which version of bcftools are you using?

ADD REPLYlink written 3.9 years ago by iraun3.5k

version 1.2

ADD REPLYlink written 3.9 years ago by blur90
7
gravatar for Evgeniia Golovina
3.9 years ago by
New Zealand
Evgeniia Golovina980 wrote:

If you are interested only in SNPs:

samtools mpileup --skip-indels -d 250 -m 1 -E --BCF --output-tags DP,DV,DP4,SP -f <reference genome.fa> -o <output.bcf> <list of input bam files>

bcftools index  <output.bcf> <indexed.bcf>

bcftools call --skip-variants indels --multiallelic-caller --variants-only  -O v <output.bcf> -o <output.vcf>

 

--skip-indels (samtools) and --skip-variants indels (bcftools) - report only SNPs

-d 250 - maximum read depth to consider per position

-m 1 - minimum number of gapped reads for an INDEL candidate

-E - recalculate BAQ on the fly, ignore existing BQ tags

--BCF - output is bcf file

--output-tags DP,DV,DP4,SP:

DP (number of reads covering position), DV (number of high-quality variant reads), DP4 (number of forward reference, reverse reference, forward non-reference and reverse non-reference alleles used in variant calling) and SP (phred-scaled strand bias P-value) tags in the output file.

--multiallelic-caller call only multi-allelic variants

--variants-only - use if you are interested to report only potential variant sites and exclude monomorphic ones (sites without alternate alleles)

 

Ask, if you have any questions. Hope, it helps you :)

ADD COMMENTlink modified 3.9 years ago • written 3.9 years ago by Evgeniia Golovina980

BTW - there is a small error in the cmd line: -o needs to be before the output file name

bcftools call --skip-variants indels --multiallelic-caller --variants-only  -O v <indexed.bcf> 
-o <output.bcf>
ADD REPLYlink written 3.9 years ago by blur90

Hmm, okey. :)

ADD REPLYlink written 3.9 years ago by Evgeniia Golovina980

Is the "DV" indicator technically the number of nucleotide that support that snp?

ADD REPLYlink written 2.7 years ago by shanasabri40

Yes, technically, DV - is a number of high-quality variant reads that support each of the reported SNP.

ADD REPLYlink written 2.7 years ago by Evgeniia Golovina980

I've tried with the above command,but it shows;

Blockquote [warning] tag DP4 functional, but deprecated. Please switch to ADF and ADR in future. Could not parse tag "AF" in "DP,DV,DP4,SP,AF"

i also need the allele frequency in my vcf file. Thank you.

ADD REPLYlink written 3 months ago by ashaneev070
2
gravatar for Pierre Lindenbaum
3.9 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum117k wrote:

with 1.2 the workflow should involve "bcftools call" http://www.htslib.org/workflow/

samtools mpileup -go <study.bcf> -f <ref.fa> <sample1.bam> <sample2.bam> <sample3.bam>
bcftools call -vmO z -o <study.vcf.gz> <study.bcf>
ADD COMMENTlink modified 3.9 years ago • written 3.9 years ago by Pierre Lindenbaum117k

the bcf file is created by the -go parameters?

ADD REPLYlink written 3.9 years ago by blur90
0
gravatar for blur
3.9 years ago by
blur90
European Union
blur90 wrote:

would you be so kind as to explain what the cmd line means exactly? thank you

ADD COMMENTlink written 3.9 years ago by blur90

Thank you all for your help - 

I managed to do the first part (mpileup) 

but when I did index for the bcf file, I got an indexed bcf file, not a vcf file:

304.bcf.csi

No vcf file was created...

ADD REPLYlink written 3.9 years ago by blur90

Ahh, yes, thanks. I edited command line!

ADD REPLYlink written 3.9 years ago by Evgeniia Golovina980

thanks! this makes more sense :)

 

after the last command line - the file containing the SNPs in the <out.bcf>?

is it the same as the first bcf file?

ADD REPLYlink written 3.9 years ago by blur90

You are welcome! :)

No, they are not the same. Read more there --> A: What's the difference between mpileup output and bcftools call ?

ADD REPLYlink modified 3.9 years ago • written 3.9 years ago by Evgeniia Golovina980

thank you! you are a saint! :)

ADD REPLYlink written 3.9 years ago by blur90
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: 1883 users visited in the last hour