I have a question on the impact on quality of SNP calls using samtools mpileup: calling each sample individually vs calling all samples altogether
I basically follow the samtools URL:
http://samtools.sourceforge.net/mpileup.shtml and
http://samtools.sourceforge.net/samtools.shtml (except using -D 2000 or -d 10 for filtering)
the commands I used basically like below:
Call all samples together (each bam file is from exome-seq of one sample )
samtools mpileup -ugf ref.fa aln1.bam aln2.bam aln3.bam ... *.bam| bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -d 10 > var.flt.vcf
Call one samples a time (each bam file is from one sample)
for loop (for each bam file x.bam)
samtools mpileup -ugf ref.fa x.bam| bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -d 10 > var.flt.vcf
end of loop
We observed that if we use samtools to call each sample (bam file) individually one at a time or call all samples (multiple bam files) together, then we seem get a better SNP calls if we call individually than if we call them altogether (of course we get one vcf file for each sample/bam file if call individually, whereas just one single vcf file for all SNPs of all samples for calling them altogether. We do lose homo reference SNPs in vcf file if calling individually, since for single sample, we only see hetero or homo variant, no homo reference unless we call multiple sample altogether, then some samples would have homo reference calls for sites that other samples were called as hetero or homo variants)
Does any one have similar observation or any comments or insight for such observation (get better SNP calls if we call individually than if call together)? Since I know most people would call all bam files altogether and would not think that calling each sample individually would get a better quality SNPs (using samtools mpileup SNP call). Any interpretation from algorithm side of samtools? We all thought it shall not matter, but we do have bench mark data to show that is true at least in our hand. We use -d 10 or -D 2000 in the filtering step for both cases.
I basically follow the samtools URL:
http://samtools.sourceforge.net/mpileup.shtml and
http://samtools.sourceforge.net/samtools.shtml (except using -D 2000 or -d 10 for filtering)
the commands I used basically like below:
Call all samples together (each bam file is from exome-seq of one sample )
samtools mpileup -ugf ref.fa aln1.bam aln2.bam aln3.bam ... *.bam| bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -d 10 > var.flt.vcf
Call one samples a time (each bam file is from one sample)
for loop (for each bam file x.bam)
samtools mpileup -ugf ref.fa x.bam| bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -d 10 > var.flt.vcf
end of loop
We observed that if we use samtools to call each sample (bam file) individually one at a time or call all samples (multiple bam files) together, then we seem get a better SNP calls if we call individually than if we call them altogether (of course we get one vcf file for each sample/bam file if call individually, whereas just one single vcf file for all SNPs of all samples for calling them altogether. We do lose homo reference SNPs in vcf file if calling individually, since for single sample, we only see hetero or homo variant, no homo reference unless we call multiple sample altogether, then some samples would have homo reference calls for sites that other samples were called as hetero or homo variants)
Does any one have similar observation or any comments or insight for such observation (get better SNP calls if we call individually than if call together)? Since I know most people would call all bam files altogether and would not think that calling each sample individually would get a better quality SNPs (using samtools mpileup SNP call). Any interpretation from algorithm side of samtools? We all thought it shall not matter, but we do have bench mark data to show that is true at least in our hand. We use -d 10 or -D 2000 in the filtering step for both cases.
Comment