RAD-Seq

[westeros]

Hi all,

I am running a RAD-Seq analysis and here is how things are right now:

1) We ran a RAD-Seq with the parental individuals with a deep coverage (paired-end);
2) There is a reference genome available;
3) 100 F2 individuals with not-so-deep coverage, but still good (single-end).

What I managed to do using novoalign/samtools:
1) Mapping and SNP calling in the parental species using the reference genome;
2) Mapping and SNP calling in the 100 individuals.
3) The problem: only 4% of the sites overlap and I guess it is because of the lower coverage of the F2 samples.

What I would like to do, but haven't found a way yet: Take my selected SNPs from the parents and use just "zoom in" each individual to check how is that position and classify its genotype.

I have decent programming skills, but was wondering if:
1) Is this correct to annotate my F2 using this strategy ?
2) Is there some tool/script available to do that or do I have to code it ?

Thanks in advance.

[westeros]

Any help here ? Thanks.

[sarvidsson]

Did you have a look at Stacks (http://creskolab.uoregon.edu/stacks/)?

samtools (I guess you use samtools mpileup) sounds like a bad plan for genotyping such data, especially with that overlap (you'll get many false het calls). So you may want to try to use freebayes (I've had some luck using that with similar data), possibly with some cleanup afterwards (low coverage sites), unless you're successful with Stacks (might fail if the overlap is that bad).

How does the partial overlap look like, e.g. how many variable sites are covered in lets say 50, 67 and 85 % of the samples? You'll get very, very few sites covering 100 % of the samples, but that is not a big problem, you can impute much of the missing data since you have a crossing population and reasonable distances between well-genotyped sites.

[SNPsaurus]

You should first limit the SNPs found in the parents to sequence space assayed by the SE reads in the F2s to get a sense how much of the non-overlap is due to coverage type issues and how much is just impossible to get overlap (e.g. a SNP in parent 200 bp from cut site).

Is this RAD and RAD-PE or paired-end ddRAD? The issue is the same either way.

What I would do is identify polymorphic loci in the parents that have SNPs in the first 100 bp (so therefore can be identified in the F2 sequencing if using 100 bp reads). Then create a pseudo-reference fasta of just those polymorphic loci and align the F2 reads to that reference. Freebayes as sarviddson mentioned uses population data in calling SNPs so could rescue low-coverage if that is the problem (what kind of coverage are you talking about?). Or just go with a heuristic approach and just look, for example, for 3 reads and at least 10% of reads containing the allele of the parent to call it.

One issue you may be having is if you haven't removed repetitive reads from the SNP calling. High-abundance RAD reads will have many sequence-error variants that can look like alleles, so the high-abundance reads and the derived sequences should be removed.

If you have a read depth of 10 or so in each F2 then I would be really surprised if a SNP in the parent in the right position to be assayed in the F2 would not be called. Are you creating a vcf? You should enable allele depth reporting and examine the polymorphic sites.