I'm using samtools multisample SNP calling to call SNPs from different samples. The command I'm using is,
samtools mpileup -uf ref.fa *.bam | bcftools view -vcg - > out.vcf
Then, I'm using vcf-to-tab tool to convert vcf to tab file.
The problem is, No calls (missed in sequencing) are called as homozygous reference calls. How to differentiate No calls from homozygous reference calls.
I appreciate any sort of help.
Edit: My vcf file looks like this:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samp1.bam samp2.bam 000000F 634700 . T C 3.55 . DP=1;AF1=1;AC1=4;DP4=0,0,0,1;MQ=60;FQ=-27.4 GT:PL:GQ 0/1:31,3,0:5 0/1:0,0,0:3
•
written
by
second_exon
1 answer
No calls (missed in sequencing) are called as homozygous reference calls. How to differentiate No calls from homozygous reference calls.
According to the bcftools documentation:
-g, --gvcf INT output also gVCF blocks of homozygous REF calls. The parameter INT is the minimum per-sample depth required to include a site in the non-variant block.
Looking at your command line, I notice that you haven't specified a minimum per-sample depth for the -g parameter. Try setting that parameter to a sane, non-zero value and you should be able to differentiate between no-call and homozygous REF.
Log in to answer this question.
-
VCF file REF column only contains N
written by genie66 •
-
VCF file REF Column outputs all 'N'
written by mhasa006
-
Hom vs Het Variant calls
written by hellbio
-
SNPs only output from bcftools
written by biogirl
-
Why Is The Variant File Produced By Samtools So Large? It Is About ~34G.
written by Jordan
-
Extracting Specific Reads From Vcf File
written by Bnfoguy •
-
Strand Information On Vcf
written by Jirapong •
-
Samtools Mpileup Vs Pileup: Multiple Samples And Gt=0/0
written by Pierre Lindenbaum
-
Mpileup By Samtools And How To Get The Counts Per Allele Per Position In The Vcf
written by Deniz •
-
How To Filter The Snps/Indel Results From Mpileup?
written by Haiping •
Generally, you'd do that by depth. If the depth for that sample is below your criteria for safely considering something homozygous-ref at a specific position, then consider it a no-call.
Thanks, Brain. That's a very good idea. But, raw read depth is not printed for individual samples, how to force bcftools to print raw read depth?
Ah, that's kind of unfortunate. I wonder if there's an option for that somewhere? I'm not all that familiar with bcftools, but for some reason depth and allele fraction do not seem to be default fields for each sample. Maybe you're supposed to guess from the PL field ("phred-scaled genotype likelihoods"); since they are all zero for the second sample, the implication is that the depth is zero, even though it was called 0/1 (het). Well, that doesn't make much sense!
One alternative is to use a different tool to calculate coverage at that location for that sample. Bedtools probably allows that, and BBMap's pileup.sh does as well. Not a convenient option, though, for sure.