Question: I'm currently trying to use the GATK pipeline for SNP calling to call SNPs. At the step where you use GATK's 'BaseRecalibrator' tool, you have to provide a set of known SNPs. Unfortunately, for my organisms, these are not annotated.

However, these organisms do have a reference genome sequenced. In a nature paper from Prado et al, they created a...

> "set of training SNP positions that are independent of Illumina short read data based on capillary sequence traces obtained for each species from the NCBI trace archive based on the whole genome shotgun (WGS) reads produced for the species reference assemblies. We mapped capillary reads using ssaha2 and called SNPs using the neighborhood quality score criteria implemented in ssahaSNP."
(Prado-Martinez et al Nature 2013, supplementary information)

So looking at the NCBI [trace archive][1], it seems like there is a pretty nice ftp site (and even a java applet that you can use to download things). However, I'm a bit confused on how this will help me call SNPs. On the ssahaSNP page, they say that:

> "ssahaSNP is a polymorphism detection tool. It detects homozygous SNPs and indels by aligning shotgun reads to the finished genome sequence."
http://www.sanger.ac.uk/science/tools/ssahasnp

How can you call 'homozygous' SNPs from the finished genome sequence using shotgun reads? Shouldn't the finished genome be the consensus of those reads, and therefore homozygous SNPs would just be the base called in the reference?

Any ideas?

Thanks!

[1]: ftp://ftp-private.ncbi.nlm.nih.gov/pub/TraceDB/

Answer: BaseRecalibrator looks at ALL sequenced and aligned bases for mismatches to the reference genome.

It is assumed that these mismatches are as a result of sequencing error - unless, of course, they are SNPs. But since the whole point of using GATK is to *find* SNPs, obviously this part is not the most important part of the algorithm. Fortunately, if there is a SNP in the genome, it isn't going to effect BQSR's model *that* much if we erroneously include it when doing BQSR. This is because compared to the hundreds of thousands if not millions of single-base-sequencing-errors, single-base-SNPs are:

1. no where near as frequent,
2. should not significantly effect the metrics BQSR looks for, like bias in read position, etc.

So my advice is just run BQSR without a SNP list. You should be fine. If you have the time/patience to do your due dilligence, use a very simple SNP caller first - one that only calls SNPs if there are at least 3 or more non-duplicate reads supporting it - and use that output when running BQSR. It is highly unlikely that a simple sequencing error will be the cause of such called homo/hetero SNPs, which will improve BQSR's chances of creating an accurate model for your sequencer. Obviously this has problems too - but its a good tradeoff I think.

Regarding your second question about how can you have a homozygous SNP and a reference genome - the answer is really simple. An individual contains two copies of every chromosome, and the sequencing may be from multiple individuals. Thus, after sequencing just 2 individuals, you could have up to 4 different alleles for that 1 position. The reference sequence is decided to be the most common base at that position (in a species sense, not a read-depth sense), and occasionally for poorly sequenced genomes there are ties and someone has to make an arbitrary decision. So in your case, I imagine the shotgun sequencing comes from another individual - that's all.