DiscoSnp-RAD: de novo detection of small variants for RAD-Seq population genomics

Abstract

Restriction site Associated DNA Sequencing (RAD-Seq) is a technique characterized by the sequencing of specific loci along the genome that is widely employed in the field of evolutionary biology since it allows to exploit variants (mainly Single Nucleotide Polymorphism-SNPs) information from entire populations at a reduced cost. Common RAD dedicated tools, such as STACKS or IPyRAD, are based on all-vs-all read alignments, which require consequent time and computing resources. We present an original method, DiscoSnp-RAD, that avoids this pitfall since variants are detected by exploiting specific parts of the assembly graph built from the reads, hence preventing all-vs-all read alignments. We tested the implementation on simulated datasets of increasing size, up to 1,000 samples, and on real RAD-Seq data from 259 specimens of Chiastocheta flies, morphologically assigned to seven species. All individuals were successfully assigned to their species using both STRUCTURE and Maximum Likelihood phylogenetic reconstruction. Moreover, identified variants succeeded to reveal a within-species genetic structure linked to the geographic distribution. Furthermore, our results show that DiscoSnp-RAD is significantly faster than state-of-the-art tools. The overall results show that DiscoSnp-RAD is suitable to identify variants from RAD-Seq data, it does not require time-consuming parameterization steps and it stands out from other tools due to its completely different principle, making it substantially faster, in particular on large datasets.

Keywords: Deletions; Insertions; RAD-seq; Reference-free; SNPs; Variants.

Figure 1. Examples of bubbles detected by SNPs in a toy de Bruijn graph, with k = 4.

In (A) the bubble is complete: this corresponds to a bubble detected by DiscoSnp++. In (B), the bubble is symmetrically truncated: it is composed of a branching node (“ACTG”) whose two successors lead to two distinct paths that both have the same length and such that their last two nodes have no successor. The graph (C) shows an example of two bubbles from the same locus. The leftmost bubble contains two symmetrically branching crossroads.

Figure 2. Recall (A), precision (B), time (C) and space (D) evolution on simulated data with different sampling sizes.

For the sampling of 100 samples, five parameter sets were tested for IPyRAD and STACKS (see “Material and Methods” for details).

Figure 3. Recall and precision on simulated data of 100 samples using DiscoSnp-RAD with respect to (A).

k-mer sizes, (B) maximal number of authorized SNP per bubble, (C) maximal number of authorized substitutions while mapping reads on predicted variants sequences, and (D) maximal number of symmetrically branching crossroads. Dashed vertical line represents on each plot the chosen default value.

Figure 4. (A) RAxML phylogeny realized on all variants predicted by DiscoSnp-RAD.

Bootstrap node supports > 80 are shown denoted by gray dots, bootstrap node supports > 90 are shown denoted by black dots. (B) STRUCTURE results obtained with SNP only and all variants on the seven Chiastocheta species. (C) Plot of the two first PC from a multivariate analysis on C. lophota samples and (D) their geographic distribution (figure made with Natural Earth Contributors (2020)).