A cost-effective sequencing method for genetic studies combining high-depth whole exome and low-depth whole genome

Whole genome sequencing (WGS) at high-depth (30X) allows the accurate discovery of variants in the coding and non-coding DNA regions and helps elucidate the genetic underpinnings of human health and diseases. Yet, due to the prohibitive cost of high-depth WGS, most large-scale genetic association studies use genotyping arrays or high-depth whole exome sequencing (WES). Here we propose a cost-effective method which we call “Whole Exome Genome Sequencing” (WEGS), that combines low-depth WGS and high-depth WES with up to 8 samples pooled and sequenced simultaneously (multiplexed). We experimentally assess the performance of WEGS with four different depth of coverage and sample multiplexing configurations. We show that the optimal WEGS configurations are 1.7–2.0 times cheaper than standard WES (no-plexing), 1.8–2.1 times cheaper than high-depth WGS, reach similar recall and precision rates in detecting coding variants as WES, and capture more population-specific variants in the rest of the genome that are difficult to recover when using genotype imputation methods. We apply WEGS to 862 patients with peripheral artery disease and show that it directly assesses more known disease-associated variants than a typical genotyping array and thousands of non-imputable variants per disease-associated locus.

shows average values and standard errors.Supplementary Figure 22.Precision of variant calls in targeted regions in ten replicates of HG002 control sample in PAD study when using WEGS8P,4X stratified by depth of coverage (DP).The figure represents SNV and InDel calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The HG002 sample was used as a control during the sequencing of PAD samples, which resulted in 10 replicates.We calculated the variant precision rates for each HG002 replicate against the GIAB reference data.Shapes represent the average values, while vertical lines represent standard errors.For InDels, DP at the start position of InDel was used.

Supplementary Table 1. Precision and recall of variant calls in ten replicates of HG002
control sample in PAD study when using WEGS8P,4X.The HG002 sample was used as a control during the sequencing of PAD samples, which resulted in 10 replicates.We calculated the variant recall and precision rates for each HG002 replicate against the GIAB reference data.14.Precision and recall rates of GLIMPSE-imputed and called variants combined when using WEGS.P -precision.R -recall.This table reports the average number for each sample.Each sample was sequenced 4 times using WEGS4P,2X.HG002 and HG004 were sequenced 5 times using WEGS8P,5X.HG003 was sequenced 6 times using WEGS8P,5X.The WES target regions are based on the Agilent V7 capture.Only variants inside the GIAB high-confidence regions were used for estimating precision and recall rates.The percent of missed true variants equals (1 -recall) * 100.The local imputation reference panel combines haplotypes from the 1000 Genomes Project and Human Genome Diversity Project (see Methods).The precision and recall rates for imputed variants when using low depth WGS data only are in Supplementary Tables 20 and 21.

Supplementary Figure 2 .
Comparison of average depths of coverage (DP) in individual target regions in experiments without and with multiplexing.We computed the average depth of coverage (DP) for each target region in Agilent V7 capture in each sequenced individual.DP includes only paired mapped reads and base pairs with minimal Phred-scaled mapping and base qualities of 20.Each point corresponds to a single target region.In total, there were 208,817 non-overlapping autosomal regions.The X-axis shows the median of average DPs in the target region across all individuals sequenced in 4-plex (panel A) and 8plex (panel B) experiments.The Y-axis shows the median of average DPs in the target region across all individuals sequenced without multiplexing.The darker color represents the higher density of the points.The vertical bar on the right of each panel shows the number of points corresponding to each color.The dotted black line corresponds to the 1:1 ratio between DP in the experiments.Supplementary Figure 4. Number of reads in autosomal chromosomes in sequencing experiments with and without multiplexing.The figure shows only paired reads in autosomal chromosomes, excluding reads that are non-primary or supplementary alignments or failed platform/vendor quality checks.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 ⨉ IQR beyond the box.The horizontal line within the box indicates the median value.Open circles are data points corresponding to the sequenced individual exomes.A) Number of paired reads in millions in sequencing experiments without sample multiplexing and when simultaneously sequencing four (4-plex) and eight (8-plex) samples.B) Percent of reads flagged as PCR or optical duplicates.C) Percent of unmapped reads.D) Average Phred-scaled base quality score across all reads in a sequenced sample.Supplementary Figure 5.The number of paired reads in autosomal chromosomes stratified by library preparation batch.The figure shows only paired reads in autosomal chromosomes, excluding reads that are non-primary or supplementary alignments or failed platform/vendor quality checks.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 ⨉ IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the number of paired reads across individual exome in batches 1 and 2, respectively.A) The number of paired reads is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The number of paired reads in the first library preparation batch.C) The number of paired reads in the second library preparation batch.Supplementary Figure 6.Percent of reads flagged as PCR or optical duplicates in autosomal chromosomes stratified by library preparation batch.The figure shows only paired reads in autosomal chromosomes, excluding reads that are non-primary or supplementary alignments or failed platform/vendor quality checks.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 ⨉ IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the percent of duplicated reads across individual exomes in batches 1 and 2, respectively.A) The percent of duplicated reads is stratified by the library preparation batch in experiments without multiplexing, with 4plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The percent of duplicate reads in the first library preparation batch.C) The percent of duplicate reads in the second library preparation batch.Supplementary Figure 7. Percent of unmapped reads in autosomal chromosomes stratified by library preparation batch.The figure shows only paired reads in autosomal chromosomes, excluding reads that are non-primary or supplementary alignments or failed platform/vendor quality checks.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 ⨉ IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the percent of unmapped reads across individual exome in batches 1 and 2, respectively.A) The percent of unmapped reads is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The pvalues above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The percent of unmapped reads in the first library preparation batch.C) The percent of unmapped reads in the second library preparation batch.Supplementary Figure 8.The average quality of reads in autosomal chromosomes stratified by library preparation batch.The figure shows only paired reads in autosomal chromosomes, excluding reads that are non-primary or supplementary alignments or failed platform/vendor quality checks.The average read quality was computed as the average of Phred-scaled base qualities.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 ⨉ IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the average read quality across individual exome in batches 1 and 2, respectively.A) The average read quality is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The average read quality in the first library preparation batch.C) The average read quality in the second library preparation batch.Supplementary Figure 9. Average depths of coverage across all targeted regions in autosomal chromosomes processed with and without UMI-aware deduplication.The average depth of coverage (DP) was computed across target regions in Agilent V7 capture using paired mapped reads and counting only base-pairs with minimal Phred-scaled mapping and base qualities of 20.The box bounds the IQR and Tukey-style whiskers extend to a maximum of 1.5 × IQR beyond the box.The horizontal line within the box indicates median value.Open circles, up-pointing and down-pointing triangles are data points corresponding to the average DP across individual exome processed without, with LocatIt and GATK's UmiAwareMarkDuplicatesWithMateCigar UMI-aware deduplication, respectively.Supplementary Figure 10.Recall and precision of the SNVs and InDels called in sequencing experiments without and with multiplexing.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.Open circles are data points corresponding to the sequenced individual exomes.A) Recall rates of the called SNVs.B) Precision of the called SNVs.C) Recall rates of the called InDels.D) Precision of the called InDels.Supplementary Figure 11.The recall of the SNVs and InDels called in sequencing experiments without and with multiplexing stratified by library preparation batch.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the recall across individual exome in batches 1 and 2, respectively.A) The recall of SNVs is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The recall of SNVs in the first library preparation batch.C) The recall of SNVs in the second library preparation batch.D) The recall of InDels is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.E) The recall of InDels in the first library preparation batch.F) The recall of InDels in the second library preparation batch.Supplementary Figure 12.The precision of the SNVs and InDels called in sequencing experiments without and with multiplexing stratified by library preparation batch.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to the 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the precision across individual exome in batches 1 and 2, respectively.A) The precision of SNVs is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.B) The precision of SNVs in the first library preparation batch.C) The precision of SNVs in the second library preparation batch.D) The precision of InDels is stratified by the library preparation batch in experiments without multiplexing, with 4-plexing and 8-plexing experiments.The p-values above each experiment pair correspond to the one-tailed Wilcoxon rank-sum test.E) The precision of InDels in the first library preparation batch.F) The precision of InDels in the second library preparation batch.Supplementary Figure 13.The number of SNV and InDel calls in sequencing experiments without and with multiplexing.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The solid black line corresponds to the linear regression line, and the dashed black lines correspond to 95% confidence interval.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.A) Number of true positive (TP) SNV calls in sequencing experiments without sample multiplexing and when simultaneously sequencing four (4-plex) and eight (8-plex) samples.B) Number of false positive (FP) SNV calls.C) Number of false negative (FN) SNV calls.D) Number of TP InDel calls.E) Number FP InDel calls.F) Number of FN InDel calls.Supplementary Figure 14.Recall and precision of the single nucleotide variations (SNVs) in sequencing experiments without and with UMI-aware read deduplication.The figure represents SNV calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.Open circles, up-pointing and down-pointing triangles are data points corresponding to the recall and precision in individual exomes processed without, with LocatIt and GATK's UmiAwareMarkDuplicatesWithMateCigar UMI-aware deduplication, respectively.The solid black lines connect pairs of individual exomes with the same underlying sequencing data (i.e.same sequenced sample) but different deduplication approaches.The p-values above experiments with varying levels of multiplexing correspond to the one-tailed Wilcoxon signed-rank test between UMI agnostic and UMI-aware deduplication.A) Recall rates of the called SNVs without UMI-aware compared to UMIaware deduplication using LocatIt.B) Precision of the called SNVs without UMI-aware compared to UMI-aware deduplication using LocatIt.C) Recall rates of the called SNVs without UMI-aware compared to UMI-aware deduplication using GATK's UmiAwareMarkDuplicatesWithMateCigar.D) Precision of the called SNVs without UMI-recall and precision rates in WES experiments with multiplexing before and after adding 2X WGS data.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.Open circles and up-pointing triangles are data points corresponding to the recall and precision in individual multiplexed WES before and after adding 2X WGS data, respectively.The solid black lines connect pairs of individual exomes with the same underlying WES data (i.e.same sequenced sample).The p-values above experiments with varying levels of multiplexing correspond to the one-tailed Wilcoxon signed-rank test between UMI agnostic and UMIaware deduplication.A) Recall rates of the called SNVs with and without 2X WGS.B) Precision rates of the called SNVs with and without 2X WGS.C) Recall rates of the called InDels with and without 2X WGS.D) Precision rates of the called InDels with and without 2X WGS.Supplementary Figure 16.Variant recall and precision rates in WES experiments with multiplexing before and after adding 5X WGS data.The figure represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.Open circles and up-pointing triangles are data points corresponding to the recall and precision in individual multiplexed WES before and after adding 5X WGS data, respectively.The solid black lines connect pairs of individual exomes with the same underlying WES data (i.e.same sequenced sample).The p-values above experiments with varying levels of multiplexing correspond to the one-tailed Wilcoxon signed-rank test between UMI agnostic and UMIaware deduplication.A) Recall rates of the called SNVs with and without 5X WGS.B) Precision rates of the called SNVs with and without 5X WGS.C) Recall rates of the called InDels with and without 5X WGS.D) Precision rates of the called InDels with and without 5X WGS.Supplementary Figure 17.SNVs calling precision and recall rates in no-plexing WES compared to WEGS stratified by library preparation batch.The figure represents SNV calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the individual WES and WEGS in batches 1 and 2, respectively.The p-values above each pair of batches or sequencing methods correspond to the one-tailed Wilcoxon rank-sum test.A) Precision rates of the called SNVs in batches 1 and 2. B) Precision rates of the called SNVs in batch 1. C) Precision rates of the called SNVs in batch 2. D) Recall rates of the called SNVs in batches 1 and 2. E) Recall rates of the called SNVs in batch 1. F) Recall rates of the called SNVs in batch 2. Supplementary Table 7 shows average values and standard errors.Supplementary Figure 18.InDel calling precision and recall rates in no-plexing WES compared to WEGS stratified by library preparation batch.The figure represents InDel calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The box bounds the IQR, and Tukey-style whiskers extend to 1.5 × IQR beyond the box.The horizontal line within the box indicates the median value.Open rectangles and diamonds are data points corresponding to the individual WES and WEGS in batches 1 and 2, respectively.The p-values above each pair of batches or sequencing methods correspond to the one-tailed Wilcoxon rank-sum test.A) Precision rates of the called InDels in batches 1 and 2. B) Precision rates of the called InDels in batch 1. C) Precision rates of the called InDels in batch 2. D) Recall rates of the called InDels in batches 1 and 2. E) Recall rates of the called InDels in batch 1. F) Recall rates of the called InDels in batch 2. Supplementary Table

Table 2 . The average changes in read properties after UMI-aware read deduplication steps relative to the UMI agnostic approach.
Supplementary Table3.Variant calling in whole exome sequencing experiments with and without multiplexing.The table represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.TP -true positives, FP -false positives, FN -false negatives.

Table 4 . The average number of SNVs missed in multiplexing experiments but correctly identified across all no-plexing experiments.
For each multiplexing experiment, we computed the number of false negative (FN) SNV calls that were true positive (TP) in all three no-plexing experiments for the corresponding individual.

Table 5 . The average changes in SNV calling in whole exome sequencing experiments with UMI-aware read deduplication relative to the UMI agnostic approach.
The table represents SNV calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The star symbols represent statistically significant differences when using a one-tailed Wilcoxon signed-rank test: * -P-value < 0.05, ** -P-value < 0.01, *** -Pvalue < 0.001.

Table 6 . The average changes in SNVs and InDels calling in whole exome sequencing experiments when adding additional whole genome sequencing reads relative to pure whole exome sequencing experiments.
The table represents SNV and InDel calls inside the target regions in Agilent V7 capture and the GIAB highconfidence regions.The star symbols represent statistically significant differences when using a one-tailed Wilcoxon signed-rank test: * -P-value < 0.05, ** -P-value < 0.01, *** -P-value < 0.001.

Table 7 . Average variant recall and precision rates in no-plexing WES and WEGS.
The table represents variant calls inside the target regions in Agilent V7 capture and the GIAB high-confidence regions.The star symbols represent statistically significant differences between WES and WEGS when using a one-tailed Wilcoxon rank-sum test: * -Pvalue < 0.05, ** -P-value < 0.01, *** -P-value < 0.001.WEGS values in bold font are higher than the corresponding values in WES.
Supplementary Table8.Average variant recall and precision rates in no-plexing WES and WEGS stratified by library preparation batch.The star symbols represent statistically significant differences between WES and WEGS in the same batch when using a one-tailed Wilcoxon rank-sum test: * -P-value < 0.05, ** -P-value < 0.01, *** -P-value < 0.001.WEGS values in bold font are higher than the corresponding values in WES in the same batch.

Supplementary Table 10. Average genome-wide variant recall and precision rates in 30X WGS, 2X WGS, 5X WGS, and WEGS.
The table represents variant calls in genetic regions overlapping with the GIAB high-confidence regions genome-wide.

Supplementary Table 13. Imputed variants, their allele frequencies, and overlap with true positive (TP) variants in WEGS outside WES target regions.
The arrows ⇧ and ⇩ denote the increase and decrease in AF fold-change (AF ASJ / AF TOPMed) compared to variants where the number of imputed alleles matched the number of true alleles.

Supplementary Table 16. Precision and recall rates of variants imputed using the Minimac4 method genome-wide.
The percent of missed true variants equals (1 -recall) * 100.The local imputation reference panel combines haplotypes from the 1000 Genomes Project and Human Genome Diversity Project (see Methods).

Supplementary Table 17. Overview of genome-wide significant loci associated with peripheral artery disease (PAD) in the 862 WEGS sequenced patients.
Abbreviations: chrchromosome; alt-alternative; ref-reference; freq-frequency; EUR-Europeans; AA-African Americans; EAS-East Asians; DP-average depth; GSA-global screening array (24v3).This table reports the allele frequency and average depth of known genome-wide significant peripheral artery disease loci in the WEGS.The ancestry specific frequency was derived from gnomAD (v3.1.2).

Table 18 . Genome-wide significant loci associated with peripheral artery disease (PAD) and the number of variants within the loci present in WEGS and TOPMed.
For each locus, we counted the number of variants surrounding the lead variant (rsid) within ±500 kilobase (kb) distance.

Table 19. Precision and recall rates of variants imputed from 2X and 5X WGS using the GLIMPSE method
. P -precision.R -recall.Each sample was sequenced one time using 5X WGS with reads split into two lanes.The WES target regions are based on the Agilent V7 capture.Only variants inside the GIAB high-confidence regions were used for estimating precision and recall rates.The percent of missed true variants equals (1 -recall) * 100.The local imputation reference panel combines haplotypes from the 1000 Genomes Project and Human Genome Diversity Project (see Methods).

Supplementary Table 20. Precision and recall rates of GLIMPSE-imputed and called variants combined when using 2X and 5X WGS
. P -precision.R -recall.Each sample was sequenced one time using 5X WGS with reads split into two lanes.The WES target regions are based on the Agilent V7 capture.Only variants inside the GIAB high-confidence regions were used for estimating precision and recall rates.The percent of missed true variants equals (1 -recall) * 100.The local imputation reference panel combines haplotypes from the 1000 Genomes Project and Human Genome Diversity Project (see Methods).