Snp calling between samples

[Rainbird]

There is a reference genome for subspecies A available. We did a resequencing for subspecies B and subspecies C using Hiseq 2000 and we'd like to know the SNP diversity(difference in allele frequency) between subspecies B and subspecies C.

We know there are substantial differences among subspecies A, B and C. So, what is the best way to find the SNP diversity between B and C? Should we do snp calling for subspecies B or C separately based on subspecies A reference genome (I used samtools) and then merge the results ? It seems to me we need a more efficient way to do the job but I don't know much (cortex ?).

Thanks in advance.

[swbarnes2]

Align both samples to your best reference, then use samtools mpileup on both .bams together.

[Rainbird]

Thanks swbarnes2

Could you explain a little more why "use samtools mpileup on both .bams together" will work? In that case, we still need the reference genome from subspecies A, right?

Another thing is that: if there are more than 1 non-reference allels reported, the samtools only gives out the depth of the 1st non-reference allel (as listed in DP4). Also, although the 1/1 indicates homozygous alternate, I don't understand the meaning of the PL value which is "131,59,26,91,0,85" (as shown below). How can we get the depths and other information for the 2nd alternate ?

chr2 213263 . A C,T 72 . DP=14;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,9,4;MQ=56;FQ=-60 GT:PL:GQ 1/1:131,59,26,91,0,85:63

[Rainbird]

Anyone can help?

[Khen]

Calling species B and C against the reference together just saves space. And yes, you will still need to use the reference genome. The output is slightly different however, so what you will get is an extra GT field:

chr2	213263	.	A	C,T	72	.	DP=14;VDB=0.0355;AF1=1;AC1=2;DP4=0,0,9,4;MQ=56;FQ=-60	GT:<Genotype of A>	GT:<Genotype of B>

I'm pretty sure that the PL field is reporting the quality score of all of the allelic possibilities, which is why you see six of them. You will have to consult the documentation for how to get multiple sample depth information.
Also, I find that the Broad Institute does a much better job documentation than does sourceforge or 1000genomes.org. Since samtools and gatk both use VCF as the standard output, you might want to start with the GATK documentation if not just switch to GATK altogether.

Hope this helps.

[Rainbird]

Thanks Khen.
If I undertand correctly, samtools is designed for diploid genome. If there are 2 alleles in your sample other than the allele in the reference genome (for example, the reference genome has a T, and you have a A and a G in your sample), samtool might not work well.

Is there any tool specifically designed to find allele frequency in your own samples regardless what is in the reference genome?