Restriction site-associated DNA sequencing (RAD-seq) has become a powerful and widely used tool in molecular ecology studies as it allows to cost-effectively recover thousands of polymorphic sites across individuals of non-model organisms. However, its successful implementation in population genetics relies on correct data processing that would minimize potential loci-assembly biases and consequent genotyping error rates. RAD-seq data processing when no reference genome is available involves the assembly of hundreds of thousands high-throughput sequencing reads into orthologous loci, for which various key parameter values need to be selected by the researcher. Previous studies exploring the effect of these parameter values found or assumed that a larger number of recovered polymorphic loci is associated with a better assembly. Here, using three RAD-seq datasets from different species, we explore the effect of read filtering, loci assembly and polymorphic site selection on number of markers obtained and genetic differentiation inferred using the Stacks software. We find (i) that recovery of higher numbers of polymorphic loci is not necessarily associated with higher genetic differentiation, (ii) that the presence of PCR duplicates, selected loci assembly parameters and selected SNP filtering parameters affect the number of recovered polymorphic loci and degree of genetic differentiation, and (iii) that this effect is different in each dataset, meaning that defining a systematic universal protocol for RAD-seq data analysis may lead to missing relevant information about population differentiation.

I would like to suggest you to follow :Article Genomic signatures of environmental selection despite near-p...

Article Amplification Biases and Consistent Recovery of Loci in a Do...

Hi Ann,

I don't understand what do you mean by 'substantial amount of loci that have >10 variant sites'? In your VCF every variant site is at a particular locus with its position assigned according to the reference (Chromosome and position). If you could elaborate more on this that will be helpful.

If the problem is that you have many variable site/loci witing 250bp, then it is normal but for your analysis you should only choose one SNP per read. The reason for not taking very closely occurring SNPs is that they might be linked.

But this will depend upon what analysis you are going to do.

To do this, you can use --thin function of VCFtools and can specify what should be the minimum distance between two SNPs.

I will be happy to discuss it more and hope it helps. Thanks.

Hi Abhinav Tyagi !

Thanks a lot for trying to help.

So I try to be more precise: Yes, all variant sites are assigned a position (CHROM and POS), however, I expect my tags to be between 200 - 250 bp long. In one tag (=CHROM is the same, just the number for POS changes) I have many variants much more than 10. Clearly some of these variants or these loci are not reliable. Thank you for the hint with the workaround. It just becomes difficult as soon as I wish to use haplotype-based analysis.

Also, I spend weeks to find the best filtering strategy to get the most out of my data. What I don't understand is: how is it possible, that I do all the quality filtering to remove unreliable data and I still end up with a high number of loci with >10.

Also then, how do I deal with that in a reasonable and publication-friendly way?

Thanks again for your input.

best, ac

Thank you for the clarification Ann-Christin Honnen .

It is normal to get more then one variant/SNP in one RAD tag/ or read.

Whether to keep all of them or not, depends on the kind of downstream analysis you will be doing (If you think linked variable sites can effect your results and inferences, you remove them).

Also, if you are not satisfied and have doubts about the data, you can always go back and redo the analysis (Also I suggest using a different variant caller, just to check). Following is the basic work flow I follow to analyze my data:

1. FastQC of raw fastq files

2. Trimming the adapters and mapping the raw reads to the reference

3. After obtaining the bam files sort and index them

4. Perform indel realignment

5. Variant calling using indel realigned bam files

6. Filter variants based on mapping quality (minQ), Genotype quality (minGQ), Depth, missing data for individuals and sites, minor allele frequency or count (whatever is suitable), thinning for variants within 500bp window, along with the filters you have already mentioned in the question.

If required I will be happy to provide the software details. What variant caller you are using?

I also find the following filtering tutorial helpful, see if it is of any help.

https://www.ddocent.com/filtering/

Hope this is of some help. Thanks

Dear Abhinav Tyagi ,

thanks a lot for your reply. I have been quite busy so I initially missed your reply. Essentially I did all of the proposed step. Except the indel-realignment. I will check that. And thanks also for the link to the ddocent page.

For variant calling I used stacks and BCFTools. To keep stuff streamlined with my collaborators we decided to use stacks.

Thanks a lot again, this is really helpful for me.

cheers, ac