extfieldsofstudy:
Computer Science
Medicine

text:
Establishing analytical validity of BeadChip array genotype data by comparison to whole-genome sequence and standard benchmark datasets

Background Clinical use of genotype data requires high positive predictive value (PPV) and thorough understanding of the genotyping platform characteristics. BeadChip arrays, such as the Global Screening Array (GSA), potentially offer a high-throughput, low-cost clinical screen for known variants. We hypothesize that quality assessment and comparison to whole-genome sequence and benchmark data establish the analytical validity of GSA genotyping. Methods To test this hypothesis, we selected 263 samples from Coriell, generated GSA genotypes in triplicate, generated whole genome sequence (rWGS) genotypes, assessed the quality of each set of genotypes, and compared each set of genotypes to each other and to the 1000 Genomes Phase 3 (1KG) genotypes, a performance benchmark. For 59 genes (MAP59), we also performed theoretical and empirical evaluation of variants deemed medically actionable predispositions. Results Quality analyses detected sample contamination and increased assay failure along the chip margins. Comparison to benchmark data demonstrated that > 82% of the GSA assays had a PPV of 1. GSA assays targeting transitions, genomic regions of high complexity, and common variants performed better than those targeting transversions, regions of low complexity, and rare variants. Comparison of GSA data to rWGS and 1KG data showed > 99% performance across all measured parameters. Consistent with predictions from prior studies, the GSA detection of variation within the MAP59 genes was 3/261. Conclusion We establish the analytical validity of GSA assays using quality analytics and comparison to benchmark and rWGS data. GSA assays meet the standards of a clinical screen although assays interrogating rare variants, transversions, and variants within low-complexity regions require careful evaluation. Supplementary Information The online version contains supplementary material available at 10.1186/s12920-022-01199-8.

Aim and design of study
This study defines an analytical validation framework for detecting and limiting genotyping error in GSA data (Fig. 1). To minimize platform specific genotyping biases, internally generated genotype data from independent platforms were paired and compared with publicly available genotype datasets.

Samples and datasets
To generate a reference genotype cluster file for the GSA, 664 DNA samples were purchased from the Coriell Institute for Medical Research, Camden, NJ (https:// www. corie ll. org) and 460 samples were selected from the Sanford Biobank. Individuals with biobank samples were enrolled in protocol number 03-11-061 approved by the Sanford Research Institutional Review Board. These samples were selected to cover different ethnicities (14 Coriell diversity panels, Additional file 1: Table S2) and the technical variability of the DNA extraction methods (460 samples from the Sanford Biobank). To capture the technical variability of the Infinium ® HTS Assay protocol (Illumina Inc.), all samples were genotyped in triplicate (by different technicians, robot-instrument configurations, reagent lots, and days) using the Infinium Global Screening Array-24 v.1.0 BeadChip. The resulting data were loaded into GenomeStudio v2.0.2 and used to generate the genotype cluster files per manufacturer recommendations (https:// www. illum ina. com/ Docum ents/ produ cts/ techn otes/ techn ote_ infin ium_ genot yping_ data_ analy sis. pdf ). Of the 1104 samples used in cluster file generation 72 were also included among the 263 samples used to define analytical validity (Supplementary Sects. "Methods" and "Aim and design of study"). Two hundred sixty-three DNA samples from Coriell were selected as representative of individuals from the 1000 Genomes Project Consortium (n = 258) and from the Genome in a Bottle Consortium (GiAB) [33] (n = 5) (Additional file 1: Table S1). Additionally, they were selected to assess assays genotyping alleles with ≥ 1% minor allele frequency (MAF) in the general population (Additional file 1 Table S3). These 263 DNA samples were resequenced with whole genome sequencing (rWGS) and genotyped in triplicate (263 × 3) with the GSA. These data were compared to 1KG and to publicly available Whole Genome Sequence (pWGS) data (1KG phase 3; downloaded: June 2018). This defined 4 genotype datasets for the 263 samples: (i) triplicate GSA genotypes (ii) pWGS, (iii) rWGS, and (iv) 1KG (Additional file 1: Table S4, Fig. S2, Sect. S3). All analyses including mapping, alignment, and genotyping were performed using HumanG1Kv37 (Genome Reference Consortium Human build 37).

Data processing GSA quality control (QC)
Laboratory QC Genotype clusters for the variants used for clinical reporting were manually curated to ensure accurate variant calling. Other variants were automatically curated using Illumina-recommended filters (https:// www. illum ina. com/ Docum ents/ produ cts/ techn otes/ techn ote_ infin ium_ genot yping_ data_ analy sis. pdf ). Using the data of DNA samples from 1104 individuals run on the GSA in triplicate, the cluster file analyses of each GSA assay found that 610,771 (92%) assays passed and 50,355 (8%) assays failed clustering quality control. Those that failed were excluded and marked as no-calls (./.) in the VCF files.

Data comparisons Principal component analysis
Principal component analysis (PCA) was used to test for intact super-population structure as a corollary for absence of batch and technical artifacts in the genotyping datasets. PCA structure derived from GSA data was compared to the super-population structure derived from 1KG data. Aggregate quality control analysis of the GSA data. A Principal Component Analysis (PCA) plots of 1KG data and GSA genotype data. red: African (AFR), yellow-green: Admixed Americans (AMR), dark-green: East Asian (EAS), blue: European (EUR), purple: South Asian (SAS). B Heatmaps of BeadChip array quality control analysis of call-rate (left), p10GC (middle), and estimated DNA contamination (right). Color gradient scales for the three panels are as follows: call-rate (orange < 0.94-blue > 0.99), p10GC (yellow < 0.50-blue > 0.60) and estimated DNA contamination (rainbow gradient: purple ~ 1%, blue ~ 2%, green ~ 3%, orange/red ~ >4%). C Heatmaps of reproducibility quality control analysis using replicate data as measured by call rate, estimated DNA contamination, number of assays with no genotype calls, and heterozygote to homozygote ratio. Color gradient scales for these four heatmaps are as follows: No genotype calls (blue < 166,000-orange > 400,000), and rainbow gradient for call rate (purple > 0.99-red < 0.94), estimated DNA contamination (purple < 1%-red > 4%), and heterozygote/homozygote ratio (purple > 2.25-red < 1.25), respectively

Whole genome sequence data quality control (QC) Bioinformatics QC
For bioinformatics quality control of rWGS data (n = 263), central tendency and anomalous outlier data points were assessed for (i) total processed reads, (ii) discordant reads, (iii) mapq0 reads, (iv) unmapped reads, (v) mapped reads, and (vi) average depth of sequencing (Additional file 1: Tables S4 and S5). On average > 95% of processed reads per sample (731,227,993/767,540,183 reads) mapped to the reference sequence. Because the concordance of two rWGS datasets (HG00111 and HG00257) with the 1KG data were 0.870 and 0.622, they were dropped from our GSA analyses leaving a total of 261 samples in the rWGS dataset. Variation data for 260 samples (SNVs and short indels ≥ 20× coverage and a Phred score ≥ 30) were submitted to the dbSNP database (Additional file 1: Sect. S8).

Performance metrics Genotype concordance, sensitivity, specificity and positive predictive value (PPV)
GSA and rWGS genotypes were compared to each other and to 1KG genotypes using the following performance metrics: (i) genotype concordance (C), (ii) sensitivity (S), (iii) specificity (P), and (iv) positive predictive value (PPV). We used the following definitions of genotype classification to label genotypes as positive [true positive (tp), false positive (fp)], negative [true negative (tn), false negative (fn)], or discordant (x) ( Table 1): Given the above definitions of true/false positive and negative and discordant genotypes (Table 1), we computed the performance metrics as follows: Genotype concordance (C) Genomic complexity of variation locus (low-complexity regions) To categorize SNVs based on the genomic complexity of the GSA assay locus, we used the UCSC genome browser bed-file definitions to define simple-repeats, micro-satellite regions, and low-complexity regions. The SimpRep, Microsatellites, and RepeatMasker bedfiles were downloaded from the UCSC Genome Browser FTP site and intersected with the GSA manifest file. Across the HumanG1Kv37 reference sequence, there were 962,715 simple repeat, 41,573 microsatellite, and 5,298,131 RepeatMasker regions.

GSA panels Medically actionable predispositions (MAP) 59 gene panel
GSA assays targeting potentially disease-associated variants in MAP59 genes [35] were selected in a multistep process (Table 2). Firstly, GSA assays that interrogated positions within 1000-bases upstream and downstream of the transcript start and end in HumanG1Kv37 were selected for the RefSeq transcript chosen for each gene. Secondly, alleles were annotated with their respective ClinVar classifications, and those that had at least one classification of pathogenic or likely pathogenic were selected. Thirdly, these assays were curated by clinical and laboratory staff to define a managed variant list (MVL) of 1883 assays appropriate for clinical reporting.

Statistics and compute infrastructure
Statistical analyses and data visualization were performed using R (version 3.4.3). Data analysis was done on a Linux Operating System with the following configuration: x86_64, 32 CPUs, 2.8 GHz AMD Opteron Processor 6320. AWS EC2 instances were spun-up for large compute jobs. All NGS and GSA data were archived on AWS S3. In-house software and data processing code and scripts were written primarily in Perl, Ruby, awk, and bash.

Data summary
DNA samples from 263 individuals were purchased from Coriell and genotyped in triplicate (n = 789) with the GSA. Genotypes and data for each replicate were saved to a VCF file. The GSA data were grouped and summarized as replicate datasets 1, 2, and 3. Of the 263 samples, 258 were present in the 1KG. Of the other 5 samples, 3 were from the Personal Genomes Project (PGP) [36], and 2 were from the NIGMS Human Genetic Cell Repository. The 263 × 3 data were compared with the 1KG data and with two WGS datasets, the resequenced WGS data (n = 261; rWGS = 37×) and the downloaded public WGS data (n = 24; pWGS = 51×) (Additional file 1: Sect. S3).

Principal component analysis defines the same population structure in GSA data and 1KG data
Principal component analysis (PCA) on each replicate of autosomal GSA data identified 5 major super populations conserved across replicates. PCA of the 1KG autosomal genotype data from the same loci generated a similar population structure ( Fig. 2A). This suggested that the GSA data did not have confounding technical factors skewing the PCA plot. To determine if fewer GSA genotypes were sufficient for this test, we randomly sub-sampled close to 10,000 genotypes; these recapitulated the population structure (Additional file 1: Sect. S4).

GSA triplicate data analysis shows data reproducibility in the majority of samples and no detectable stochastic QC failure
Given that PCA did not detect major technical confounders within the GSA genotypes, we analyzed the 263 × 3 data for quality and reproducibility [10] (Table 3; Fig. 3). Data were stratified by BeadChip identifiers and sample location on the BeadChip (row, column). Additionally, samples were grouped by replicates, and each replicate sample was evaluated for (i) genotype callrate (n = 610,771 assays), (ii) p10GC, and (iii) estimated DNA sample contamination (Additional file 1: Sect. S5). Aggregate quality control analysis showed a lower p10GC in higher numbered rows on the BeadChip (Fig. 2B); excluding contaminated samples, p10GC ranged from 0.56-0.61 (mean = 0.60, SD = 0.0085) in row 1 and from 0.50-0.61 (mean = 0.55, SD = 0.03) in row 12. Over 99% (782/789) of samples had a call rate of > 0.98. 3 samples in the third replicate dataset were contaminated, and 2 of these 3 samples had a call rate < 0.98 (0.93 and 0.94, Fig. 2C).
To test if call-rates were reproducible across replicates, we measured deviations from expectation and dispersion. The first approach, a Z-score method, computes the number of standard deviations a replicate Table 2 Selection process for GSA assays targeting genotypes considered medically actionable predispositions  sample call-rate is from the expected as defined by the global dataset average and standard deviation. The second approach computes the average call-rate of all replicates for a given sample and then computes variation around the average. Using the Z-score method, 7 samples had a Z-score ≤ − 4. With a more conservative cut-off (Z-score < −3), 11 samples deviated from expectation (Additional file 1: Sect. S6; Fig. S12). When analyzed relative to the BeadChip row and column, outlier Z-scores occurred for wells on the edge of the Illumina BeadChip-R12C01 or R11C01; the only exceptions were two contaminated samples that were in wells R01C01 and R01C02. Dispersion metrics calculated for call rates across each set of three replicates (Table 4) identified higher relative dispersion for the same samples detected by the Z-score method. Replicate pairwise concordance was calculated to assess the stochastic nature of sample genotyping quality and these were plotted as a 3D scatter-plot: [R1 vs. R2 (x-axis), R2 vs. R3 (y-axis), and R1 vs. R3 (z-axis)] (Fig. 3). The data along the diagonal of the cube are correlated data values across triplicates for all measured GSA genotypes for a given DNA sample. 260 of 263 samples in the triplicate dataset (262/263 R1 vs. R2; 260/263 R2 vs. R3; 261/263 R1 vs. R3) had concordance greater than 0.999 between replicates suggesting high reproducibility. Offdiagonal points, i.e., those with poor call rates (< 0.98) (Fig. 2B), were along the edge of Illumina chip or contaminated; we did not observe random occurrence of poor call rates.

Grouping GSA assays by variation type shows that SNVs have > 0.99 performance relative to the benchmark dataset 1KG across all metrics
Of the 263 samples with GSA data, 258 had corresponding 1KG genotype data for computing performance metrics of concordance, sensitivity, specificity, and PPV. Each GSA assay was grouped according to the type of nucleotide change assessed: (a) single nucleotide variant (SNV), (b) multi-allelic variant (MAV), (c) insertion, and (d) deletion (Table 5). SNVs accounted for 99.3% (656,601/661,126); 610,771 of these passed cluster file quality control, and 594,361 detected genotypes present in the 1KG. Among the MAV assays, 526 of 616 passed cluster file QC; however, because only 3 of these had genotypes present in the 1KG, we excluded MAVs from further analysis. Among insertion assays, 1044 of 1110 passed cluster file QC, and 36 of these had genotypes present in the 1KG. Among deletion assays 2677 of 2799 assays passed cluster file QC, and 95 of these had genotypes present in the1KG. Using the three replicate . The data is plotted as correlation across triplicates for all measured GSA genotypes for a given DNA sample. Note that most samples had concordance greater than 0.999 between replicates suggesting high reproducibility. A few samples had off-diagonal points, i.e., those with poor call rates or reproducibility. The color rainbow gradient is from blue (< 0.996) to dark red (1.00) GSA genotype datasets, the performance metrics of SNV assays were > 0.99. In contrast, insertion assays had highly variable concordance with the 1KG, and deletion assays had poor performance metrics (Fig. 4A).

GSA assays for transitions perform better than do those for transversions
Classifying the GSA-detected SNVs as transitions (purine-to-purine OR pyrimidine-to-pyrimidine) or transversions (purine-to-pyrimidine or vice versa)  (Fig. 4B). The assays for transversions between complementary nucleotides (i.e., A>T, T>A, C>G, G>C; Additional file 1: Sect. S7) had lower sensitivity (< 0.99) and lower cluster file QC pass rate (66-73%; Table 6) than did those for other transversions.    . 2011), the average performance metrics for GSA assays passing cluster file QC in each bin showed that PPV and sensitivity suffered when the alternate allele frequency was < 5%, whereas specificity and concordance declined as the alternate allele frequency increased (Fig. 4C). Among the 594,346 GSA SNV assays with IKG genotypes, 476,707 had a PPV equal to 1 (zero false positives) based on concordance with the 1KG and rWGS genotypes (Table 8). We observed that 81% of GSA assays in the [0-0.1%] and 28% of GSA assays in the (0.1-1%] bins had PPV < 1 (Fig. 5; Table 9), whereas other allele frequency bins had an average of 13% (8-17%) with a PPV < 1 (Table 8). These results are consistent with prior observations showing that accurate calling of rare alleles (MAF < 0.01) by genotyping arrays is compromised by low genotype frequencies and an absence of the homozygous alternate alleles needed for construction of cluster files [22,37].

GSA assays interrogating low-complexity genomic regions perform poorer than other assays
To determine assay performance characteristics within repetitive regions of the genome, we intersected GSA assays with annotated low complexity regions (LCRs) including simple repeats, microsatellites, and repeat masked (RepeatMasker-defined) regions in the human genome. Of a total of 594,346 assays passing QC and present in the 1KG, 203,901 (~ 34%) assessed a variant within one of the three annotation classes. 201,579 GSA assays mapped within the RepeatMasker class. Overlapping partially with the other two classes, 431 GSA assays mapped within the simple repeat class. GSA assays targeting genotypes within each LCR class had poorer performance metrics than did assays interrogating genotypes outside of these regions (Fig. 4D).

rWGS performed better than GSA relative to the benchmark dataset 1KG
rWGS data corresponding to GSA assays passing QC were extracted from the rWGS gVCF files and compared to the 1KG. Restricting the analyses to GSA assays for which > 90% of rWGS samples had genotype data defined 602,582 assays and excluded 38,093 GSA assays. An additional 9642 assays on the chromosome X were excluded due to discrepancies in genotype representation in comparison datasets. For the remaining 592,940 autosomal assays, the rWGS genotypes with ≥ 20 × coverage and  a Phred score ≥ 30 were used for calculation of performance metrics. These analyses, i.e., GSA vs. 1KG and rWGS vs. 1KG, showed consistent average metrics and small standard deviations among datasets (Table 9). For the 256 Coriell samples with 1KG data, we observed that rWGS performed better than GSA across all 4 performance metrics ( Fig. 6A; Table 9). Overall average concordance, sensitivity, and specificity for rWGS vs. 1KG were 0.9981, 0.9981 and 0.9991, respectively, whereas for GSA vs. 1KG, they were 0.9932, 0.9927, and 0.9957, respectively. PPV was 0.9977 for rWGS vs. 1KG and was 0.9892 for GSA vs. 1KG (Table 9).

Over 82% of all GSA assays have a PPV = 1
We compared the GSA and rWGS genotypes to the 1KG and computed the PPV. As shown in Fig. 6B, over 82% (476,828) of assays had a PPV of 1 for both the GSA and rWGS. Approximately 1.5% (8710) of rWGS assays had a PPV of 1 when GSA was 0, whereas only 0.12% (699) of GSA assays did when rWGS was 0.

GSA MAP59 secondary findings validated using rWGS, pWGS, and 1KG
Given that > 80% of GSA assays have a PPV = 1, we assessed rare variation detection within the 59 medically actionable predisposition genes (MAP59) defined by the American College of Medical Genetics (ACMG) [35]. Given the expected secondary finding rate of 1-2% [38][39][40] and the limited genomic space profiled by the GSA, we hypothesized 2-3 or fewer samples with GSA-detectable variants in the 261 cohort. Additionally, we hypothesized that comparison of these data to the 1KG and the rWGS data identifies false negative and false positive variants as well as pathogenic variation undetected by the GSA. Focusing on nucleotides with ≥ 20 × rWGS coverage (Fig. 7), we found that an average of 6347 (± 88) sites were genotyped by both rWGS and GSA in any given DNA sample. The GSA vs. rWGS average concordance, sensitivity, specificity, and PPV were 0.99897, 0.99367, 0.99962, and 0.9946, respectively.
For clinically reportable rare variants curated into the managed variant list (MVL), the GSA and rWGS were concordant for a heterozygous variant (MUTYH p.(Gly368Asp); rs36053993) in three samples and across GSA replicates. Two of the 3 samples had 1KG data Fig. 6 Scatter-plot comparison of performance metrics of whole genome sequencing (rWGS) and GSA using 1KG as the benchmark dataset.
A Scatter plots show sample-level performance metrics of rWGS and GSA relative to 1KG reference data. Plots are concordance (top left; blue), sensitivity (top right; orange), specificity (bottom left; green) and positive predictive value (PPV) (bottom right; maroon) respectively. Each dot represents a single sample's performance metric value. B Density scatterplot of each GSA assay's positive predictive value computed for GSA (y-axis) vs. rWGS (x-axis) using 1KG as the benchmark dataset. Each square represents PPV measured for GSA and rWGS relative to 1KG benchmark dataset, and the color indicates the number of assays within each square. The color gradient of each square ranges from 1 assay (dark purple) to 476,828 assays (yellow); therefore, the color on the scatterplot indicates the density of data-points in 2 dimensions and were concordant; one of these two had pWGS data that was also concordant. Highlighting the potential for false positives, rWGS and 1KG data refuted a GSA call of PKP2 p.(Arg355Ter) (rs754912778) in one sample. Conversely, highlighting the potential for false negatives, rWGS and 1KG detected two variants that were not detected by GSA: RB1 p.(Arg661Trp) (rs137853294), which the GSA called homozygous reference in triplicate, and MUTYH p.(Pro391Leu) (rs529008617), which the GSA called "no-call" in triplicate. In summary, the GSA identified 1 pathogenic variant (true positive), 1 false positive, and 2 false negatives (2 assayed and missed) among the MAP59.
The rWGS rate of detection of rare pathogenic variants in the MAP59 genes was 0.034 (3.4%); 7 variants in Fig. 7 Plot of the average percentage of bases within each MAP59 gene covered by whole genome sequencing (rWGS) to a read depth of A ×10 or more (gte10x) B ×20 or more (gte20x) among the 263 samples. Each rWGS nucleotide was required to have a Phred-based quality score of greater than 30 to be considered for this analysis 9 samples from a population of 261. Removing the 3 variants that were not independently confirmed by the 1KG due to lack of 1KG data gives 4 pathogenic variants in 5 individuals from a population of 261 or a rate of 0.019 (1.9%). This range (0.019-0.034) of pathogenic variants in the MAP59 genes is consistent with the published discovery rate [38,39,41,42].

Discussion
We report an approach to analytical validation of the GSA through quality analyses and through assessment of performance by comparison to benchmark datasets and independent whole-genome sequencing data. To the best of our knowledge, this is the first comprehensive analytical validation of the GSA for clinical genotyping. Our findings support and extend recently reported research studies assessing the utility of the GSA for genetic screening in primary immunodeficiency [43], for population-based genomic screening for rare and medically relevant variation [44], and for detecting rare and clinically relevant markers in multiethnic Indian populations [45].
In our study we used call rate and sample contamination as preliminary parameters of quality control for genotype analysis. Call rate is a primary quality control parameter in all genotyping studies [12,13]. A high threshold for call rate not only ensures inclusion of samples with high quality genotype data but also allows, independent of sample DNA quality, for detection of assays that perform poorly. Additionally, sample contamination detection [14] is key in preventing return of false positive genotypes and is demonstrated by our results. While more advanced quality control methods such as Hardy-Weinberg Equilibrium (HWE) test [15], likelihood of error [19], departure from Mendelian inheritance, and pedigree information are used in various research studies [4,20], they are implemented in analyses that follow genotype generation and are dependent on what analyses are subsequently performed using the genotype data. HWE is used to detect genotypes that deviate from the expectation of HWE, and it is typically applied to variants with a MAF of greater than 0.05 [12]. Consequently, because of our interest in variants of lower MAF, we did not implement this QC metric; however, HWE might be useful within certain cut-offs for MAF as implemented by Suratannon et al. [43] and Narang et al. [45]. Similarly, Mendelian inheritance and pedigree information quality control are critical for linkage and segregation analyses and did not apply to our individual-focused assay.
This evaluation of GSA data is consistent with previous studies that demonstrated the utility of sample data quality metrics like genotype call-rate, p10GC, and DNA contamination detection [11,22]. By analysis of replicates, we show that the majority of the GSA data are highly reproducible. Outliers arose either from positioning along the edges of the Illumina BeadChip or from contamination. Characterization of each GSA assay by variation class, type, genomic DNA complexity, and alternate allele frequency showed that the GSA has the highest performance for SNVs and transition nucleotide changes in genomic regions of high complexity. In contrast, assays interrogating low-complexity regions, rare alleles, or transversions performed poorly. Transversions between complementary nucleotides likely performed poorly because of the characteristics of the assays for these particular transversions (Additional file 1: Sect. S7). Also, consistent with previous reports [46][47][48], assays for rare alleles (< 0.001) had lower performance and might be improved by using algorithms for rare variant detection [10,31,32] or joint-calling [22] rather than the default genotype caller (GenCall). These should be considered in the future to improve detection of rare variants by genotyping chips.
The analytical framework implemented in this study followed a three-way analysis (GSA-rWGS-1KG) to assess the strengths and limitations of individual GSA assays. Unlike many published analyses in which WGS is the test dataset and the BeadArray genotypes are the truth [25][26][27]30], our study had the BeadArray as the test dataset and WGS as the truth. The reversal of test and truth datasets is a major challenge for comparing our results to the published literature. To overcome this challenge, we ensured that the rWGS data had performance metrics (concordance = 0.9981) comparable to that previously published (concordance = 0.9984 [25]). The three-way analysis framework also allowed detection of false positive and false negative genotypes on the GSA platform. Though not evaluated in the current study, the three-way comparison framework in our analysis allows for modeling of genotyping-error specific to variation classes and categories triaged during characterization of the GSA.
Over 82% of assays on the GSA returned genotypes with a high positive predictive value (PPV). The GSA detected some pathogenic variation (MAP59) in the test dataset of 261 Coriell samples, and these variants were independently validated by either the 1KG data or the rWGS/pWGS data or both. Although we attempted to compare GSA results to other chip results, the comparison to previous work was impeded by differences in probe content and density as well as chip design (e.g., 610 k assays on GSA, vs. 247 k assays on HumanExome chip). Some of the pros and cons of using the GSA are summarized in Table 10 below.
The test characteristics of the GSA compared to WGS clearly show that the GSA is not a diagnostic genomic test for individuals with rare disorders. As shown by our MAP59 results and recent research studies [43,44], the GSA lacks robustness for genotyping rare variants as well as probes for detection of private familial disease variants. On the other hand, we show that the GSA has the analytical robustness to serve as a clinical screen for genotypes for which one can establish robust cluster files for the AA, AB, and BB genotypes. This is most easily accomplished for more common genotypes that contribute to polygenic predispositions to disease, particularly common diseases. Screening of an asymptomatic population to assess the likelihood of predisposition to a disease is well established within medicine, and examples include newborn screening for inborn errors of metabolism, mammography for breast cancer, and cholesterol levels for coronary artery disease [50,51]. A major objective of screening tests is to reduce morbidity and mortality in the subject population through risk stratification to target surveillance, early detection, and treatment. With the characterization of genomic risk for drug responsiveness and predisposition to various cancers and cardiovascular disease [52][53][54], we propose that the GSA offers a potential clinical tool for genomic screening.

Limitations of our study
Our comparison of BeadChip arrays to NGS and benchmark datasets has some limitations. Firstly, we evaluated our dataset using accepted algorithms. This did not take into account the benefits of consensus genotyping by multiple algorithms for GSA or NGS data; Hwang et al. found that consensus genotyping minimized false findings [47]. Secondly, cell-line derived variation or lowlevel somatic variation might also have contributed to differences between datasets [25]. Thirdly, we did not analyze variants close to or overlapping other variation in the same location, e.g., insertions/deletions and copy number variation, because these loci are eukaryotic mutation hotspots [55]. Fourthly, our analysis would benefit from comparison to variant benchmark datasets defined in more recent publications [47] and to NIST/ GiAB datasets.

Conclusions
We established the analytical validity of the GSA via a systematic approach utilizing benchmark and rWGS data to evaluate the performance of each assay. We highlight that the GSA assays interrogating rare variants, transversions, and variants within low-complexity regions need careful evaluation. GSA assays can be analytically validated to clinically screen for common genotypes predisposing to disease.


s2fieldsofstudy: Biology Chemistry

text:
Ultra-Large-Scale Screening of Natural Compounds and Free Energy Calculations Revealed Potential Inhibitors for the Receptor-Binding Domain (RBD) of SARS-CoV-2

The emergence of immune-evading variants of SARS-CoV-2 further aggravated the ongoing pandemic. Despite the deployments of various vaccines, the acquired mutations are capable of escaping both natural and vaccine-induced immune responses. Therefore, further investigation is needed to design a decisive pharmacological treatment that could efficiently block the entry of this virus into cells. Hence, the current study used structure-based methods to target the RBD of the recombinant variant (Deltacron) of SARS-CoV-2, which was used as a model variant. From the virtual drug screenings of various databases, a total of four hits were identified as potential lead molecules. Key residues were blocked by these molecules with favorable structural dynamic features. The binding free energies further validated the potentials of these molecules. The TBE for MNP was calculated to be −32.86 ± 0.10 kcal/mol, for SANC00222 the TBE was −23.41 ± 0.15 kcal/mol, for Liriodenine the TBE was −34.29 ± 0.07 kcal/mol, while for Carviolin the TBE was calculated to be −27.67 ± 0.12 kcal/mol. Moreover, each complex demonstrated distinct internal motion and a free energy profile, indicating a different strategy for the interaction with and inhibition of the RBD. In conclusion, the current study demands further in vivo and in vitro validation for the possible usage of these compounds as potential drugs against SARS-CoV-2 and its variants.

Introduction
The suffering of humanity due to coronavirus disease 2019 (COVID- 19), which is characterized by severe respiratory complexities, is linked to an offensive human pathogen discovery, we used the Deltacron variant as a model virus, which is considered to be the current global health issue. The ability of the Deltacron variant to evade the host immune or vaccine-induced immune response is enabled by its acquired mutations. To date, proteomics-and genomics-based vaccine technologies are employed to design potential vaccine candidates. However, the assimilated unique alterations in the amino acid sequences of the SARS-CoV-2 spike protein provide it with a different ability to evade the immune responses instigated by these vaccines. In this regard, the development of novel drugs or small-molecule inhibitors is of great interest to avoid such problems, i.e., immune evasion. To accelerate the search for potential lead candidates, advanced computational algorithms can fast-track this process. These computational methods have been previously used to design potential inhibitors for SARS-CoV-2. For instance, molecular modelling and simulation-based approaches have been used to target the RBD of the Omicron variant [13]. Similarly, drug repurposing and molecular screening-based approaches also identified potential small molecules for other proteins such as 3CLpro, PLpro, RdRp, S protein, and NRP1 [16,17]. Thus, considering the wide range of applicability of these methods, we also used a large-scale drug screening approach for the RBD of the Deltacron variant. We screened several databases such as the South African natural compounds database (SANCDB), East African natural compounds database, North-East African natural compounds database, and Marine natural compounds database to identify potential smallmolecule inhibitors for SARS-CoV-2. The identified potential molecules were validated by using molecular simulation approaches with which the stability, compactness, flexibility, dynamic motions, and binding free energy calculations for each complex were performed. For this purpose, mutations in the RBD of the Deltacron variant were retrieved from the literature. The structure of the RBD was modeled using AlphaFold 2.0. The methodological flow of the study is presented in Figure 1.
The emergences of ancestral SARS-CoV-2 and now its new variants have generated a distressing situation that has challenged global health, the effectiveness of the marketed therapeutics, and the severity of the disease. For instance, many variants are reported to increase infection rates and frequency. In order to provide a framework for new drug discovery, we used the Deltacron variant as a model virus, which is considered to be the current global health issue. The ability of the Deltacron variant to evade the host immune or vaccine-induced immune response is enabled by its acquired mutations. To date, proteomics-and genomics-based vaccine technologies are employed to design potential vaccine candidates. However, the assimilated unique alterations in the amino acid sequences of the SARS-CoV-2 spike protein provide it with a different ability to evade the immune responses instigated by these vaccines. In this regard, the development of novel drugs or small-molecule inhibitors is of great interest to avoid such problems, i.e., immune evasion. To accelerate the search for potential lead candidates, advanced computational algorithms can fast-track this process. These computational methods have been previously used to design potential inhibitors for SARS-CoV-2. For instance, molecular modelling and simulation-based approaches have been used to target the RBD of the Omicron variant [13]. Similarly, drug repurposing and molecular screening-based approaches also identified potential small molecules for other proteins such as 3CLpro, PLpro, RdRp, S protein, and NRP1 [16,17]. Thus, considering the wide range of applicability of these methods, we also used a large-scale drug screening approach for the RBD of the Deltacron variant. We screened several databases such as the South African natural compounds database (SANCDB), East African natural compounds database, North-East African natural compounds database, and Marine natural compounds database to identify potential smallmolecule inhibitors for SARS-CoV-2. The identified potential molecules were validated by using molecular simulation approaches with which the stability, compactness, flexibility, dynamic motions, and binding free energy calculations for each complex were performed. For this purpose, mutations in the RBD of the Deltacron variant were retrieved from the literature. The structure of the RBD was modeled using AlphaFold 2.0. The methodological flow of the study is presented in Figure 1. Methodological flow of the whole study. The first step is the modelling of the structure using AlphaFold. In the second step, database retrieval, preparation, and screening are involved. In the third step, the re-evaluation of the top hits from each database is performed, while in the fourth and fifth steps, the validation of the top hits using molecular simulation is performed. Methodological flow of the whole study. The first step is the modelling of the structure using AlphaFold. In the second step, database retrieval, preparation, and screening are involved. In the third step, the re-evaluation of the top hits from each database is performed, while in the fourth and fifth steps, the validation of the top hits using molecular simulation is performed.
The structure of the whole S protein with its domain's organization is shown in Figure 2A, while the RBD structure with its respective mutations is shown in Figure 2B. We also superimposed the structure of the wild-type RBD with that of the Deltacron RBD, which demonstrated an RMSd difference of 0.82 Å, thus informing us of the significant The structure of the whole S protein with its domain's organization is shown in Figure 2A, while the RBD structure with its respective mutations is shown in Figure 2B. We also superimposed the structure of the wild-type RBD with that of the Deltacron RBD, which demonstrated an RMSd difference of 0.82 Å, thus informing us of the significant structural variations caused by the acquired mutations. The superimposed structures of the wild-type and Deltacron RBDs are shown in Figure 2C.

Molecular Screening of Natural Compounds Databases
We screened ~55,000 compounds from various databases to identify potential lead molecules for the inhibition of the Deltacron variant of SARS-CoV-2. Among these, 27 compounds reported good docking scores from all the databases, collectively. Among the total twenty-seven compounds, five compounds from the Marine natural compounds database, six compounds from the South African natural compounds database, nine compounds from the North African natural compounds database, one compound from the North-East African natural compounds database, and six compounds from the East African natural compounds database were reported to have competitive docking scores. These top compounds were subjected to further validations by performing another round of screening using AutoDock Vina. The rationale behind the re-evaluation of these compounds was to remove false predictions with smina. The Autodock Vina docking scores ranged from −6.60 to −5.40 kcal/mol. The molecules included 2- (2,3-dihydroxyphenyl [1,3,5]-trien]-4-one (SANC0022) from the SANCDB with docking score −6.32 kcal/mol, and Liriodenine from the EANCD and Carviolin from the NANCDB with reported docking scores of −5.90 kcal/mol, respectively. These four compounds were observed to be the best among the re-evaluated 27 compounds and were subjected to further investigations such as interactions analysis, dynamic stability, protein's packing assessment, flexibility indexing, hydrogen bonding, and binding free energy calculations.
The investigation of the binding mode of each of these compounds demonstrated significant results. As shown in Figure 3A, the MNC compound established several hydrogen bonds with the essential mutated residues in the RBD. The residue Arg403 was reported to be involved in the salt-bridge, while Tyr453, Ty495, Ty501, and His505 were

Molecular Screening of Natural Compounds Databases
We screened~55,000 compounds from various databases to identify potential lead molecules for the inhibition of the Deltacron variant of SARS-CoV-2. Among these, 27 compounds reported good docking scores from all the databases, collectively. Among the total twenty-seven compounds, five compounds from the Marine natural compounds database, six compounds from the South African natural compounds database, nine compounds from the North African natural compounds database, one compound from the North-East African natural compounds database, and six compounds from the East African natural compounds database were reported to have competitive docking scores. These top compounds were subjected to further validations by performing another round of screening using AutoDock Vina. The rationale behind the re-evaluation of these compounds was to remove false predictions with smina. The Autodock Vina docking scores ranged from −6.60 to −5.40 kcal/mol. The molecules included 2-(2,3-dihydroxyphenyl)-N-[2-(1H-imidazol-5-yl)ethyl]-1,3-thiazole-4-carboxamide (MNC) from the CMNPD with docking score −6.60 kcal/mol, (R)-3 ,5,7-Trihydroxy-4 -methoxyspiro [2H-1-benzopyran-3(4H),7bicyclo [4.2.0]-octa [1,3,5]-trien]-4-one (SANC0022) from the SANCDB with docking score −6.32 kcal/mol, and Liriodenine from the EANCD and Carviolin from the NANCDB with reported docking scores of −5.90 kcal/mol, respectively. These four compounds were observed to be the best among the re-evaluated 27 compounds and were subjected to further investigations such as interactions analysis, dynamic stability, protein's packing assessment, flexibility indexing, hydrogen bonding, and binding free energy calculations.
The investigation of the binding mode of each of these compounds demonstrated significant results. As shown in Figure 3A, the MNC compound established several hydrogen bonds with the essential mutated residues in the RBD. The residue Arg403 was reported to be involved in the salt-bridge, while Tyr453, Ty495, Ty501, and His505 were reported to have established hydrophobic contacts with the RBD of the Deltacron variant. The MNC compound also established seven hydrogen bonds with the residues Arg403, Gln409, Asn417, Tyr453, and Ser496. Interestingly, most of the targeted residues were the mutated residues of the Deltacron variant, which help to bind with the host receptor more robustly. Moreover, previous studies targeting the RBD of the Omicron variant also reported the same site for inhibition. Using structural modeling approaches, Abbas et. al reported that Tyr453, Arg493, Ser494, Ser496, Tyr501, and His505 are essential residues targeted by the shortlisted compounds and consequently abrogate the RBD-ACE2 complex in the case of the Omicron variant (B.1.1.529) [13]. Moreover, the role of Arg403, Asn417, Tyr501, and His505 in the higher infectivity of other variants has also been revealed, thus showing promising results for the inhibition of the PPI complex [18]. On the other hand, the SANC00222 compound, with two hydrophobic interactions with Tyr501 and His505 and seven hydrogen bonds, also demonstrated a similar interaction pattern to MNC. The hydrogen bonds involved Asn417, Tyr453, Tyr501, and His505 residues. Herein, Tyr453, which is a novel mutation in the Deltacron variant, was observed to be involved in the interaction with the SANC00222 compound. This shows that these two compounds target the key mutated residues by establishing an essential bonding network and preventing the host receptor from interacting with the RBD, which consequently results in no viral entry and infection. The results for SANC00222 are shown in Figure 3B.
reported to have established hydrophobic contacts with the RBD of the Deltacron variant. The MNC compound also established seven hydrogen bonds with the residues Arg403, Gln409, Asn417, Tyr453, and Ser496. Interestingly, most of the targeted residues were the mutated residues of the Deltacron variant, which help to bind with the host receptor more robustly. Moreover, previous studies targeting the RBD of the Omicron variant also reported the same site for inhibition. Using structural modeling approaches, Abbas et. al reported that Tyr453, Arg493, Ser494, Ser496, Tyr501, and His505 are essential residues targeted by the shortlisted compounds and consequently abrogate the RBD-ACE2 complex in the case of the Omicron variant (B.1.1.529) [13]. Moreover, the role of Arg403, Asn417, Tyr501, and His505 in the higher infectivity of other variants has also been revealed, thus showing promising results for the inhibition of the PPI complex [18]. On the other hand, the SANC00222 compound, with two hydrophobic interactions with Tyr501 and His505 and seven hydrogen bonds, also demonstrated a similar interaction pattern to MNC. The hydrogen bonds involved Asn417, Tyr453, Tyr501, and His505 residues. Herein, Tyr453, which is a novel mutation in the Deltacron variant, was observed to be involved in the interaction with the SANC00222 compound. This shows that these two compounds target the key mutated residues by establishing an essential bonding network and preventing the host receptor from interacting with the RBD, which consequently results in no viral entry and infection. The results for SANC00222 are shown in Figure 3B. Next, we also evaluated the binding mode of Liriodenine from the EANDB and Carviolin from the NANCDB. The results for these compounds are shown in Figure 4A, B. The Liriodenine compound from the EANDB established a single hydrophobic interaction with Phe497 and four hydrogen bonds involving Tyr453, Ser496, Tyr501, and His505. Interestingly, Tyr453 and Phe497 are among the novel mutations acquired by the Deltacron variant; hence, this shows the potential activity of this compound against the RBD of SARS-CoV-2. The binding mode of Liriodenine from the EANDB is shown in Figure 4A. On the other hand, the Carviolin compound from the NANDB reported a similar docking score to Liriodenine; however, variations in the binding pattern were observed. For the two hydrophobic interactions involving the Tyr495 and Tyr501 residues, five hydrogen bonds were reported, including two with Tyr453 and a single hydrogen bond with Ser496, while His505 also established two hydrogen bonds with this compound. The binding mode of Carviolin with the RBD is shown in Figure 4B. In conclusion, the current findings show promising results by targeting the key essential residues required for direct interaction with the host receptor, thus preventing viral entry into the host cell and consequently Next, we also evaluated the binding mode of Liriodenine from the EANDB and Carviolin from the NANCDB. The results for these compounds are shown in Figure 4A,B. The Liriodenine compound from the EANDB established a single hydrophobic interaction with Phe497 and four hydrogen bonds involving Tyr453, Ser496, Tyr501, and His505. Interestingly, Tyr453 and Phe497 are among the novel mutations acquired by the Deltacron variant; hence, this shows the potential activity of this compound against the RBD of SARS-CoV-2. The binding mode of Liriodenine from the EANDB is shown in Figure 4A. On the other hand, the Carviolin compound from the NANDB reported a similar docking score to Liriodenine; however, variations in the binding pattern were observed. For the two hydrophobic interactions involving the Tyr495 and Tyr501 residues, five hydrogen bonds were reported, including two with Tyr453 and a single hydrogen bond with Ser496, while His505 also established two hydrogen bonds with this compound. The binding mode of Carviolin with the RBD is shown in Figure 4B. In conclusion, the current findings show promising results by targeting the key essential residues required for direct interaction with the host receptor, thus preventing viral entry into the host cell and consequently reducing the viral infection. The detailed interaction patterns of these compounds along with their docking scores are presented in Table 1. reducing the viral infection. The detailed interaction patterns of these compounds along with their docking scores are presented in Table 1.

Structural and Dynamic Stability Assessment of the Top-Scoring Complexes
The computation of RMSd in an apo or holo system demonstrates essential information regarding the function of a protein. It is essential in order to determine the binding stability of a small molecule in a complex with a receptor that is usually a protein/druggable target. For this purpose, estimation of the root-mean-square deviation or, shortly, RMSd is necessary to comprehend the binding association between the two molecules. The RMSd informs us of the stability of a system in a dynamic environment solvated with water molecules and neutralized by other factors in order to essentially understand the effects of the interacting partners. The determination of the dynamic stabilities of small molecules in particular is important in order to gauge their inhibitory potentials and therapeutic efficacies. A small molecule may also destabilize the structural integrity of a receptor to hinder the normal functionality of that particular receptor. Hence, considering the utmost significance of dynamic stability, we herein also used RMSd by executing the simulation trajectories to demonstrate the stability of each complex over the course of the simulation time. The calculated RMSd for each complex is shown in Figure 5A-D. The apo-RBD shown in Figure 5A reveals a similar RMSd pattern as the RMSd pattern of the wild type during 1-180 ns, and then abruptly increased. Each complex demonstrated comparatively more stable dynamic behaviors throughout the simulation compared with the apo system. The snapshots at 190 ns, 230 ns, 250 ns, and 300 ns from the apo system were

Structural and Dynamic Stability Assessment of the Top-Scoring Complexes
The computation of RMSd in an apo or holo system demonstrates essential information regarding the function of a protein. It is essential in order to determine the binding stability of a small molecule in a complex with a receptor that is usually a protein/druggable target. For this purpose, estimation of the root-mean-square deviation or, shortly, RMSd is necessary to comprehend the binding association between the two molecules. The RMSd informs us of the stability of a system in a dynamic environment solvated with water molecules and neutralized by other factors in order to essentially understand the effects of the interacting partners. The determination of the dynamic stabilities of small molecules in particular is important in order to gauge their inhibitory potentials and therapeutic efficacies. A small molecule may also destabilize the structural integrity of a receptor to hinder the normal functionality of that particular receptor. Hence, considering the utmost significance of dynamic stability, we herein also used RMSd by executing the simulation trajectories to demonstrate the stability of each complex over the course of the simulation time. The calculated RMSd for each complex is shown in Figure 5A-D. The apo-RBD shown in Figure 5A reveals a similar RMSd pattern as the RMSd pattern of the wild type during 1-180 ns, and then abruptly increased. Each complex demonstrated comparatively more stable dynamic behaviors throughout the simulation compared with the apo system. The snapshots at 190 ns, 230 ns, 250 ns, and 300 ns from the apo system were retrieved and superimposed onto the native conformation. The comparative structural investigation revealed that three regions, i.e., 333-350, 470-494, and 510-526, demonstrated significant structural perturbation, which is a major deviation in the RMSd. More interestingly, the deviation that occurs between 470 and 494 signifies the important region required for direct interaction with the host receptor. The two important loops reside in this region that is essential for binding. Furthermore, the transition of essential secondary structure elements was also observed. The superimposed structures at these time intervals are shown in Supplementary Figure S1. The RMSd of the MNC compound initially increases gradually and attains equilibrium at 50 ns. During the 300 ns simulation, no significant deviation was recorded; however, a minor structural perturbation can be seen between 75 and 85 ns, and then a flattened graph continues until 300 ns. The minor perturbation could be due to the loops distributed in the RBD, which was also witnessed in other studies targeting the hACE2-RBD interface [13]. The RMSd is seen to show variations until 175 ns. However, after this, till the end of the simulation (300 ns), a flat graph can be seen, which shows the stabilized binding of the MNC compound. Moreover, the SANC00222 complex with the RBD reports a comparatively higher RMSd but is stable. The graph shows an initial increase in the RMSd during the first 50 ns, and then an increase/decrease pattern until 135 ns is observed. Although a fluctuation between 51 and135 ns can be seen, it does not demonstrate instability, as it is very minor and the lower RMSd is regained afterward. The RMSd after 135 ns decreases continuously and follows the same trend until 200 ns. This trend continues until 300 ns with no significant structural perturbation, thus demonstrating a comparatively unstable dynamic compared with MNC. This, consequently, shows that although an increase in the RMSd is observed, overall, the complex reports stable dynamics and thus produces an inhibitory potential to abstain from the hACE2-RBD interaction. Interestingly, the two other complexes, i.e., Liriodenine-and Carviolin-RBD, demonstrated more stable dynamics than the first two complexes. The Liriodenine complex reaches equilibrium at 25 ns, followed by a smaller deviation at 45-50 ns, then continues to follow a uniform pattern until 115 ns. The RMSd then increases a little, although with no structural perturbation, and then decreases back down after reaching 165 ns. A stable and straight graph is then observed until 200 ns. Moreover, the Carviolin complex reaches stability at 1.5 Å and attains equilibrium at 15 ns. Except for a minor deviation between 70 and 75 ns, this complex demonstrates stable behavior until 200 ns. This complex was observed to be the most stable among the four complexes. For all the complexes, the structures reach stability at 1.5-2.5 Å, and the average RMSd was calculated to be 1.5-2.5 Å. The Liriodenine and Carviolin complexes during the last 100 ns (201-300 ns) further continue the stability pattern and report stable dynamics. The RMSds of all the complexes are shown to converge with the wild type, which shows that all the complexes attain a similar conformational configuration, which consequently shows the effect of the ligand binding in the later part of the simulation. This demonstrates that the ligand-bound complexes had already reached stability and the 300 ns simulation was enough to reveal the pharmacological properties of these compounds. Overall, these findings show the stable binding and dynamic behaviors of each molecule in a complex with RBD, thus favoring the pharmacologically essential features for the inhibition of SARS-CoV-2.

Assessing the Structural Packing in a Dynamic Environment
Next, we evaluated the structural compactness of each protein in a dynamic solution as the radius of gyration (Rg). Rg may provide essential information regarding the kinetics, inhibition, and thermodynamics of a protein that may generate new insights into the biological processes. It may, alternatively, provide information regarding the binding and unbinding events that occurred during the simulation. Interrelationships between some structural, thermodynamic, and kinetic parameters of proteins have been reported in recent years [19]. Keeping the high importance of Rg in the biological context herein, we also calculated Rg as a function of time using 15,000 structural frames from the 300 ns simulation. The Rg pattern for the wild type demonstrates a uniform graph that converges with other complexes, and the average Rg was reported to be 18.60 Å. As presented in Figure  6A, the Rg graph follows the same pattern as the RMSd graph where a gradual increase over the time was observed. The initial Rg increases gradually between 0 and35 ns, and then until 100 ns, smaller deviations are observed. Afterward, the Rg slightly increases with time, thus implying that binding and unbinding events occurred during the simulation. For the MNP-RBD complex, the average Rg was calculated to be 18.50 Å. On the other hand, the SANC00222 Rg results also strongly corroborate with the RMSd results. The Rg during the first 50 ns demonstrates a uniform pattern fluctuating between 18.40 and18.60 Å and then increases abruptly, reaching 19.0 Å. The Rg decreases back down and continues to follow the same pattern with minor deviations between 51 and135 ns. Further decline in the Rg value is observed particularly between 140 and160 ns and then, again, reaches the average level of 18.50 Å. The potential rationale behind the increase and decrease in the protein size is that it is due to the free movement of the three loops that are present at the interface of the RBD required for direct interaction with the host receptor ACE2. The Rg graph for the SANC00222-RBD complex is shown in Figure 6B. Interestingly, the Rg graphs for Liriodenine and Carviolin also present the same pattern as those of the RMSds. These two complexes demonstrate stable Rg patterns and lower average Rg values. The Liriodenine-RBD complex starts at 18.40 Å, following the trend till 40 ns, and then minor deviations between 41 and60 ns are recorded. Subsequently, the Rg is recorded as decreasing until 110 ns, and then a gradual increase is observed. Though the Rg values slightly increase, no significant deviation is experienced until the end of the simulation. On the other hand, the Carviolin-RBD complex reports higher Rg values between 40 and80 ns and between 170 and185 ns. The graph is observed to be straight. For these two complexes, the average Rg values were calculated to be 18.60 Å and 18.73 Å, respectively. The

Assessing the Structural Packing in a Dynamic Environment
Next, we evaluated the structural compactness of each protein in a dynamic solution as the radius of gyration (Rg). Rg may provide essential information regarding the kinetics, inhibition, and thermodynamics of a protein that may generate new insights into the biological processes. It may, alternatively, provide information regarding the binding and unbinding events that occurred during the simulation. Interrelationships between some structural, thermodynamic, and kinetic parameters of proteins have been reported in recent years [19]. Keeping the high importance of Rg in the biological context herein, we also calculated Rg as a function of time using 15,000 structural frames from the 300 ns simulation. The Rg pattern for the wild type demonstrates a uniform graph that converges with other complexes, and the average Rg was reported to be 18.60 Å. As presented in Figure 6A, the Rg graph follows the same pattern as the RMSd graph where a gradual increase over the time was observed. The initial Rg increases gradually between 0 and 35 ns, and then until 100 ns, smaller deviations are observed. Afterward, the Rg slightly increases with time, thus implying that binding and unbinding events occurred during the simulation. For the MNP-RBD complex, the average Rg was calculated to be 18.50 Å. On the other hand, the SANC00222 Rg results also strongly corroborate with the RMSd results. The Rg during the first 50 ns demonstrates a uniform pattern fluctuating between 18.40 and 18.60 Å and then increases abruptly, reaching 19.0 Å. The Rg decreases back down and continues to follow the same pattern with minor deviations between 51 and 135 ns. Further decline in the Rg value is observed particularly between 140 and 160 ns and then, again, reaches the average level of 18.50 Å. The potential rationale behind the increase and decrease in the protein size is that it is due to the free movement of the three loops that are present at the interface of the RBD required for direct interaction with the host receptor ACE2. The Rg graph for the SANC00222-RBD complex is shown in Figure 6B. Interestingly, the Rg graphs for Liriodenine and Carviolin also present the same pattern as those of the RMSds. These two complexes demonstrate stable Rg patterns and lower average Rg values. The Liriodenine-RBD complex starts at 18.40 Å, following the trend till 40 ns, and then minor deviations between 41 and 60 ns are recorded. Subsequently, the Rg is recorded as decreasing until 110 ns, and then a gradual increase is observed. Though the Rg values slightly increase, no significant deviation is experienced until the end of the simulation. On the other hand, the Carviolin-RBD complex reports higher Rg values between 40 and 80 ns and between 170 and 185 ns. The graph is observed to be straight. For these two complexes, the average Rg values were calculated to be 18.60 Å and 18.73 Å, respectively. The Rg graphs for Liriodenine and Carviolin are shown in Figure 6C,D. In the extended simulation, no significant deviation was observed till the end of the simulation, which shows that these complexes had already reached stability, and the simulation time was enough to conclude the findings. Overall, these results show that either due to the larger movement of the interface loops or the binding/unbinding of these small molecules, they experienced different interaction paradigms during the simulation. Rg graphs for Liriodenine and Carviolin are shown in Figure 6C, D. In the extended simulation, no significant deviation was observed till the end of the simulation, which shows that these complexes had already reached stability, and the simulation time was enough to conclude the findings. Overall, these results show that either due to the larger movement of the interface loops or the binding/unbinding of these small molecules, they experienced different interaction paradigms during the simulation.

Figure 6. Determination of structural packing in a dynamic environment for the apo-RBD and top complexes. (A) shows the Rg patterns for apo-RBD and MNP-RBD complex (B) shows the Rg patterns for apo-RBD and SANC00222RBD complex, (C) shows the Rg graphs for apo-RBD and
Liriodenine-RBD complex, while (D) represents the Rg graphs for apo-RBD and Carviolin-RBD complex over the course of the 300 ns simulation time.

Indexing the Residual Flexibility
The residual flexibility has an important association with the kinetics and functionality of biological molecules. This approach can be used to determine essential residues that play a significant role in molecular recognition, bio-catalytic mechanisms, manual protein designs, and in the coupling of two or more proteins/interacting partners. In this regard, conformational variations that span a wide variety of amplitude scales are typically linked to protein function. Protein dynamics have been shown to be crucial; thus, knowing about protein flexibility is just as important as knowing about protein structure when it comes to understanding a protein's function and to improving drug development [20]. Herein, we indexed the flexibility of each residue during the course of the simulation. As presented in Figure 7A, all the complexes show a similar pattern of flexibility for most of the regions. A minimal fluctuation for most of the residues is observed. However, regions such as 358-375 and 473-500 report higher fluctuations. The region 510-525 also reports a higher fluctuation. The residues 470-500 comprise the three important loops required for direct interactions with the host receptor ACE2 and demonstrate higher fluctuations. To determine the flexibility of these loops, the residue flexibility was calculated and is shown in Figure 7B. It can be seen that all the regions, particularly 473-490, exhibit higher fluctuations which, consequently, differentiates the dynamic behavior of each complex induced by the binding of each small molecule.

Indexing the Residual Flexibility
The residual flexibility has an important association with the kinetics and functionality of biological molecules. This approach can be used to determine essential residues that play a significant role in molecular recognition, bio-catalytic mechanisms, manual protein designs, and in the coupling of two or more proteins/interacting partners. In this regard, conformational variations that span a wide variety of amplitude scales are typically linked to protein function. Protein dynamics have been shown to be crucial; thus, knowing about protein flexibility is just as important as knowing about protein structure when it comes to understanding a protein's function and to improving drug development [20]. Herein, we indexed the flexibility of each residue during the course of the simulation. As presented in Figure 7A, all the complexes show a similar pattern of flexibility for most of the regions. A minimal fluctuation for most of the residues is observed. However, regions such as 358-375 and 473-500 report higher fluctuations. The region 510-525 also reports a higher fluctuation. The residues 470-500 comprise the three important loops required for direct interactions with the host receptor ACE2 and demonstrate higher fluctuations.
To determine the flexibility of these loops, the residue flexibility was calculated and is shown in Figure 7B

Analysis of the Hydrogen Bonding during Simulation
The inhibitory properties of the interacting compounds can be comprehended better by quantifying the hydrogen bonds between the protein and drug molecules. One of the critical elements that define the binding drug's inhibitory actions is hydrogen bonding. The binding affinity of the interacting molecules changes as the number of hydrogen bonds fluctuates. The estimation of the hydrogen-bonding paradigm has previously been utilized in order to demonstrate how the emerging SARS-CoV-2 variants preferentially impact the binding with the hACE2 receptor by changing the hydrogen-bonding network Moreover, it has also been used in the discovery of novel drugs that target the receptorbinding domain (RBD) of the Omicron variant [13]. Firstly, we calculated the average number of hydrogen bonds in each trajectory to see how the bonding was affected throughout the simulation time. As presented in Figure 8A-D, the hydrogen bonding demonstrates changes over the course of the simulation time in each complex. In each complex, the average numbers of hydrogen bonds were calculated to be 75 in apo-RBD, 80 in MNP-RBD, 78 in SANC00222-RBD, 802 in Liriodenine-RBD, and 79 in Carviolin-RBD. This shows that the additional bonds were probably due to the binding of ligands with the RBD interface.

Gibbs Free Energy Calculation to Re-Evaluate the Best Conformations
The pharmacological effects of any particular drug/small molecule can be pre-clinically validated by predicting the accurate binding conformation and re-evaluating the binding free energy using the MM/GBSA or MM/PBSA methods. These methods are the most popular and reliable methods for re-ranking a potent inhibitor against a particular receptor. They are virtually reasonable and have the advantage over other rational methods. Hence, in view of its accuracy and applicability, we also considered using MM/GSBA to calculate the binding free energy for each complex. In the case of the MNP-RBD complex, the vdW energy was computed to be −41.11 ± 0.11 kcal/mol, the electrostatic energy was computed to be −10.34 ± 0.10 kcal/mol, while the total binding energy was estimated to be −32.86 ± 0.10 kcal/mol. On the other hand, for the SANC00222-RBD complex, the vdW was calculated to be −29.96 ± 0.10 kcal/mol, the electrostatic energy was estimated to be −13.31 ± 0.13 kcal/mol, while the total binding free energy was computed to be −23.41 ± 0.15 kcal/mol. Moreover, the Liriodenine-RBD complex demonstrated a vdW of −31.73 ± 0.08 kcal/mol, an electrostatic energy of −8.12 ± 0.08 kcal/mol, while the total binding free energy was estimated to be −34.29 ± 0.07 kcal/mol. Finally, the total binding free energy was reported to be −27.67 ± 0.12 kcal/mol for the Carviolin-RBD complex, while the vdW and electrostatic energies were calculated to be −28.98 ± 0.13 and −6.54 ± 0.09 kcal/mol, respectively. This shows that the four hits could inhibit the RBD more strongly and require further in vivo and in vitro validation for their possible usage as potential drugs against SARS-CoV-2 and its variants. The binding free energy results are shown in Table 2.

Clustering of Protein Motion in Trajectories
For each complex, a PCA analysis was conducted to group similar trajectories and conformations. This method has been widely applied to filter large concerted motions from the set of structures produced by experiments and molecular computations. As shown in Figure 10A-E, each drug reports varied dynamic motions of 47%, 45%, 48%, 39%, and 44% motion of the total internal motion by the first three eigenvectors of each complex. In the illustration, the two colors orange and blue stand in for the two conformations, while purple stands in for the two transition phases. In conclusion, this demonstrates that the final hits exhibit distinct internal motions, indicating their different strategies for the interaction with and inhibition of the RBD.

Free Energy Landscape Analysis
Next, the conformational alterations induced by the binding of these small molecules were further explored by utilizing the free energy landscape. The constructed FEL maps for all the complexes are presented in Figure 10F-J. For each complex, the conformational states are separated by a subspace. For each complex, the free energy profiles are altered, and the overall results reflect the conformational and dynamic changes that favor the bindings of the identified hits to the respective sites.

Conclusions
In conclusion, the current study used structure-based methods to target the RBD of the recombinant variant (Deltacron) of SARS-CoV-2. From the virtual drug screenings of various databases, four hits were identified as potential lead molecules. Key residues were blocked by these molecules with favorable structural dynamic features. The binding free energies further validated the potentials of these molecules. Further validations such as binding free energy calculations, PCAs, and free energy profiles confirmed the pharmacological potentials of these top hits. The current study demands further in vivo and in vitro validation for the possible usage of these compounds as potential drugs against SARS-CoV-2 and its variants.

Structures Retrieval and Molecular Screening
The accession ID 6M0J was used to obtain the 3D structure of the SARS-CoV-2 spike glycoprotein from the Protein Databank [21]. AlphaFold 2.0 in Google colab was used to model the structure of the RBD of the Deltacron variant. The South African natural compounds database (SANCDB), East African natural compounds database, North-East African natural compounds database, Chinese medicine database, and Marine natural compounds database were searched for the retrieval of small molecules, and then the Open Babel and FAF-drugs4 web tools were used for the preparation and estimation of the Lipinski's rule violation [22,23]. The compounds which did not follow the R5 rules were removed, and the remaining compounds were processed for molecular screening. The Smina docking tool with the AutoDock scoring function was utilized for the initial virtual

Free Energy Landscape Analysis
Next, the conformational alterations induced by the binding of these small molecules were further explored by utilizing the free energy landscape. The constructed FEL maps for all the complexes are presented in Figure 10F-J. For each complex, the conformational states are separated by a subspace. For each complex, the free energy profiles are altered, and the overall results reflect the conformational and dynamic changes that favor the bindings of the identified hits to the respective sites.

Conclusions
In conclusion, the current study used structure-based methods to target the RBD of the recombinant variant (Deltacron) of SARS-CoV-2. From the virtual drug screenings of various databases, four hits were identified as potential lead molecules. Key residues were blocked by these molecules with favorable structural dynamic features. The binding free energies further validated the potentials of these molecules. Further validations such as binding free energy calculations, PCAs, and free energy profiles confirmed the pharmacological potentials of these top hits. The current study demands further in vivo and in vitro validation for the possible usage of these compounds as potential drugs against SARS-CoV-2 and its variants.

Structures Retrieval and Molecular Screening
The accession ID 6M0J was used to obtain the 3D structure of the SARS-CoV-2 spike glycoprotein from the Protein Databank [21]. AlphaFold 2.0 in Google colab was used to model the structure of the RBD of the Deltacron variant. The South African natural compounds database (SANCDB), East African natural compounds database, North-East African natural compounds database, Chinese medicine database, and Marine natural compounds database were searched for the retrieval of small molecules, and then the Open Babel and FAF-drugs4 web tools were used for the preparation and estimation of the Lipinski's rule violation [22,23]. The compounds which did not follow the R5 rules were removed, and the remaining compounds were processed for molecular screening. The Smina docking tool with the AutoDock scoring function was utilized for the initial virtual molecular search. It has certain advantages, such as its supporting multi-ligand files, various ligand molecular formats, and additional energy terms, such as desolvation and electrostatics energies [24,25]. Moreover, the automatic box creation option makes it more feasible for users. The active site was identified on the basis of the previous literature, and the residues were crosschecked for accuracy. The 50 top hits were re-docked for the removal of any false-positive predictions.

Molecular Dynamics Simulation of Protein-Ligand Complexes
The inhibitory features of potential lead molecules and the validation of dynamic behaviors were assessed with the help of molecular simulation. The complexes of compounds and the RBD (Deltacron variant) were subjected to the AMBER21 package for estimation of tier flexibility, stability, and packing. The coordinates, topology, and solvation of the drugs were achieved with an antechamber algorithm and optimal point charge (OPC) water model, respectively. The addition of Na+ counter ions was essential for neutralization purposes. For improved and optimized structures, the minimization of each complex was performed for a cumulative 20,000 steps. After improvement, the complexes were subjected to heating and equilibration for 50 ns. The particle mesh Ewald (PME) algorithm was employed for long-range contacts, whereas van der Waals (vdW) interactions for contacts within a range of 1.4 nm and Columbic interactions for short-range contacts were considered. A total run of 1500 ns and a complex-production run of 300 ns were performed with a time step of 2 fs. CPPTRAJ and PTRAJ were used to obtain the trajectories of each simulation [26]. A root-mean-square deviation (RMSD) analysis for the function of time was performed to check the structural stability by using the following equation: where d denotes the position difference between atoms, and i is used for the original and superimposed structure. On the basis of the B-factor, root-mean-square fluctuation (RMSF) was calculated and it is the most essential constraint for the calculation of flexibility for all protein residues. RMSF values can be obtained with the following equation: Thermal f actor or B-f actor = [(8π × 2)/3] (ms f )

The Binding Free Energy Calculations
The complexes' binding free energies were estimated using the MMPBSA.PY script [27]. It calculates each and every energy term that contributes to the total binding energy. We used the whole simulation trajectories for the binding free energy calculations. For all the complexes, free energy was calculated with the following equation: Free binding energy is represented by "∆G bind " in the equation. The specific energy terms that contributed to the total free energy were obtained with the following equation: G = G bond + G electrostatic + G vdW + G polar + G non-polar (4)

Principal Component Analysis (PCA)
A PCA was performed to capture the motion trajectories with the help of the CPPTRAJ package. With this method, high-amplitude motions were mapped by considering the eigenvalues of the principal motions [28]. From the obtained trajectories, the PC1 and PC2 high-amplitude fluctuations were selected for plotting against each other.

Free Energy Landscape (FEL)
The two PCs (PC1 and PC2) were used for the contour map construction in order to obtain the lowest energy state (metastable state) and also show the lowest conformations obtained by each complex. The subspaces between the states represent the transition states [29]. The Gromacs module g_sham was used for the free energy landscape (FEL). PC1 and PC2 were used for the FEL calculations by using the following equation: ∆G(PC1, PC2) = K B ln P(PC1, PC2) Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/molecules27217317/s1, Figure S1: