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,
2 answers
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> </code><code>-o </code><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 :)
•
written
by
Evgeniia Golovina
Hmm, okey. :)
Is the "DV" indicator technically the number of nucleotide that support that snp?
Yes, technically, DV - is a number of high-quality variant reads that support each of the reported SNP.
I've tried with the above command,but it shows;
[warning] tag DP4 functional, but deprecated. Please switch to
ADFandADRin future. Could not parse tag "AF" in "DP,DV,DP4,SP,AF"
i also need the allele frequency in my vcf file. Thank you.
•
written
by
ashaneev07 •
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>
Log in to answer this question.
-
Failed to open var.raw.bcf: unknown file type
written by Ric
-
Failed to open var.raw.bcf: unknown file type
written by Ric
-
Error in bcftools view -bvcg command
written by hafiz.talhamalik
-
samtools mpileup to bcftools mpileup
written by Ric
-
Empty VCF file with bcftools call
written by AP
-
How can I filter SNP for 30 samples in vcf files?
written by sun.nation
-
Reducing sequencing errors in .vcf file
written by lcc1844 •
-
bcftools: invalid option -- b was showed error
written by bigbang •
-
Snp Calling With Bcftools Using Pileup File Format Instead Of Bcf
written by Justin
-
Homozygous Ref Calls And No Calls
written by win
Which version of bcftools are you using?
version 1.2
would you be so kind as to explain what the cmd line means exactly? thank you
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...Ahh, yes, thanks. I edited command line!
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?
You are welcome! :)
No, they are not the same. Read more here --> What's the difference between mpileup output and bcftools call ?
thank you! you are a saint! :)