Recently, I'm using samtools mpileup to call SNPs using low coverage (average 1x) resequencing data of 600+ populations. But I encountered a few problems.
My command line is as follows:
samtools mpileup -E -C50 -Q20 -q20 -DSuf ref.fasta -r chr4:21000001-24000000 -b bam_list | bcftools view -Nbcvg - > partE-chr4: 21000001-24000000.bcf
The individuals sequenced are of inbred lines. So we want to get homoSNPs. HeteroSNPs shold be wrong calling.
1) The SNPs printed by samtools mpileup seemed weird:
SNP1: 0/1:26,3,0:1:0:10 SNP2: 0/1:26,3,0:1:0:9
Both SNP1 and SNP2 have PLs with (26,3,0), but samtools- bcftools printed them as hetero.
lh3 have give a answer at the post: Why does samtools/bcftools give incorrect genotypes and innacurate quality scores?
but finally which genotype should I give? First, I think I should compare the three values of PLs, and choose the maximum value as the final genotype. If two of the values were very close, how to give the genotype? I'm not sure of this. Can someone give me some suggestions and some explanations?
2) For large population, how to deal with QUAL columns (the sixth column)? How to choose the cutoff?
This post was the same as in seqanswer (http://seqanswers.com/forums/showthread.php?t=25354).