The evaluation of different combinations of enzyme set, aligner and caller in GBS sequencing of soybean

The evaluation of different combinations of enzyme set, aligner and caller in GBS sequencing of soybean | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article The evaluation of different combinations of enzyme set, aligner and caller in GBS sequencing of soybean Aleksei Zamalutdinov, Stepan Boldyrev, Cécile Ben, Laurent Gentzbittel This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-5821852/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 06 Aug, 2025 Read the published version in Plant Methods → Version 1 posted 9 You are reading this latest preprint version Abstract Background Genotype-by-sequencing (GBS) is a cost-effective method for large-scale genotyping, widely used across various species, particularly those with large genomes. A critical aspect of GBS lies in the selection of restriction enzymes for genome digestion and the optimization of data analysis pipelines. However, few studies have comprehensively examined the combined effects of enzyme choice and pipeline configuration. Results In this study, we created GBS libraries using three commonly used restriction enzyme combinations ( HindIII - NlaIII , PstI - MspI , and ApeKI ) and assessed multiple SNP-calling pipelines in 15 soybean varieties. We tested four aligners (BWA-MEM, Bowtie2, BBMap, and Strobealign) and seven SNP callers (Bcftools, Stacks, DeepVariant, FreeBayes, VarScan, BBCallVariants, and GATK). Our finding reveal that enzyme choice significantly influences the number of identified SNP, gene localization preferences, and accuracy. Furthermore, the performance of SNP callers varied markedly in terms of SNP count, precision, recall, and false discovery rate (FDR). DeepVariant exhibited the highest accuracy, with 76.0% of its SNPs intersecting with whole-genome sequencing (WGS)-derived SNPs and an FDR of 0.0095, compared to FreeBayes, which had 47.8% intersection and an FDR of 0.6321. Conclusions Our findings underscore the importance of optimizing both enzyme selection for sequencing libraries and data analysis pipelines to ensure robust and reproducible results. This study provides a general framework for designing large-scale genotyping experiments aimed to specific quality and quantity requirements in various plant species. Genotype-by-sequencing (GBS) ddRAD SNP calling soybean comparison Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Background Genotyping is a crucial step in plant and animal breeding, as well as ecological population studies. It can be carried out in a number of ways, including microarray chips, WGS and reduced representation sequencing. Microarrays represent a large-scale and cost-effective approach that is widely used in human genetics and animal breeding. However, array design is a costly process [ 1 ], it works with a fixed set of polymorphisms, cannot identify new genetic variation in diverse populations and is susceptible to structural variation [ 2 ]. Furthermore, microarray chips are proprietary tools, which places researchers in a position of dependency with regard to the manufacturer and supplier, as well as their sales policies. WGS is the most informative approach in sequencing as it provides the genome-wide data for organism. WGS allows not only to retrieve information about polymorphisms in the entire genome but also to benefit from available pangenome assemblies and to efficiently study dispensable fractions of the genome [ 3 ]. However, it requires high coverages (up to 30x) to obtain accurate results, which can be prohibitively costly for large eukaryotic genomes or in case of genotyping large collections of accessions. As an alternative, several reduced representation sequencing (RRS) approaches have been developed to reduce the amount of required reads by enriching, separating, or eliminating a portion of the genome prior to sequencing. The majority of these approaches are costly (such as exome sequencing) and require the use of specialized commercial kits and high-quality annotation of the reference genome [ 4 ]. An alternative approach is to reduce genome complexity through the use of restriction enzymes. Restriction-site associated (RAD) [ 5 ] or Genotype-by-sequencing (GBS) [ 6 ] families of methods employ the use of one or more restriction enzymes for the digestion of genome DNA, followed by the adapter ligation. This approach is easy to implement due to the availability of published protocols and restriction enzymes at a low cost. A key question for the GBS method is the selection of an appropriate enzyme set. This is a complex problem with different parameters to consider such as differences in restriction site frequencies in different genomes, the rate of fragment occurrence in repetitive regions and restriction site methylation. It has been demonstrated that the use of methylation-sensitive enzymes allows for a better reduction of genome complexity and better coverage of gene-containing regions [ 6 , 7 ]. Nevertheless, it is worth noting that methylation-insensitive enzyme pairs are also used in GBS studies such as HindIII - NlaIII [ 8 ]. An important area for potential improvement is the bioinformatics pipeline for data analysis. Variations in read quality and genome characteristics can influence the approaches employed for analysis. There are numerous tools available for read filtering, mapping (Bushnell, 2014; Langmead & Salzberg, 2012; Li, 2013) and polymorphism calling [ 9 – 15 ]. Noteworthy, it has been demonstrated that the set of obtained SNPs is highly dependent on the caller that is employed [ 16 ]. In this study, three libraries of 15 Glycine max varieties were prepared using three different popular enzyme combinations: HindIII - NlaIII , PstI - MspI and ApeKI , and also subjected to 25x WGS. The libraries were tested with multiple pipelines using three mapping tools (BWA, Bowtie2 and BBMap) and seven SNP callers. The quality of the obtained SNPs was evaluated in different tool combinations by comparison with each other and with respect to WGS-based SNPs dataset identified within the 15 varieties, and using the most recent genome annotation. We show that the choice of enzyme combination significantly affects the number and localization of the SNPs obtained. The influence of aligner is not significant in case of total SNPs numbers and intersection SNPs set detected using WGS data. It is also of minor significance in genotyping accuracy. Callers differ in all evaluated criteria and DeepVariant showed the highest genotyping accuracy and lowest false positive rate. Methods Genotyping We selected a total of 15 soybean varieties of different world regions and with realized economic potential. All accessions were genotyped using GBS and WGS approaches (Additional file 1). GBS libraries were prepared for three enzyme combination of restriction enzymes according to established protocols ( ApeKI [ 6 ], HindIII - NlaIII [ 17 ], PstI - MspI [ 18 ]). The resulting libraries were then sequenced using the Illumina HiSeq4000 in a 2*150bp paired-end mode with 18.2, 31.5 and 19.5 million of reads per library, respectively. WGS was conducted at an average 25x coverage in a 2*150bp paired-end mode using the Illumina HiSeq4000. WGS-based SNP calling As a “gold-standard” for quality evaluation, WGS-based set of SNPs was used. A data quality check was conducted using FastQC [ 19 ]. Trimming and filtering were performed with the bbduk tool [ 13 ], using the following parameters: "k = 31 ref = artifacts,phix,adapters ordered cardinality qtrim = r trimq = 20 maq = 25 minlen = 50". The read mapping was performed using BWA-MEM [ 20 ] with default settings. The most recent Williams 82 assembly was used as the reference (Wm82.a4, https://www.soybase.org ). Conversion between SAM and BAM formats was performed using samtools 1.14 [ 14 ]. SNP calling of WGS set was conducted using GATK 4.2.6.1 HaplotypeCaller with the default settings [ 15 ]. This caller was used as it is reported to be the optimal choice for SNP calling using WGS data in plant species [ 21 ]. The output of this tool in gVCF format was then utilized as the input for joint genotyping performed by two additional tools: GATK GenotypeGVCFs and GLnexus [ 22 ]. In case of GenotypeGVCFs, all gVCF files were merged to one database per chromosome using GenomicsDBImport. GenotypeGVCFs was employed to obtain VCF files from each database. All VCF files were concatenated using bcftools 1.16. In the case of GLnexus, all gVCF files were merged in one step with the "--config gatk" parameter. The two resulting VCF files were compared and the intersection of SNPs was used in following analysis. The FORMAT and INFO fields were retained from the GATK VCF file. According to GATK guidelines we used following filters to remove low-quality SNPs: QD > 2, QUAL > 30, SOR < 3, FS 40. Then we applied additional filtering for depth (> 8 and < 80), MAF 0.05 and fraction of missing < 0.3. Pipeline for GBS data analysis A data quality check was performed using FastQC. In order to scale libraries by their read number, a random sampling of 18 million reads was performed using the reformat.sh tool from BBTools. Following this, reads were demultiplexed using the process_radtags tool from Stacks 2.66 [ 12 ]. Then, trimming and filtering were conducted using the bbduk tool with following parameters: " k = 31 ref = artifacts,phix,adapters,lambda,pjet,mtst,kapa ordered cardinality qtrim = rl trimq = 20 maq = 25 minlen = 50 tbo mink = 11 ktrim = r minlen = 50". The mapping was performed using BWA-MEM v.0.7.17-r1198-dirty [ 20 ], Bowtie2 v. 2.5.2 [ 23 ], BBMap v. 39.06 [ 24 ] or Strobealign v.0.13 [ 25 ] with default settings. The conversion between SAM and BAM formats was performed using the samtools 1.14. SNP calling was conducted using each one of seven callers: Bcftools 1.18 mpileup + call, Stacks 2.66 ref_map.pl, DeepVariant 1.6.0 [ 11 ] with glnexus 1.4.3, freebayes 1.3.6 [ 10 ], VarScan 2.4.6 [ 9 ] with samtools 1.19 mpileup, BBcallvariants 39.01, GATK 4.5.0.0 HaplotypeCaller. The default settings for each caller were employed, with the exception of BBcallvariants (border = 0, ploidy = 2), DeepVariant (--model_type = WGS) and VarScan (--min-coverage 2). The resulting vcf files were subjected to filtering using Bcftools 1.18, with genotypes with depth < 3 set to missing. Subsequently, biallelic single-nucleotide polymorphisms (SNPs) with an average depth per genotype of greater than five, a minor allele frequency (MAF) of 0.05, and a fraction of missing data of less than 0.4 were used for further analysis. Evaluation of SNP numbers using different-sized libraries In order to call SNP sets with different numbers of reads per library (ranging from 3 to 18 million paired reads), we performed random sampling of reads in triplicate using reformat.sh from BBTools. The obtained reads were then analyzed using BWA-MEM aligner and VarScan caller (--min-coverage 2) with SNP filtering described before. Metrics calculation The SNP chromosome and positions detected using GBS methods and the various combinations of aligners and callers were used for comparison with the reference SNP set computed using the WGS data. A pairwise comparison of the SNP sets was performed by Bcftools stats and Pearson's r^2 on dosage (number of non-ref alleles in genotype), non-reference discordance (NRD) and SNP number were used. $$\:NRD\:=\:100\:*\frac{\left(xRR\:+\:xRA\:+\:xAA\right)}{xRR\:+\:xRA\:+\:xAA\:+\:mRA\:+\:mAA}$$ The intersection with gene annotation was performed by Bedtools 2.30.0. [ 26 ]. Intersection of SNP sets was performed using the supervenn tool v 0.5.0 ( https://github.com/gecko984/supervenn ). Coverage estimation Fraction of covered regions was calculated using Bedtools 2.30.0 Data analysis All analyses were performed using R Statistical Software 4.3.1 [ 27 ]. dplyr, tidyr and ggplot2 were used from tidyverse [ 28 ]. Statistical computations were performed using packages car [ 29 ] and rstatix [ 30 ]. Results Enzyme set comparison With an equal number of reads per library (18 million paired-end), we ran full pipelines and compared libraries between each other. Firstly, we compared libraries by the total number of obtained SNPs, their presence in genes and the fraction of the SNP set identified using WGS data. We observed that libraries vary in number of genes covered (Fig. 1 A), consistent with proportion of covered regions in genome (Additional file 2). A comparison of the SNP sets with the SNPs identified using WGS data reveals a significant degree of variation within each of the evaluated enzyme sets (Fig. 1 B). This result is influenced by the accuracy of the SNP caller and will be discussed below. With regard to the location of the SNPs, the number of SNPs in genes correlates with the total number of SNPs for ApeKI and PstI - MspI libraries, but not for HindIII - NlaIII library (Fig. 1 C). Figure 1 D illustrates that observed case with the HindIII - NlaIII library. Despite this library is exhibiting the highest number of SNPs, the majority of them are located outside annotated genes. We also evaluated the same metrics using the full amount of reads, i.e. not scaled to 18 million per library, and found that the total number of SNPs and number of SNPs in genes are sensitive to genome coverage (Additional file 3), while both fraction metrics remained similar. In order to evaluate the dependency, we called SNPs using the BWA-VarScan pipeline with different read coverage per library, generated via sampling reads from the full set of reads (Fig. 2 ). The PstI - MspI library reaches the plateau in the total number of SNPs quickly than other libraries and further increase in read number does not affect the other metrics greatly. The ApeKI library is close to the plateau and the HindIII - NlaIII library exhibits an exponential growth. It is interesting to note that the fraction of SNPs in genes remains constant in different read numbers. Aligners performance A comparative evaluation was conducted of three popular aligner tools (BBMap, Bowtie2, BWA-MEM) and recent Strobealign across all pipelines, with the different enzyme sets and callers. The amount of called SNPs after filtering and the ratio of SNPs presented in WGS were calculated (Fig. 3 ). The results demonstrated that difference in aligner performances is not significant according to ANOVA analysis. Callers performance We compared performance of seven popular SNP callers in terms of the total number of SNPs and the presence of the corresponding called SNPs in WGS data (Fig. 4 ). The results demonstrated a significant disparity in the performance of the callers, with the outcomes being largely consistent across all three enzyme combinations. DeepVariant consistently produces a low number of SNPs (28,8% of maximum number of SNPs on average), that have a high presence in WGS-based SNPs (0.76 on average). In contrast, freebayes outputs a considerable number of unique SNPs that are less frequently presented in WGS. Comparison of SNP sets obtained using different enzymatic combinations and different callers Nevertheless, a direct comparison of GBS data with WGS is an open question, as the quality of the latter also depends on the caller used. Consequently, we have calculated recall, precision and false discovery rate (FDR) with respect to GBS data (Additional file 4). Intersections have been made between all datasets in each Enzyme set–Aligner combination to calculate the quantities of SNPs. For illustrative purposes, we present such intersection in Fig. 5 . For all pipelines, we defined false positive (FP) as SNP positions that were identified by one or two callers out of seven, and true positive (TP) as SNP positions were called by six or seven callers. SNPs detected using WGS data are also influenced by the 'caller effect,' so we deliberately avoid considering the WGS-based SNP set as the 'gold standard' reference. Instead, we rely exclusively on GBS data for editing results. We then compared the obtained data with respect to the enzyme set used and aligners (Fig. 6 ). In terms of accuracy DeepVariant demonstrated the best performance in Precision and FDR (on average 0.957 and 0.0095, respectively). Conversely, freebayes exhibited the lowest (0.365 precision and 0.6321 FDR). However, the recall of DeepVariant was the lowest achieved, with an average of 0.623. In order to estimate both the correct position identification and the genotype called for each sample, we performed pairwise estimation of R^2 dosage and NRD for each obtained SNP set. The results are presented in Additional file 5. We focused on comparisons within enzyme sets and calculated summary statistics of genotyping accuracy between with respect to enzyme set, aligner and caller used (Fig. 7 ). Genotyping accuracy with respect to enzyme set (Fig. 7 A) is correlated with the total number of SNPs (Fig. 1 A). The HindIII – NlaIII set demonstrates the lowest genotyping accuracy (average R^2 dosage 0.956). The highest accuracy was demonstrated by DeepVariant and Stacks with average R^2 dosage 0.988 and 0.982, respectively. Additionally, the range of R^2 dosage values was the narrowest among the seven callers. Discussion The choice of enzyme set and the design of SNP calling pipelines represent both crucial aspects of GBS genotyping. However, they are poorly discussed especially in addressing possible optimum of combining both aspects. In this study, we sought to address this gap in literature by genotyping several libraries prepared with different enzyme sets. Our findings demonstrate that the choice of enzyme and caller can significantly influence the final number of SNPs and their quality and we may provide some recommendations for soybean genotyping. Choice of enzyme set The selection of an appropriate enzyme set is the most important step as it cannot be modified subsequently by switching between tools and re-analyses. We conducted an evaluation of three popular enzyme sets, comparing their respective metrics. It is clear that the total number of SNPs is correlated with genome coverage. The greater the number of covered regions, the greater the number of SNPs that are identified. As illustrated in Additional file 2, there is a considerable disparity in the fraction of the genome covered by three libraries. A similar variation is observed in Fig. 1 A. From this perspective, it is therefore preferable to achieve wider coverage, which can be obtained through the use of enzyme sets like HindIII - NlaIII . However, the GBS method was developed with the objective of reducing genotyping costs, which means that the use of narrow coverage enzyme set may be more appropriate in cases where limited resources are available. As expected, the number of SNPs and genome coverage are dependent on the read number, but only to a certain extent. At some point, it is no longer possible to cover any additional regions and only the quality of the genotypes is improving due to increasing depth of sequencing. In our experimental data, we clearly observed this phenomenon in the PstI - MspI enzyme set, which can be considered an additional point for cost optimization. Upon further investigation of SNP metrics, it becomes evident that they also differ in other parameters, such as gene localization and the occurrence of SNPs in WGS (Fig. 1 ). On the one hand, ApeKI and PstI - MspI libraries demonstrate effective enrichment of gene sequences (Fig. 1 D). On the other hand, ApeKI library shows a smaller intersection with WGS than other enzyme sets, which may decrease our confidence in identified SNPs. If our main interest relies in LD between SNPs markers and QTLs, as it is the key for GWAS and development of molecular markers, enrichment of SNPs in gene sequences is not of upmost importance. An additional point to consider is the fact that three evaluated enzyme sets result in often non-overlapping sets of SNPs (Additional file 5, Fig. 5 ). This means that they cover different regions of the genome and data from distinct accessions obtained using different enzyme sets cannot be concatenated. Choice of aligner Alignment is the intermediate but important step in the SNP calling process. Accurate read positioning is critical for high-quality SNP calling. It appears that the role of the aligner is frequently underestimated. Bowtie2 or BWA are frequently employed in a variety of pipelines without undergoing a comprehensive evaluation. Our results suggest that the choice of aligner has a limited effect on SNP quantity, and their quality and performance is similar. Nevertheless, we notice that the Strobealign aligner, while maintaining similar performance, is several times faster than other aligners tested. This leads to a better use of computational resources and thus to a cost reduction. Choice of caller The role of caller is probably the most crucial in the pipeline, as it provides the final output for many other analyses. Several callers that have been created recently or updated regularly, which have been used in human genetics. However, they have not been widely used in plant genetics, particularly in the context of GBS. In this study, we compared six popular callers and the Stacks pipeline which has been designed for GBS. We also used both classical callers and deep learning based caller. To perform analyses such as GWAS, it is necessary to deal with a reliable and consistent set of SNPs. In order to assess this, we called SNPs from the same read sets using different callers. The results demonstrated that the sets of SNPs obtained with using the different callers varied in total number of SNPs as well as their occurrence in the WGS-based SNP set. These results demonstrate that callers may discover SNPs from the same data in a different way, which affects the final result. This can be explained by the fundamentals of the SNP calling process: the caller should efficiently process the aligned reads, which may contain errors, and should be able to take into account genomes with different characteristics. To take into account sources of errors and genome parameters, the callers apply different comprehensive statistical models. However, it is not yet perfect and we again illustrate here that callers may produce SNPs that are unique to a particular caller. The ability of deep learning to capture underlying unnoticed dependencies is widely used in data science and has been shown to be applicable to variant calling. Despite the low number of SNPs called by DeepVariant, they can be found in WGS with a high probability. Because identification of SNPs from WGS data can be also affected by the caller effects, we calculated precision, recall and FDR based only on comparing GBS data. In this instance, we again observe that DeepVariant outperforms other callers in terms of precision and FDR, but recall is not so good. With regard to genotyping accuracy, DeepVariant also exhibits the highest R^2 dosage and lower variation in pairwise comparison. Our analyses raise the question of which caller is preferable. The answer of course depends on the purpose, species and resources available. In some cases, it may be reasonable to obtain the most reliable set in order to be confident in the SNP set and not to work with FP SNPs that do not exist in reality. However, in this way we may miss some valuable SNPs due to the stricter caller procedure. We postulate that the described workflow can be easily adapted to other plant species, with the aim of identifying an optimal combination of molecular experiments and bioinformatics pipelines for a modest cost, as a preliminary step before large-scale analyses. Conclusions The choice of enzyme set and the design of SNP calling pipelines represent are both crucial aspects of GBS genotyping. However, they are poorly discussed especially in addressing possible optimum of combining both aspects. Our results show that the choice of enzymes significantly influences the number of SNPs identified, gene localisation preferences and accuracy. Furthermore, the performance of SNP callers varied significantly in terms of SNP number, precision, recall and false discovery rate (FDR). These results highlight the importance of optimising both enzyme selection for sequencing libraries and data analysis pipelines to ensure robust and reproducible results. Declarations Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and materials The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. Competing interests The authors declare that they have no competing interests Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript. Authors' contributions AZ performed the analysis and drafted the manuscript. SB prepared sequencing libraries. CB guided the research. LG guided the research and revised the manuscript. All authors read and approved the final manuscript. Acknowledgements Not applicable References Thomson MJ (2014) High-Throughput SNP Genotyping to Accelerate Crop Improvement. Plant Breed Biotechnol 2:195–212 Moragues M, Comadran J, Waugh R, Milne I, Flavell AJ, Russell JR (2010) Effects of ascertainment bias and marker number on estimations of barley diversity from high-throughput SNP genotype data. Theoretical and Applied Genetics 120:1525–1534 Tranchant‐Dubreuil C, Rouard M, Sabot F (2019) Plant Pangenome: Impacts on Phenotypes and Evolution. In: Annual Plant Reviews online. Wiley, pp 453–478 Warr A, Robert C, Hume D, Archibald A, Deeb N, Watson M (2015) Exome sequencing: Current and future perspectives. G3: Genes, Genomes, Genetics 5:1543–1550 Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA (2008) Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One. https://doi.org/10.1371/journal.pone.0003376 Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE (2011) A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. https://doi.org/10.1371/journal.pone.0019379 Heffelfinger C, Fragoso CA, Moreno MA, Overton JD, Mottinger JP, Zhao H, Tohme J, Dellaporta SL (2014) Flexible and scalable genotyping-by-sequencing strategies for population studies. BMC Genomics. https://doi.org/10.1186/1471-2164-15-979 Shin MG, Ithnin M, Vu WT, Kamaruddin K, Chin TN, Yaakub Z, Chang PL, Sritharan K, Nuzhdin S, Singh R (2021) Association mapping analysis of oil palm interspecific hybrid populations and predicting phenotypic values via machine learning algorithms. Plant Breeding 140:1150–1165 Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK (2012) VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res 22:568–576 Garrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. Poplin R, Chang P-C, Alexander D, et al (2018) A universal SNP and small-indel variant caller using deep neural networks. Nat Biotechnol 36:983–987 Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA (2013) Stacks: an analysis tool set for population genomics. Mol Ecol 22:3124–3140 Bushnell B, Rood J, Singer E (2017) BBMerge – Accurate paired shotgun read merging via overlap. PLoS One 12:e0185056 Danecek P, Bonfield JK, Liddle J, et al (2021) Twelve years of SAMtools and BCFtools. Gigascience. https://doi.org/10.1093/gigascience/giab008 Poplin R, Ruano-Rubio V, DePristo MA, et al (2018) Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv 201178 Taniguti CH, Taniguti LM, Amadeu RR, et al (2023) Developing best practices for genotyping-by-sequencing analysis in the construction of linkage maps. Gigascience. https://doi.org/10.1093/gigascience/giad092 Gupta SK, Baek J, Carrasquilla-Garcia N, Penmetsa RV (2015) Genome-wide polymorphism detection in peanut using next-generation restriction-site-associated DNA (RAD) sequencing. Molecular Breeding. https://doi.org/10.1007/s11032-015-0343-0 Poland J (2011) PstI-­MspI GBS Genotyping-­by-­sequencing Protocol PstI-­MspI. Andrews S. (2010) FastQC: a quality control tool for high throughput sequence data. Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Wu X, Heffelfinger C, Zhao H, Dellaporta SL (2019) Benchmarking variant identification tools for plant diversity discovery. BMC Genomics. https://doi.org/10.1186/s12864-019-6057-7 Yun T, Li H, Chang PC, Lin MF, Carroll A, McLean CY (2020) Accurate, scalable cohort variant calls using DeepVariant and GLnexus. Bioinformatics 36:5582–5589 Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9:357–359 Bushnell B (2014) BBMap: a fast, accurate, splice-aware aligner. Sahlin K (2022) Strobealign: flexible seed size enables ultra-fast and accurate read alignment. Genome Biol 23:260 Quinlan AR, Hall IM (2010) BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26:841–842 R Core Team (2021) R: A Language and Environment for Statistical Computing. Wickham H, Averick M, Bryan J, et al (2019) Welcome to the Tidyverse. J Open Source Softw 4:1686 John Fox and Sanford Weisberg (2019) An R Companion to Applied Regression, third. Sage Alboukadel Kassambara (2023) rstatix: Pipe-Friendly Framework for Basic Statistical Tests. Additional Declarations No competing interests reported. Supplementary Files Additionalfile1.csv Additional file 1.csv List of soybean varieties used in the study Additionalfile2.png Additional file 2.png Cumulative genome coverage by different enzyme sets. X-axis represents fraction of the genome covered by particular number of reads and Y-axis represents how many such regions in genome plus previous fractions (cumulative). 1.2 mln of reads were used for coverage estimation and this number of reads corresponds to average value per accession for a library of 15 pooled accessions, and containing 18 mln reads in total. Additionalfile3.png Additional file 3.png Comparison of different enzyme combinations in GBS using unscaled number of reads. Each pipeline combination was tested once and plot represent data of these tests. A – Total number of SNPs after filtering depending on used enzyme set. B – fraction of SNPs that were observed also in WGS data for the same varieties. C – Number of SNPs located in genes according to genome annotation available. D – fraction of SNPs located in genes. Additionalfile4.csv Additional file 4.csv Confusion matrix for all tested pipelines. Recall, precision and false discovery rate (FDR) were calculated with respect to GBS data. Intersections have been made between all datasets in each Enzyme set–Aligner combination Additionalfile5.png Additional file 5.png R^2, NRD and the number of evaluated SNPs in each pair. Upper left part represents R^2 and the number of evaluated SNPs in each pair. Blue gradient represents R^2 where the lightest blue corresponds to maximum value. Lower right part represents NRD and the lightest grey corresponds to minimum NRD. Cite Share Download PDF Status: Published Journal Publication published 06 Aug, 2025 Read the published version in Plant Methods → Version 1 posted Editorial decision: Revision requested 29 Apr, 2025 Reviews received at journal 27 Apr, 2025 Reviewers agreed at journal 15 Apr, 2025 Reviews received at journal 17 Feb, 2025 Reviewers agreed at journal 10 Feb, 2025 Reviewers invited by journal 05 Feb, 2025 Editor assigned by journal 17 Jan, 2025 Submission checks completed at journal 17 Jan, 2025 First submitted to journal 13 Jan, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-5821852","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":403723175,"identity":"4a0bb5c7-c666-4696-a233-5efd3e027f79","order_by":0,"name":"Aleksei Zamalutdinov","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABLUlEQVRIie2QMWuDQBTHnwi6vOB6pSH5CicFQ1HwqxwEziVDuxQ6hSDYxblQ6Jfo3sFyYBbpLNwiBJwymM2hhZ6GdqjaucP94HgPHr/7vzsAjeYfYtnnulXHyABM7JoKMph30zHFMc+1LxmwXjGpahEnlIv4l9JHk78UKsycnBKwnEdRiVnrz1d26t23RYChHdeV8SqGirV2nxJAIjkVyCK8TgtPYskRMV9RqEcUdA+zBAhI1ikCacm5hEYgkk234UAJhdOIzwToUkbNj3LbdsryOKqom42DkQCjcvOdss4Byy4FJxTryk3fCXPl5ubtmUdICyEusVBvKbhHWBYNlH1ck/Yu2C5k9FIdAz+k+93u1ObBwn4QNWkyf+yjwVALTMCmBh9TA41Go9EovgA+OW/eIljl8gAAAABJRU5ErkJggg==","orcid":"","institution":"Skolkovo Institute of Science and Technology","correspondingAuthor":true,"prefix":"","firstName":"Aleksei","middleName":"","lastName":"Zamalutdinov","suffix":""},{"id":403723176,"identity":"2bbbb985-eb34-447c-b86f-4e0ff31868d3","order_by":1,"name":"Stepan Boldyrev","email":"","orcid":"","institution":"Skolkovo Institute of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Stepan","middleName":"","lastName":"Boldyrev","suffix":""},{"id":403723177,"identity":"e422e221-fcce-4107-bb70-0abb8096c5b3","order_by":2,"name":"Cécile Ben","email":"","orcid":"","institution":"Skolkovo Institute of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Cécile","middleName":"","lastName":"Ben","suffix":""},{"id":403723178,"identity":"e4cdd071-a311-4064-ab4e-5fc5e17346f7","order_by":3,"name":"Laurent Gentzbittel","email":"","orcid":"","institution":"Skolkovo Institute of Science and Technology","correspondingAuthor":false,"prefix":"","firstName":"Laurent","middleName":"","lastName":"Gentzbittel","suffix":""}],"badges":[],"createdAt":"2025-01-13 17:23:21","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-5821852/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-5821852/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1186/s13007-025-01410-8","type":"published","date":"2025-08-06T15:57:18+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":74238714,"identity":"d6e2c8d0-4ced-4131-b6e7-f231f6995b6c","added_by":"auto","created_at":"2025-01-20 09:15:21","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":150680,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of different enzyme combinations in GBS.\u003c/strong\u003e Each pipeline combination was tested once and plot represent data of these tests. A – Total number of SNPs after filtering depending on used enzyme set. B – fraction of SNPs that were observed also in WGS data for the same varieties. C – Number of SNPs located in genes according to genome annotation available. D – fraction of SNPs located in genes. \u003cbr\u003e\nFor B and D, ANOVA was performed using arcsin square root transformed data.\u003c/p\u003e","description":"","filename":"Figure1.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/0d4c7de1e640a977552d16db.jpg"},{"id":74239372,"identity":"1eb9618c-1d2b-4c3e-9384-31409e023d99","added_by":"auto","created_at":"2025-01-20 09:23:21","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":232031,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDependency of SNP metrics using different amounts of read per library.\u003c/strong\u003e A – Total number of SNPs after filtering. B – fraction of SNPs that were observed also in WGS data for the same varieties. C – Number of SNPs located in genes according to genome annotation available. D – fraction of SNPs located in genes.For this analysis, only the BWA-VarScan pipeline was used.\u003c/p\u003e","description":"","filename":"Figure2.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/614d0f085d43069932fa8656.jpg"},{"id":74239416,"identity":"f12594a5-5674-493c-8cb7-450992a7bff0","added_by":"auto","created_at":"2025-01-20 09:23:22","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":196122,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of different aligners in all enzyme sets.\u003c/strong\u003e Each pipeline combination was tested once and plot represent data of these tests. Differences between aligner performances are not significant. ANOVA was performed on arcsin square root transformed data. A – Amount of SNPs after filtering. B – fraction of SNPs that were observed also in WGS data for the same varieties.\u003c/p\u003e","description":"","filename":"Figure3.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/cb638f077a7e4956f0930646.jpg"},{"id":74238730,"identity":"41589f4b-eb9a-4181-a12d-354233e6cb57","added_by":"auto","created_at":"2025-01-20 09:15:22","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":170353,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eComparison of callers’ performance.\u003c/strong\u003e Each pipeline combination was tested once and plot represent data of these tests. Tukey post hoc test was performed to distinguish treatment groups. ANOVA was performed on arcsin square root transformed data. A – Amount of SNPs after filtering. Values represent number of SNPs divided by maximum number of obtained SNPs per enzyme set. ANOVA(6,77)=30.4, p\u0026lt;0.0001, B – fraction of SNPs that were observed also in WGS data for the same varieties ANOVA(6,77)=15.76, p\u0026lt;0.0001.\u003c/p\u003e","description":"","filename":"Figure4.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/a2c5ea01e9732b45949eee18.jpg"},{"id":74238717,"identity":"0373fbe7-1eb5-4e46-8e2d-c2091c01d485","added_by":"auto","created_at":"2025-01-20 09:15:21","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":407035,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eVisualisation of intersections in SNP sets between 7 callers. \u003c/strong\u003eBWA was used as the common aligner for all callers. It is venn-like diagram, where first row shows how many callers called the same SNPs, seven next rows show which caller (identified by different colours) participate in this intersection, the last row contains the number of SNPs in intersection. Last column contains total number of SNPs identified by particular caller. A – ApeKI enzyme set. B – HindIII–NlaIII enzyme set. C – PstI–MspI enzyme set.\u003c/p\u003e","description":"","filename":"Figure5.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/ddaf9e675258974b9d101464.jpg"},{"id":74238741,"identity":"81c7bcf5-d9b1-4312-bcdd-e082f8891ab0","added_by":"auto","created_at":"2025-01-20 09:15:22","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":144202,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eEvaluation of precision, recall and FDR of different callers.\u003c/strong\u003e For each combination of aligner and enzyme set TP and FP were defined independently. Each pipeline combination was tested once and plot represent data of these tests. Tukey post hoc test was performed to distinguish treatment groups. ANOVA was performed on arcsin square root transformed data. A – precision Anova(6,77)=26.75,p\u0026lt;0.0001, B – recall Anova(6,77)=17.83,p\u0026lt;0.0001, C – FPR Anova(6,77)=101.7,p\u0026lt;0.0001.\u003c/p\u003e","description":"","filename":"Figure6.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/42a095df0daac49b2536adcb.jpg"},{"id":74238745,"identity":"6bffccc4-5458-4673-a5a6-64c7737b3d26","added_by":"auto","created_at":"2025-01-20 09:15:25","extension":"jpg","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":205862,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDependency of R^2 dosage on enzyme set, aligner and caller used.\u003c/strong\u003e R^2 dosage was calculated by comparison of genotypes called by different pipelines. We compared enzyme sets which were obtained in one enzyme set, not between different enzyme sets. A – R^2 dosage depending on enzyme set used. B – R^2 dosage depending on aligner used. C – R^2 dosage depending on caller used. 6 most significant pairwise differences are shown.\u003c/p\u003e","description":"","filename":"Figure7.tiff.jpg","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/9af633dd69a5d8e1b238ce52.jpg"},{"id":88814110,"identity":"f14890f4-8dd2-4e74-92df-6300cd0d6d1a","added_by":"auto","created_at":"2025-08-11 16:06:59","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2092083,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/cb00ff00-ca9e-4ac5-8071-6fc160dd4359.pdf"},{"id":74239371,"identity":"8e444890-2e5f-464c-bc83-1569d7985b6b","added_by":"auto","created_at":"2025-01-20 09:23:21","extension":"csv","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":217,"visible":true,"origin":"","legend":"\u003cp\u003eAdditional file 1.csv \u003cstrong\u003eList of soybean varieties used in the study\u003c/strong\u003e\u003c/p\u003e","description":"","filename":"Additionalfile1.csv","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/f771c7b70e5b91d48e799102.csv"},{"id":74238705,"identity":"040e1427-dd4b-44f7-b8b9-e3b46d9e08bf","added_by":"auto","created_at":"2025-01-20 09:15:20","extension":"png","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":136216,"visible":true,"origin":"","legend":"\u003cp\u003eAdditional file 2.png \u003cstrong\u003eCumulative genome coverage by different enzyme sets\u003c/strong\u003e. X-axis represents fraction of the genome covered by particular number of reads and Y-axis represents how many such regions in genome plus previous fractions (cumulative). 1.2 mln of reads were used for coverage estimation and this number of reads corresponds to average value per accession for a library of 15 pooled accessions, and containing 18 mln reads in total.\u003c/p\u003e","description":"","filename":"Additionalfile2.png","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/35205075e783e208c365e3e1.png"},{"id":74238704,"identity":"26ec3767-f483-4c65-9862-484571da0c14","added_by":"auto","created_at":"2025-01-20 09:15:20","extension":"png","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":93639,"visible":true,"origin":"","legend":"\u003cp\u003eAdditional file 3.png \u003cstrong\u003eComparison of different enzyme combinations in GBS using unscaled number of reads\u003c/strong\u003e. Each pipeline combination was tested once and plot represent data of these tests. A – Total number of SNPs after filtering depending on used enzyme set. B – fraction of SNPs that were observed also in WGS data for the same varieties. C – Number of SNPs located in genes according to genome annotation available. D – fraction of SNPs located in genes.\u003c/p\u003e","description":"","filename":"Additionalfile3.png","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/3be74dff98bc74bdd97cca99.png"},{"id":74238706,"identity":"d2192407-c396-49e1-b68f-fcbd089153b9","added_by":"auto","created_at":"2025-01-20 09:15:20","extension":"csv","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":9313,"visible":true,"origin":"","legend":"\u003cp\u003eAdditional file 4.csv \u003cstrong\u003eConfusion matrix for all tested pipelines\u003c/strong\u003e. Recall, precision and false discovery rate (FDR) were calculated with respect to GBS data. Intersections have been made between all datasets in each Enzyme set–Aligner combination\u003c/p\u003e","description":"","filename":"Additionalfile4.csv","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/359baaea3fd4eea24d44bd90.csv"},{"id":74238735,"identity":"98ee0c69-c4d9-4ee4-a533-7ce8d771d3de","added_by":"auto","created_at":"2025-01-20 09:15:22","extension":"png","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":17992677,"visible":true,"origin":"","legend":"\u003cp\u003eAdditional file 5.png \u003cstrong\u003eR^2, NRD and the number of evaluated SNPs in each pair.\u003c/strong\u003e Upper left part represents R^2 and the number of evaluated SNPs in each pair. Blue gradient represents R^2 where the lightest blue corresponds to maximum value. Lower right part represents NRD and the lightest grey corresponds to minimum NRD.\u003c/p\u003e","description":"","filename":"Additionalfile5.png","url":"https://assets-eu.researchsquare.com/files/rs-5821852/v1/b1c7cbf85239e421c2ad8a49.png"}],"financialInterests":"No competing interests reported.","formattedTitle":"The evaluation of different combinations of enzyme set, aligner and caller in GBS sequencing of soybean","fulltext":[{"header":"Background","content":"\u003cp\u003eGenotyping is a crucial step in plant and animal breeding, as well as ecological population studies. It can be carried out in a number of ways, including microarray chips, WGS and reduced representation sequencing. Microarrays represent a large-scale and cost-effective approach that is widely used in human genetics and animal breeding. However, array design is a costly process [\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e], it works with a fixed set of polymorphisms, cannot identify new genetic variation in diverse populations and is susceptible to structural variation [\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e]. Furthermore, microarray chips are proprietary tools, which places researchers in a position of dependency with regard to the manufacturer and supplier, as well as their sales policies. WGS is the most informative approach in sequencing as it provides the genome-wide data for organism. WGS allows not only to retrieve information about polymorphisms in the entire genome but also to benefit from available pangenome assemblies and to efficiently study dispensable fractions of the genome [\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e]. However, it requires high coverages (up to 30x) to obtain accurate results, which can be prohibitively costly for large eukaryotic genomes or in case of genotyping large collections of accessions. As an alternative, several reduced representation sequencing (RRS) approaches have been developed to reduce the amount of required reads by enriching, separating, or eliminating a portion of the genome prior to sequencing. The majority of these approaches are costly (such as exome sequencing) and require the use of specialized commercial kits and high-quality annotation of the reference genome [\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eAn alternative approach is to reduce genome complexity through the use of restriction enzymes. Restriction-site associated (RAD) [\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e] or Genotype-by-sequencing (GBS) [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e] families of methods employ the use of one or more restriction enzymes for the digestion of genome DNA, followed by the adapter ligation. This approach is easy to implement due to the availability of published protocols and restriction enzymes at a low cost.\u003c/p\u003e \u003cp\u003eA key question for the GBS method is the selection of an appropriate enzyme set. This is a complex problem with different parameters to consider such as differences in restriction site frequencies in different genomes, the rate of fragment occurrence in repetitive regions and restriction site methylation. It has been demonstrated that the use of methylation-sensitive enzymes allows for a better reduction of genome complexity and better coverage of gene-containing regions [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e, \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e]. Nevertheless, it is worth noting that methylation-insensitive enzyme pairs are also used in GBS studies such as \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e [\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e]. An important area for potential improvement is the bioinformatics pipeline for data analysis. Variations in read quality and genome characteristics can influence the approaches employed for analysis. There are numerous tools available for read filtering, mapping (Bushnell, 2014; Langmead \u0026amp; Salzberg, 2012; Li, 2013) and polymorphism calling [\u003cspan additionalcitationids=\"CR10 CR11 CR12 CR13 CR14\" citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e–\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. Noteworthy, it has been demonstrated that the set of obtained SNPs is highly dependent on the caller that is employed [\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e].\u003c/p\u003e \u003cp\u003eIn this study, three libraries of 15 \u003cem\u003eGlycine max\u003c/em\u003e varieties were prepared using three different popular enzyme combinations: \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e, \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e and \u003cem\u003eApeKI\u003c/em\u003e, and also subjected to 25x WGS. The libraries were tested with multiple pipelines using three mapping tools (BWA, Bowtie2 and BBMap) and seven SNP callers. The quality of the obtained SNPs was evaluated in different tool combinations by comparison with each other and with respect to WGS-based SNPs dataset identified within the 15 varieties, and using the most recent genome annotation. We show that the choice of enzyme combination significantly affects the number and localization of the SNPs obtained. The influence of aligner is not significant in case of total SNPs numbers and intersection SNPs set detected using WGS data. It is also of minor significance in genotyping accuracy. Callers differ in all evaluated criteria and DeepVariant showed the highest genotyping accuracy and lowest false positive rate.\u003c/p\u003e \u003cp\u003e\u003c/p\u003e "},{"header":"Methods","content":"\u003cp\u003eGenotyping\u003c/p\u003e\u003cp\u003eWe selected a total of 15 soybean varieties of different world regions and with realized economic potential. All accessions were genotyped using GBS and WGS approaches (Additional file 1). GBS libraries were prepared for three enzyme combination of restriction enzymes according to established protocols (\u003cem\u003eApeKI\u003c/em\u003e [\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e], \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e [\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e], \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e [\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e]). The resulting libraries were then sequenced using the Illumina HiSeq4000 in a 2*150bp paired-end mode with 18.2, 31.5 and 19.5\u0026nbsp;million of reads per library, respectively. WGS was conducted at an average 25x coverage in a 2*150bp paired-end mode using the Illumina HiSeq4000.\u003c/p\u003e\u003cp\u003eWGS-based SNP calling\u003c/p\u003e\u003cp\u003eAs a “gold-standard” for quality evaluation, WGS-based set of SNPs was used. A data quality check was conducted using FastQC [\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e]. Trimming and filtering were performed with the bbduk tool [\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e], using the following parameters: \"k = 31 ref = artifacts,phix,adapters ordered cardinality qtrim = r trimq = 20 maq = 25 minlen = 50\". The read mapping was performed using BWA-MEM [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e] with default settings. The most recent Williams 82 assembly was used as the reference (Wm82.a4, \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://www.soybase.org\u003c/span\u003e\u003cspan address=\"https://www.soybase.org\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). Conversion between SAM and BAM formats was performed using samtools 1.14 [\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e].\u003c/p\u003e\u003cp\u003eSNP calling of WGS set was conducted using GATK 4.2.6.1 HaplotypeCaller with the default settings [\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e]. This caller was used as it is reported to be the optimal choice for SNP calling using WGS data in plant species [\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e]. The output of this tool in gVCF format was then utilized as the input for joint genotyping performed by two additional tools: GATK GenotypeGVCFs and GLnexus [\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e]. In case of GenotypeGVCFs, all gVCF files were merged to one database per chromosome using GenomicsDBImport. GenotypeGVCFs was employed to obtain VCF files from each database. All VCF files were concatenated using bcftools 1.16. In the case of GLnexus, all gVCF files were merged in one step with the \"--config gatk\" parameter. The two resulting VCF files were compared and the intersection of SNPs was used in following analysis. The FORMAT and INFO fields were retained from the GATK VCF file.\u003c/p\u003e\u003cp\u003eAccording to GATK guidelines we used following filters to remove low-quality SNPs: QD \u0026gt; 2, QUAL \u0026gt; 30, SOR \u0026lt; 3, FS \u0026lt; 60, MQ \u0026gt; 40. Then we applied additional filtering for depth (\u0026gt; 8 and \u0026lt; 80), MAF 0.05 and fraction of missing \u0026lt; 0.3.\u003c/p\u003e\u003cp\u003ePipeline for GBS data analysis\u003c/p\u003e\u003cp\u003eA data quality check was performed using FastQC. In order to scale libraries by their read number, a random sampling of 18\u0026nbsp;million reads was performed using the reformat.sh tool from BBTools. Following this, reads were demultiplexed using the process_radtags tool from Stacks 2.66 [\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e]. Then, trimming and filtering were conducted using the bbduk tool with following parameters: \" k = 31 ref = artifacts,phix,adapters,lambda,pjet,mtst,kapa ordered cardinality qtrim = rl trimq = 20 maq = 25 minlen = 50 tbo mink = 11 ktrim = r minlen = 50\".\u003c/p\u003e\u003cp\u003eThe mapping was performed using BWA-MEM v.0.7.17-r1198-dirty [\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e], Bowtie2 v. 2.5.2 [\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e], BBMap v. 39.06 [\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e] or Strobealign v.0.13 [\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e] with default settings. The conversion between SAM and BAM formats was performed using the samtools 1.14.\u003c/p\u003e\u003cp\u003eSNP calling was conducted using each one of seven callers: Bcftools 1.18 mpileup + call, Stacks 2.66 ref_map.pl, DeepVariant 1.6.0 [\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e] with glnexus 1.4.3, freebayes 1.3.6 [\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e], VarScan 2.4.6 [\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e] with samtools 1.19 mpileup, BBcallvariants 39.01, GATK 4.5.0.0 HaplotypeCaller. The default settings for each caller were employed, with the exception of BBcallvariants (border = 0, ploidy = 2), DeepVariant (--model_type = WGS) and VarScan (--min-coverage 2).\u003c/p\u003e\u003cp\u003eThe resulting vcf files were subjected to filtering using Bcftools 1.18, with genotypes with depth \u0026lt; 3 set to missing. Subsequently, biallelic single-nucleotide polymorphisms (SNPs) with an average depth per genotype of greater than five, a minor allele frequency (MAF) of 0.05, and a fraction of missing data of less than 0.4 were used for further analysis.\u003c/p\u003e\u003cp\u003eEvaluation of SNP numbers using different-sized libraries\u003c/p\u003e\u003cp\u003eIn order to call SNP sets with different numbers of reads per library (ranging from 3 to 18\u0026nbsp;million paired reads), we performed random sampling of reads in triplicate using reformat.sh from BBTools. The obtained reads were then analyzed using BWA-MEM aligner and VarScan caller (--min-coverage 2) with SNP filtering described before.\u003c/p\u003e\u003cp\u003eMetrics calculation\u003c/p\u003e\u003cp\u003eThe SNP chromosome and positions detected using GBS methods and the various combinations of aligners and callers were used for comparison with the reference SNP set computed using the WGS data. A pairwise comparison of the SNP sets was performed by Bcftools stats and Pearson's r^2 on dosage (number of non-ref alleles in genotype), non-reference discordance (NRD) and SNP number were used.\u003c/p\u003e\u003cdiv id=\"Equa\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equa\" name=\"EquationSource\"\u003e\n$$\\:NRD\\:=\\:100\\:*\\frac{\\left(xRR\\:+\\:xRA\\:+\\:xAA\\right)}{xRR\\:+\\:xRA\\:+\\:xAA\\:+\\:mRA\\:+\\:mAA}$$\u003c/div\u003e\u003c/div\u003e\u003cp\u003eThe intersection with gene annotation was performed by Bedtools 2.30.0. [\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e]. Intersection of SNP sets was performed using the supervenn tool v 0.5.0 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/gecko984/supervenn\u003c/span\u003e\u003cspan address=\"https://github.com/gecko984/supervenn\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eCoverage estimation\u003c/p\u003e\u003cp\u003eFraction of covered regions was calculated using Bedtools 2.30.0\u003c/p\u003e\u003ch2\u003eData analysis\u003c/h2\u003e\u003cp\u003eAll analyses were performed using R Statistical Software 4.3.1 [\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e]. dplyr, tidyr and ggplot2 were used from tidyverse [\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e]. Statistical computations were performed using packages car [\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e] and rstatix [\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e].\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003eEnzyme set comparison\u003c/p\u003e \u003cp\u003eWith an equal number of reads per library (18\u0026nbsp;million paired-end), we ran full pipelines and compared libraries between each other. Firstly, we compared libraries by the total number of obtained SNPs, their presence in genes and the fraction of the SNP set identified using WGS data. We observed that libraries vary in number of genes covered (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA), consistent with proportion of covered regions in genome (Additional file 2). A comparison of the SNP sets with the SNPs identified using WGS data reveals a significant degree of variation within each of the evaluated enzyme sets (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB). This result is influenced by the accuracy of the SNP caller and will be discussed below. With regard to the location of the SNPs, the number of SNPs in genes correlates with the total number of SNPs for \u003cem\u003eApeKI\u003c/em\u003e and \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e libraries, but not for \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e library (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eC). Figure\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eD illustrates that observed case with the \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e library. Despite this library is exhibiting the highest number of SNPs, the majority of them are located outside annotated genes. We also evaluated the same metrics using the full amount of reads, i.e. not scaled to 18\u0026nbsp;million per library, and found that the total number of SNPs and number of SNPs in genes are sensitive to genome coverage (Additional file 3), while both fraction metrics remained similar. In order to evaluate the dependency, we called SNPs using the BWA-VarScan pipeline with different read coverage per library, generated via sampling reads from the full set of reads (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). The \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e library reaches the plateau in the total number of SNPs quickly than other libraries and further increase in read number does not affect the other metrics greatly. The \u003cem\u003eApeKI\u003c/em\u003e library is close to the plateau and the \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e library exhibits an exponential growth. It is interesting to note that the fraction of SNPs in genes remains constant in different read numbers.\u003c/p\u003e \u003cp\u003eAligners performance\u003c/p\u003e \u003cp\u003eA comparative evaluation was conducted of three popular aligner tools (BBMap, Bowtie2, BWA-MEM) and recent Strobealign across all pipelines, with the different enzyme sets and callers. The amount of called SNPs after filtering and the ratio of SNPs presented in WGS were calculated (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e). The results demonstrated that difference in aligner performances is not significant according to ANOVA analysis.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eCallers performance\u003c/p\u003e \u003cp\u003eWe compared performance of seven popular SNP callers in terms of the total number of SNPs and the presence of the corresponding called SNPs in WGS data (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The results demonstrated a significant disparity in the performance of the callers, with the outcomes being largely consistent across all three enzyme combinations. DeepVariant consistently produces a low number of SNPs (28,8% of maximum number of SNPs on average), that have a high presence in WGS-based SNPs (0.76 on average). In contrast, freebayes outputs a considerable number of unique SNPs that are less frequently presented in WGS.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eComparison of SNP sets obtained using different enzymatic combinations and different callers\u003c/p\u003e \u003cp\u003eNevertheless, a direct comparison of GBS data with WGS is an open question, as the quality of the latter also depends on the caller used. Consequently, we have calculated recall, precision and false discovery rate (FDR) with respect to GBS data (Additional file 4). Intersections have been made between all datasets in each Enzyme set\u0026ndash;Aligner combination to calculate the quantities of SNPs. For illustrative purposes, we present such intersection in Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFor all pipelines, we defined false positive (FP) as SNP positions that were identified by one or two callers out of seven, and true positive (TP) as SNP positions were called by six or seven callers. SNPs detected using WGS data are also influenced by the 'caller effect,' so we deliberately avoid considering the WGS-based SNP set as the 'gold standard' reference. Instead, we rely exclusively on GBS data for editing results. We then compared the obtained data with respect to the enzyme set used and aligners (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e). In terms of accuracy DeepVariant demonstrated the best performance in Precision and FDR (on average 0.957 and 0.0095, respectively). Conversely, freebayes exhibited the lowest (0.365 precision and 0.6321 FDR). However, the recall of DeepVariant was the lowest achieved, with an average of 0.623.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn order to estimate both the correct position identification and the genotype called for each sample, we performed pairwise estimation of R^2 dosage and NRD for each obtained SNP set. The results are presented in Additional file 5. We focused on comparisons within enzyme sets and calculated summary statistics of genotyping accuracy between with respect to enzyme set, aligner and caller used (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e). Genotyping accuracy with respect to enzyme set (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003eA) is correlated with the total number of SNPs (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA). The \u003cem\u003eHindIII\u003c/em\u003e\u0026ndash;\u003cem\u003eNlaIII\u003c/em\u003e set demonstrates the lowest genotyping accuracy (average R^2 dosage 0.956). The highest accuracy was demonstrated by DeepVariant and Stacks with average R^2 dosage 0.988 and 0.982, respectively. Additionally, the range of R^2 dosage values was the narrowest among the seven callers.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eThe choice of enzyme set and the design of SNP calling pipelines represent both crucial aspects of GBS genotyping. However, they are poorly discussed especially in addressing possible optimum of combining both aspects. In this study, we sought to address this gap in literature by genotyping several libraries prepared with different enzyme sets. Our findings demonstrate that the choice of enzyme and caller can significantly influence the final number of SNPs and their quality and we may provide some recommendations for soybean genotyping.\u003c/p\u003e \u003cp\u003eChoice of enzyme set\u003c/p\u003e \u003cp\u003eThe selection of an appropriate enzyme set is the most important step as it cannot be modified subsequently by switching between tools and re-analyses. We conducted an evaluation of three popular enzyme sets, comparing their respective metrics. It is clear that the total number of SNPs is correlated with genome coverage. The greater the number of covered regions, the greater the number of SNPs that are identified. As illustrated in Additional file 2, there is a considerable disparity in the fraction of the genome covered by three libraries. A similar variation is observed in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA. From this perspective, it is therefore preferable to achieve wider coverage, which can be obtained through the use of enzyme sets like \u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e. However, the GBS method was developed with the objective of reducing genotyping costs, which means that the use of narrow coverage enzyme set may be more appropriate in cases where limited resources are available. As expected, the number of SNPs and genome coverage are dependent on the read number, but only to a certain extent. At some point, it is no longer possible to cover any additional regions and only the quality of the genotypes is improving due to increasing depth of sequencing. In our experimental data, we clearly observed this phenomenon in the \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e enzyme set, which can be considered an additional point for cost optimization. Upon further investigation of SNP metrics, it becomes evident that they also differ in other parameters, such as gene localization and the occurrence of SNPs in WGS (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). On the one hand, \u003cem\u003eApeKI\u003c/em\u003e and \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e libraries demonstrate effective enrichment of gene sequences (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eD). On the other hand, \u003cem\u003eApeKI\u003c/em\u003e library shows a smaller intersection with WGS than other enzyme sets, which may decrease our confidence in identified SNPs. If our main interest relies in LD between SNPs markers and QTLs, as it is the key for GWAS and development of molecular markers, enrichment of SNPs in gene sequences is not of upmost importance.\u003c/p\u003e \u003cp\u003eAn additional point to consider is the fact that three evaluated enzyme sets result in often non-overlapping sets of SNPs (Additional file 5, Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e). This means that they cover different regions of the genome and data from distinct accessions obtained using different enzyme sets cannot be concatenated.\u003c/p\u003e \u003cp\u003eChoice of aligner\u003c/p\u003e \u003cp\u003eAlignment is the intermediate but important step in the SNP calling process. Accurate read positioning is critical for high-quality SNP calling. It appears that the role of the aligner is frequently underestimated. Bowtie2 or BWA are frequently employed in a variety of pipelines without undergoing a comprehensive evaluation. Our results suggest that the choice of aligner has a limited effect on SNP quantity, and their quality and performance is similar. Nevertheless, we notice that the Strobealign aligner, while maintaining similar performance, is several times faster than other aligners tested. This leads to a better use of computational resources and thus to a cost reduction.\u003c/p\u003e \u003cp\u003eChoice of caller\u003c/p\u003e \u003cp\u003eThe role of caller is probably the most crucial in the pipeline, as it provides the final output for many other analyses. Several callers that have been created recently or updated regularly, which have been used in human genetics. However, they have not been widely used in plant genetics, particularly in the context of GBS. In this study, we compared six popular callers and the Stacks pipeline which has been designed for GBS. We also used both classical callers and deep learning based caller.\u003c/p\u003e \u003cp\u003eTo perform analyses such as GWAS, it is necessary to deal with a reliable and consistent set of SNPs. In order to assess this, we called SNPs from the same read sets using different callers. The results demonstrated that the sets of SNPs obtained with using the different callers varied in total number of SNPs as well as their occurrence in the WGS-based SNP set. These results demonstrate that callers may discover SNPs from the same data in a different way, which affects the final result. This can be explained by the fundamentals of the SNP calling process: the caller should efficiently process the aligned reads, which may contain errors, and should be able to take into account genomes with different characteristics. To take into account sources of errors and genome parameters, the callers apply different comprehensive statistical models. However, it is not yet perfect and we again illustrate here that callers may produce SNPs that are unique to a particular caller.\u003c/p\u003e \u003cp\u003eThe ability of deep learning to capture underlying unnoticed dependencies is widely used in data science and has been shown to be applicable to variant calling. Despite the low number of SNPs called by DeepVariant, they can be found in WGS with a high probability. Because identification of SNPs from WGS data can be also affected by the caller effects, we calculated precision, recall and FDR based only on comparing GBS data. In this instance, we again observe that DeepVariant outperforms other callers in terms of precision and FDR, but recall is not so good. With regard to genotyping accuracy, DeepVariant also exhibits the highest R^2 dosage and lower variation in pairwise comparison.\u003c/p\u003e \u003cp\u003eOur analyses raise the question of which caller is preferable. The answer of course depends on the purpose, species and resources available. In some cases, it may be reasonable to obtain the most reliable set in order to be confident in the SNP set and not to work with FP SNPs that do not exist in reality. However, in this way we may miss some valuable SNPs due to the stricter caller procedure. We postulate that the described workflow can be easily adapted to other plant species, with the aim of identifying an optimal combination of molecular experiments and bioinformatics pipelines for a modest cost, as a preliminary step before large-scale analyses.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eThe choice of enzyme set and the design of SNP calling pipelines represent are both crucial aspects of GBS genotyping. However, they are poorly discussed especially in addressing possible optimum of combining both aspects. Our results show that the choice of enzymes significantly influences the number of SNPs identified, gene localisation preferences and accuracy. Furthermore, the performance of SNP callers varied significantly in terms of SNP number, precision, recall and false discovery rate (FDR). These results highlight the importance of optimising both enzyme selection for sequencing libraries and data analysis pipelines to ensure robust and reproducible results.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003eEthics approval and consent to participate\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003eConsent for publication\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e\n\u003cp\u003eAvailability of data and materials\u003c/p\u003e\n\u003cp\u003eThe datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.\u003c/p\u003e\n\u003cp\u003eCompeting interests\u003c/p\u003e\n\u003cp\u003eThe authors declare that they have no competing interests\u003c/p\u003e\n\u003cp\u003eFunding\u003c/p\u003e\n\u003cp\u003eThe authors declare that no funds, grants, or other support were received during the preparation of this manuscript.\u003c/p\u003e\n\u003cp\u003eAuthors\u0026apos; contributions\u003c/p\u003e\n\u003cp\u003eAZ performed the analysis and drafted the manuscript. SB prepared sequencing libraries. CB guided the research. LG guided the research and revised the manuscript. All authors read and approved the final manuscript.\u003c/p\u003e\n\u003cp\u003eAcknowledgements\u003c/p\u003e\n\u003cp\u003eNot applicable\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eThomson MJ (2014) High-Throughput SNP Genotyping to Accelerate Crop Improvement. Plant Breed Biotechnol 2:195\u0026ndash;212\u003c/li\u003e\n\u003cli\u003eMoragues M, Comadran J, Waugh R, Milne I, Flavell AJ, Russell JR (2010) Effects of ascertainment bias and marker number on estimations of barley diversity from high-throughput SNP genotype data. Theoretical and Applied Genetics 120:1525\u0026ndash;1534\u003c/li\u003e\n\u003cli\u003eTranchant‐Dubreuil C, Rouard M, Sabot F (2019) Plant Pangenome: Impacts on Phenotypes and Evolution. In: Annual Plant Reviews online. Wiley, pp 453\u0026ndash;478\u003c/li\u003e\n\u003cli\u003eWarr A, Robert C, Hume D, Archibald A, Deeb N, Watson M (2015) Exome sequencing: Current and future perspectives. G3: Genes, Genomes, Genetics 5:1543\u0026ndash;1550\u003c/li\u003e\n\u003cli\u003eBaird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA (2008) Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One. https://doi.org/10.1371/journal.pone.0003376\u003c/li\u003e\n\u003cli\u003eElshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, Mitchell SE (2011) A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One. https://doi.org/10.1371/journal.pone.0019379\u003c/li\u003e\n\u003cli\u003eHeffelfinger C, Fragoso CA, Moreno MA, Overton JD, Mottinger JP, Zhao H, Tohme J, Dellaporta SL (2014) Flexible and scalable genotyping-by-sequencing strategies for population studies. BMC Genomics. https://doi.org/10.1186/1471-2164-15-979\u003c/li\u003e\n\u003cli\u003eShin MG, Ithnin M, Vu WT, Kamaruddin K, Chin TN, Yaakub Z, Chang PL, Sritharan K, Nuzhdin S, Singh R (2021) Association mapping analysis of oil palm interspecific hybrid populations and predicting phenotypic values via machine learning algorithms. Plant Breeding 140:1150\u0026ndash;1165\u003c/li\u003e\n\u003cli\u003eKoboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, Miller CA, Mardis ER, Ding L, Wilson RK (2012) VarScan 2: Somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res 22:568\u0026ndash;576\u003c/li\u003e\n\u003cli\u003eGarrison E, Marth G (2012) Haplotype-based variant detection from short-read sequencing. \u003c/li\u003e\n\u003cli\u003ePoplin R, Chang P-C, Alexander D, et al (2018) A universal SNP and small-indel variant caller using deep neural networks. Nat Biotechnol 36:983\u0026ndash;987\u003c/li\u003e\n\u003cli\u003eCatchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA (2013) Stacks: an analysis tool set for population genomics. Mol Ecol 22:3124\u0026ndash;3140\u003c/li\u003e\n\u003cli\u003eBushnell B, Rood J, Singer E (2017) BBMerge \u0026ndash; Accurate paired shotgun read merging via overlap. PLoS One 12:e0185056\u003c/li\u003e\n\u003cli\u003eDanecek P, Bonfield JK, Liddle J, et al (2021) Twelve years of SAMtools and BCFtools. Gigascience. https://doi.org/10.1093/gigascience/giab008\u003c/li\u003e\n\u003cli\u003ePoplin R, Ruano-Rubio V, DePristo MA, et al (2018) Scaling accurate genetic variant discovery to tens of thousands of samples. bioRxiv 201178\u003c/li\u003e\n\u003cli\u003eTaniguti CH, Taniguti LM, Amadeu RR, et al (2023) Developing best practices for genotyping-by-sequencing analysis in the construction of linkage maps. Gigascience. https://doi.org/10.1093/gigascience/giad092\u003c/li\u003e\n\u003cli\u003eGupta SK, Baek J, Carrasquilla-Garcia N, Penmetsa RV (2015) Genome-wide polymorphism detection in peanut using next-generation restriction-site-associated DNA (RAD) sequencing. Molecular Breeding. https://doi.org/10.1007/s11032-015-0343-0\u003c/li\u003e\n\u003cli\u003ePoland J (2011) PstI-\u0026shy;MspI GBS Genotyping-\u0026shy;by-\u0026shy;sequencing Protocol PstI-\u0026shy;MspI. \u003c/li\u003e\n\u003cli\u003eAndrews S. (2010) FastQC: a quality control tool for high throughput sequence data. \u003c/li\u003e\n\u003cli\u003eLi H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. \u003c/li\u003e\n\u003cli\u003eWu X, Heffelfinger C, Zhao H, Dellaporta SL (2019) Benchmarking variant identification tools for plant diversity discovery. BMC Genomics. https://doi.org/10.1186/s12864-019-6057-7\u003c/li\u003e\n\u003cli\u003eYun T, Li H, Chang PC, Lin MF, Carroll A, McLean CY (2020) Accurate, scalable cohort variant calls using DeepVariant and GLnexus. Bioinformatics 36:5582\u0026ndash;5589\u003c/li\u003e\n\u003cli\u003eLangmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9:357\u0026ndash;359\u003c/li\u003e\n\u003cli\u003eBushnell B (2014) BBMap: a fast, accurate, splice-aware aligner. \u003c/li\u003e\n\u003cli\u003eSahlin K (2022) Strobealign: flexible seed size enables ultra-fast and accurate read alignment. Genome Biol 23:260\u003c/li\u003e\n\u003cli\u003eQuinlan AR, Hall IM (2010) BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26:841\u0026ndash;842\u003c/li\u003e\n\u003cli\u003eR Core Team (2021) R: A Language and Environment for Statistical Computing. \u003c/li\u003e\n\u003cli\u003eWickham H, Averick M, Bryan J, et al (2019) Welcome to the Tidyverse. J Open Source Softw 4:1686\u003c/li\u003e\n\u003cli\u003eJohn Fox and Sanford Weisberg (2019) An R Companion to Applied Regression, third. Sage\u003c/li\u003e\n\u003cli\u003eAlboukadel Kassambara (2023) rstatix: Pipe-Friendly Framework for Basic Statistical Tests. \u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"plant-methods","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"plme","sideBox":"Learn more about [Plant Methods](http://plantmethods.biomedcentral.com/)","snPcode":"13007","submissionUrl":"https://submission.nature.com/new-submission/13007/3","title":"Plant Methods","twitterHandle":"@PlantMethods","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"Genotype-by-sequencing (GBS), ddRAD, SNP calling, soybean, comparison","lastPublishedDoi":"10.21203/rs.3.rs-5821852/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-5821852/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003eGenotype-by-sequencing (GBS) is a cost-effective method for large-scale genotyping, widely used across various species, particularly those with large genomes. A critical aspect of GBS lies in the selection of restriction enzymes for genome digestion and the optimization of data analysis pipelines. However, few studies have comprehensively examined the combined effects of enzyme choice and pipeline configuration.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eIn this study, we created GBS libraries using three commonly used restriction enzyme combinations (\u003cem\u003eHindIII\u003c/em\u003e-\u003cem\u003eNlaIII\u003c/em\u003e, \u003cem\u003ePstI\u003c/em\u003e-\u003cem\u003eMspI\u003c/em\u003e, and \u003cem\u003eApeKI\u003c/em\u003e) and assessed multiple SNP-calling pipelines in 15 soybean varieties. We tested four aligners (BWA-MEM, Bowtie2, BBMap, and Strobealign) and seven SNP callers (Bcftools, Stacks, DeepVariant, FreeBayes, VarScan, BBCallVariants, and GATK). Our finding reveal that enzyme choice significantly influences the number of identified SNP, gene localization preferences, and accuracy. Furthermore, the performance of SNP callers varied markedly in terms of SNP count, precision, recall, and false discovery rate (FDR). DeepVariant exhibited the highest accuracy, with 76.0% of its SNPs intersecting with whole-genome sequencing (WGS)-derived SNPs and an FDR of 0.0095, compared to FreeBayes, which had 47.8% intersection and an FDR of 0.6321.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eOur findings underscore the importance of optimizing both enzyme selection for sequencing libraries and data analysis pipelines to ensure robust and reproducible results. This study provides a general framework for designing large-scale genotyping experiments aimed to specific quality and quantity requirements in various plant species.\u003c/p\u003e","manuscriptTitle":"The evaluation of different combinations of enzyme set, aligner and caller in GBS sequencing of soybean","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-01-20 09:15:15","doi":"10.21203/rs.3.rs-5821852/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-04-29T22:17:29+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-04-27T12:26:14+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"100373281184177043732781939010973069225","date":"2025-04-15T10:27:13+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-02-17T09:13:37+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"256875930723664779040035886707017804709","date":"2025-02-10T09:15:39+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-02-06T02:49:07+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-01-17T11:34:54+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-01-17T11:31:57+00:00","index":"","fulltext":""},{"type":"submitted","content":"Plant Methods","date":"2025-01-13T17:19:03+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"plant-methods","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"plme","sideBox":"Learn more about [Plant Methods](http://plantmethods.biomedcentral.com/)","snPcode":"13007","submissionUrl":"https://submission.nature.com/new-submission/13007/3","title":"Plant Methods","twitterHandle":"@PlantMethods","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"7ff06b01-89b5-491f-97f2-a4c4888f1210","owner":[],"postedDate":"January 20th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2025-08-11T16:00:19+00:00","versionOfRecord":{"articleIdentity":"rs-5821852","link":"https://doi.org/10.1186/s13007-025-01410-8","journal":{"identity":"plant-methods","isVorOnly":false,"title":"Plant Methods"},"publishedOn":"2025-08-06 15:57:18","publishedOnDateReadable":"August 6th, 2025"},"versionCreatedAt":"2025-01-20 09:15:15","video":"","vorDoi":"10.1186/s13007-025-01410-8","vorDoiUrl":"https://doi.org/10.1186/s13007-025-01410-8","workflowStages":[]},"version":"v1","identity":"rs-5821852","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-5821852","identity":"rs-5821852","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}