A comparison of BeadChip and WGS genotyping outputs using partial validation by sanger sequencing

Background Head-to-head comparison of BeadChip and WGS/WES genotyping techniques for their precision is far from straightforward. A tool for validation of high-throughput genotyping calls such as Sanger sequencing is neither scalable nor practical for large-scale DNA processing. Here we report a cross-validation analysis of genotyping calls obtained via Illumina GSA BeadChip and WGS (Illumina HiSeq X Ten) techniques. Results When compared to each other, the average precision and accuracy of BeadChip and WGS genotyping techniques exceeded 0.991 and 0.997, respectively. The average fraction of discordant variants for both platforms was found to be 0.639%. A sliding window approach was utilized to explore genomic regions not exceeding 500 bp encompassing a maximal amount of discordant variants for further validation by Sanger sequencing. Notably, 12 variants out of 26 located within eight identified regions were consistently discordant in related calls made by WGS and BeadChip. When Sanger sequenced, a total of 16 of these genotypes were successfully resolved, indicating that a precision of WGS and BeadChip genotyping for this genotype subset was at 0.81 and 0.5, respectively, with accuracy values of 0.87 and 0.61. Conclusions We conclude that WGS genotype calling exhibits higher overall precision within the selected variety of discordantly genotyped variants, though the amount of validated variants remained insufficient.

non-NGS (next-generation sequencing) techniques for genotype calling limits evaluations of accuracy.
BeadChip genotyping is an efficient and scalable way of genotype resolution, with two inherent limitations: a necessity to decompose binary alleles and confinement to a predefined list of genotyped variants. These two limitations, however, do not prevent its usefulness for a variety of clinical and non-clinical applications [3]. Predefined nature of the variants to be tested makes Bead-Chip genotyping amenable to validation by either PCR (polymerase chain reaction), or Sanger sequencing which has been recently shown to have limited utility and erroneous behaviour in the validation of NGS variants [4]. A comparison of genotyping calls made by BeadChip and WGS/WES techniques may provide an insight into the possible nature of the discordant calls observed during genotyping quality control stage. Here we attempted to identify analytical issues leading to the discordance in genotyping calls made independently by WGS and Bead-Chip techniques. Table 1 lists the sequencing statistics for the three sequenced samples. FastQC reports are available as Additional files 1, 2, 3, 4, 5 and 6 for sample_001, sample_ 002, sample_003, respectively. Percentage of reads falling into a category with averaged Phred scaled sequencing quality above 30 is shown by %Q. The absence of unpaired reads, repeatability in terms of GC content and more than 90% of bases exceeding sequencing quality of 30 was used as a mark of confidence in data quality.

Mapping statistics
All the data produced by WGS were analyzed for their depth (DOC) and breadth (BOC) of coverage using GATK 3.8 DepthOfCoverage tool [5]. The mean filtered coverage for all three samples exceeded 27x (Fig. 1a-c, Table 2), which complied with recommendations [6]. Repeatability in BOC values for each base quality interval (Fig. 1b, Table 2) and other sequencing metrics (Table 1) proved the quality of the sequencing.

Concordance metrics
Percentages of discordant calls per chromosome for each sample are shown in Fig. 2. The average fraction of discordant results among all comparable genotyping results (intersection of BeadChip and WGS obtained genotyping data) was at 0.3317% for sample_001 (male), 0.8448% for sample_002 (female) and 0.7392% for sample_003 (male) which falls within the range reported previously [2]. Both the sequencing quality and mapping quality parameters were reproduced for all three sequenced samples with an average genotype concordance estimate being higher than 99%.

Mapping of analyzed variants
For each pair of BeadChip-genotyped neighboring variants, distance intervals were extracted within both concordant and discordant group to generate pairs of distance values before and after each variant, followed by their visualization. The resulting maps with a Gaussian kernel density estimation are shown in Fig. 3. The observed approximate evenness of distribution for both concordant and discordant variants supported by the visible clusters allocation along the bisectrix of the axes shows that both concordant and discordant variants are evenly distributed across the chromosome length, with no congregation across the genome. The only observed difference in cluster allocation arises from the frequency of observation for the variants from each group. In other words, less frequent discordant group corresponds to lower mapping density and, consequently, to more considerable distances between every two neighboring variants.
Overall randomness of the locations of discordant genotypes across the genome, measured as cluster evenness and cluster center distance from the bisectrix, was high. However, because locations of discordant variants were limited to variants present on the corresponding BeadChip, the comparison with WGS was necessarily limited to the locations of BeadChip genome variants only. For this reason, patterns in discordant genotypes location throughout the genome were only detectable at locations thoroughly covered by variants presented at BeadChip, which generally are outside any complex

Concordance analysis
Calculated confusion matrices for all three analyzed DNA samples are shown in Fig. 4, with metrics values presented in Table 3. These values show high overall concordance in genotyping between both WGS and BeadChip techniques. Detected discordant results may have arisen from ambiguous BeadChip genotyping call within multiallelic sites. Nevertheless, all calculated quality metrics values surpassed the recommended thresholds [6].

Genotyping quality metrics distributions
For all three samples, the distributions of the WGS genotyping metrics were analyzed and compared within the groups formed based on the concordance, variation class and genotype zygosity criteria (data not shown). While differences in mean values were statistically significant (Welch's t-test, 0.05 p-value threshold) in several groups, distributions themselves exhibited low intergroup divergence. Moreover, for all group comparisons, observed differences were quite small, and their patterns were not uniform, with several groups being insufficiently large, potentially obscuring inter-group dissimilarity. Thus, discordance in genotyping could not be explained by the variance in these quality metrics. When quality metrics were analyzed for the BeadChip calls, some differences in parameters distributions were found. The distinction was remarkable only for the SNV group, as it included the significant amount of all genotyping results and thus could more precisely represent any actual deviations. Inspection of R, Theta and GC scores for both concordant and discordant variants revealed a pattern of discordant variants located close to the borders of the variant clusters (Fig. 5). Importantly, similar clusters are used in the Illumina GenomeStudio software to assign genotypes to variants (AA, AB and BB), and, therefore, any genotype lying far from the cluster centre may be mistakenly assigned an irrelevant designation. This explanation for the observed discordance pattern was, however, relevant only to sample_002 and sample_ 003. In sample_001, distribution of discordant variants within the clusters was more or less even. Therefore, we conclude that the observed phenomenon requires further examination. Potential finding of the clustering pattern for discordant calls may then be exploited for quality control.

Sanger sequencing
The sliding window approach performed on sample_002 resulted in the mapping of 6 regions containing from 1 to 3 variants successfully genotyped on both WGS and BeadChip platforms and exhibiting discordant genotypes between the used platforms. Locations of these regions and encompassed variants are collated in Tables 4 and 5.
A total of 12 discordant variants were selected for validation by Sanger sequencing. These variants were accompanied by 14 concordant variants located closely to discordant variants. 1. Selected SNPs were successfully genotyped using both WGS and BeadChip platforms.  2. Selected SNPs are located close to each other within 500 bp window length (reasonable limit of one Sanger read). 3. Specific primers can be selected for these SNPscontaining region (by Primer-BLAST [7]).
The list of designed primers with respective amplification parameters can be found in Table 5. The gDNA of sample_002 was used for amplifying regions of interest (ROIs) in Table 4. All ROIs, except those located on chromosomes 10 and 13, were successfully amplified and Sanger sequenced. The comparison of genotypes obtained via microarray genotyping, whole genome sequencing and Sanger sequencing of the amplified ROIs is shown in Table 6. Sanger-derived genotypes containing only one letter were the ones obtained from only one read. Plus signs denote Sanger genotypes concordant with the WGS-derived genotype, while double plus signs denote concordance of all three genotyping methods. Hashtag represents concordance with the BeadChip-derived genotype, and asteriskan absence of concordance with any of the tested methods. Sixteen Sanger-resolved diploid genotypes (forward and reverse Sanger chromatograms covered the variant location) out of the 26 listed in Table 6 were used for confusion matrices calculation using Sanger-derived calls as a "truth" call set, which resulted in WGS and BeadChip precision of 0.81 and 0.5, respectively, and the accuracy values of 0.87 and 0.61. Results of Sanger validation may possibly be explained by low complexity or repetitive genomic context which surrounds some of the validated variants, hindering accuracy and precision of either genotyping or read alignment. All the listed variants which were discordant between BeadChip and WGS were plotted (Fig. 5) to reveal possible clustering of discordant variants based on the initial BeadChip genotyping metrics. No clustering of the validated variants was observed. Thus, current analysis did not allow to make any definite conclusions about the clustering of discordant and concordant variants based on their genotyping quality metrics.

Discussion
Although genotype concordance analysis experiments using different sequencing and genotyping platforms have  been performed previously [2,8], no reports on attempts to explain the observed differences were made. Here we tried to scrutinize the underlying genotyping process proxy such as genotyping quality metrics to find a possible explanation for the discordance pattern. Unfortunately, the genotyping discordance of the observed levels (less than 1%) is usually an underestimation, and does not motivate investigations into the nature of the discordant genotype calling. However, we show that discordant genotypes tend to form clusters in a 2-dimensional space, where the most vivid dimensions are the genotyping quality and technical metrics obtained from the BeadChip genotyping pipeline. Therefore, speculation on possibly lower precision of the BeadChip genotyping platform, as compared to WGS-based pipelines, might find a new conceivable basis upon further investigation. Genotype assignment in the array pipelines is based on marker clustering in a 2-dimensional space (clusters A/A, A/B, B/B), which might happen to be erroneous due to poor cluster separation. This clusterization problem may possibly be solved by incorporating additional dimensions into the analyzed genotyping metrics space, which can be exploited for enhancement of the BeadChip genotyping pipeline upon further investigation and vast Sanger validation of the variants.

Conclusions
Here we show the presence of some parametric differences in quality metrics of genotyping performed by WGS and BeadChip. This phenomenon warrants comprehensive investigation by combining genotyping metrics produced by WGS and BeadChip pipelines and extracting patterns in the observed discordance. In clusters, Sanger validation should be performed for genotype resolution.

Materials
Three human genomic DNA samples, two males and one female, were selected for this comparison. After collection by Oragene DNA saliva-based collection device (DNA Genotek, Canada), genomic DNA was extracted as per the manufacturer protocol and stored frozen in Tris-EDTA buffer at − 20°C until genotyping was performed.

BeadChip genotyping
Infinium iSelect 24 × 1 HTS Custom Beadchip Kit (GSA-sharedCUSTOM_20018389_A2) genotyping was performed using 50 ng of genomic DNA. Microarray fluorescence was scanned using the Illumina iScan system, and the genotype calling was executed using the Illumina GenomeStudio Genotyping Module software (Illumina, USA). No imputation was implemented, and only variants successfully genotyped on the specified array were used in the analysis.

Quality metrics and validation statistics
Concordance analysis implied calculation of specificity, sensitivity, precision, accuracy, genotype concordance, non-reference genotype concordance and non-reference sensitivity of genotyping between BeadChip and WGS. The analysis implemented confusion matrices calculation (Fig. 6). As there is no baseline "truth" when a WGS call set is compared with a BeadChip one, no method can be described as "comparator". Because of this, each call set within each pair was treated as alternating "truth" or "test" call set, followed by the averaging of obtained statistics. Filling those matrices with observed numbers of class counts was followed by dimensionality reduction, which implied leaving one class in both call sets intact and combining the counts in all other classes (Fig. 6). The above-specified quality metrics for the initial and reduced matrices were calculated, as shown in Fig. 7 [2] with amendments. As each non-reduced matrix produced 6 submatrices, the calculated sensitivity and precision value for each submatrix was weighted by the fraction of analyzed elements in a nonreduced matrix (i.e., each calculated metric was then multiplied by the ratio of orange-outlined elements to all elements in the reduced matrix in Fig. 7). Such an approach was only necessary for precision and sensitivity metrics due to the fact that these metrics were calculated on non-overlapping subsets of the initial matrix and their values depended on the size of those subsets, given uneven marker distribution in the initial matrix. Other used metrics (accuracy, specificity) were calculated as mean values between 6 produced metric values (excluding NANs, which originated from calculations involving division by zero). All concordant and discordant genotyped variants were analyzed for genotyping quality metrics provided in the VCF files and GenomeStudio report files for WGS  TRUTHa call set produced by an orthogonal method (comparator), TEST a call set produced by a test method equal to the AB cluster's theta mean, or 1 if it is equal to or greater than the BB cluster's theta mean. B Allele Freq is linearly interpolated between 0 and 1.

Variant selection for sanger sequencing
An intersection of WGS-and array-genotyped markers which exhibited discordant calls between the two platforms was used for selection. A sliding window approach was utilized to find regions spanning no more than 500 bp (Sanger sequencing reasonable read length limit) and encompassing as many discordant variants from the selected set as possible. All primers for PCR amplification were designed using NCBI Primer-BLAST suite [7] with default parameters (except the melting temperature limits of 58-62°C) and the human genome reference sequence for BLAST. PCR product lengths and primer lengths were manually optimized to find the most suitable unique match. Primer synthesis, amplification and Sanger sequencing was performed by Evrogen (Russia, Moscow). Sanger chromatograms were visualized and analyzed using 4Peaks software (Nucleobytes, The Netherlands) and CodonCode Aligner (CodonCode Corp., USA) with default trimming and quality filtering parameters.