NCBench: providing an open, reproducible, transparent, adaptable, and continuous benchmark approach for DNA-sequencing-based variant calling

Datasets

We have sequenced the NA12878 sample from the genome in a bottle (GIAB) [8] project with two exome sequencing kits at varying average coverages. The genomic DNA from NA12878 was obtained from the NIGMS Human Genetic Cell Repository at the Coriell Institute for Medical Research. The Agilent Human All Exon V7 kit was used to yield a dataset with 182 million paired-end reads sequenced on an Illumina Nova Seq 6000 (211 bp mean insert size and 2 $\times$ 101 bp read length). We used random subsampling to derive two datasets from this that were used in the benchmarking, one with 37.5 million and one with 100 million paired-end reads. The Twist Human Comprehensive Exome (Twist Bioscience, San Francisco, CA, USA) sequencing kit was used according to the manufacturer’s protocol to generate 200 million paired-end reads on an Illumina NovaSeq 6000 (291 bp mean insert size and 2 $\times$ 101 bp read length). The raw reads of the two subsampled Agilent and Twist exome datasets are available via Zenodo.^11,12

Evaluation pipeline

To analyze the quality of the callsets yielded by each pipeline on the given datasets, we have developed a generic, reproducible Snakemake¹³ workflow, which conducts all steps from downloading benchmark data, preprocessing, comparison with a known ground truth, plotting, and automatic deployment of the required software stacks via Snakemake’s Conda/Mamba [9] integration: https://github.com/snakemake-workflows/dna-seq-benchmark. The workflow comes with predefined standard datasets like CHM-eval and GIAB, but can be additionally configured to use any other DNA-seq-based benchmark dataset consisting of a known set of true variants, confident regions where the reported true variants are considered to be complete (i.e. every non-variant position is assumed to homozygously have the reference allele), raw read data (as FASTQ files), and (optionally) sequenced target regions (e.g. in case of exome sequencing). The workflow uses BWA-mem,¹⁴ Picard tools [10] and Mosdepth¹⁵ for calculating the read coverage across the genome. We use Bedtools¹⁶ to limit the known true variants to the confident regions provided by the respective truth publishers and to stratify variants by coverage (see below). For interactive exploration of the results, we use Datavzrd [11] and Vega-Lite.¹⁷ The matching of calls and true variants in a haplotype-aware manner happens via RTG-tools vcfeval [12]. To ensure a fair and correct comparison of the different evaluated callsets, several key points had to be considered, which we outline below.

Read depth stratification and selection of regions of interest. The available read depth can naturally affect both the precision and recall of a pipeline. Hence, the read depth characteristics of a benchmark dataset can have an impact on the derived precision and recall, which can limit the generalizability of obtained results. In order to avoid this effect, we decided to stratify recall and precision by read depth. For any benchmark dataset, this workflow generates a quantized set of regions with low (0-9), medium (10-30), and high ($>30$) read depth using Mosdepth, while considering only reads with mapping quality (MAPQ) $\ge 60$. Notably, this means that, for example, regions in the low read depth category have either only few reads or a lot of reads with uncertain alignments (high mapping uncertainty). We intersect these regions with the confidence regions of the benchmark sample (e.g. as provided by GIAB) using Bedtools. If the given dataset was generated using a capturing approach (e.g. exome sequencing) we further restrict the regions to the captured loci according to the manufacturer. Afterwards, any given callset is split into three subsets with low, medium, and high coverage using Bedtools.

Separating genotyping from calling performance At decreasing read depth or increasing mapping uncertainty, one can expect a callset to yield a decreasing recall: with less evidence, it will become harder to find variants. This is true for both genotyping (i.e. requiring that the variant caller detects the correct genotype) as well as when just requiring the variant allele to be correctly recognized without considering whether the variant is predicted to be homo- or heterozygous (i.e. plain variant calling without genotyping). In contrast, a variant callset’s precision should ideally remain constant and unaffected by a decrease in read depth or increase in mapping uncertainty, if the method manages to correctly report the increasing uncertainty with decreasing depth or increasing mapping uncertainty. The latter behavior differs between measuring a callset’s genotyping or calling precision. In order to make these differences visible, we therefore decided to calculate precision and recall for both genotyping and calling separately.

Variant atomization Some variant callers report complex variants as replacements of longer alleles (i.e., both the reported reference and the alternative allele are longer than one base, e.g. ACCGCGT>ACGCT). While this is in general a good idea (e.g. in order to be able to properly assess the combined impact on proteins), we found this to introduce problems with vcfeval’s internal comparison approach. This resulted in spurious false positives and false negatives in callsets having such variants. Similar to the approach implemented in the hap.py pipeline [13], we solved this issue by introducing a normalization step prior to vcfeval into our analysis workflow, which uses Bcftools¹⁸ to normalize variants, in a way that indels are moved to their left-most possible location, and complex replacements are split into their atomic components—i.e. single nucleotide variants (SNVs), insertions or deletions (indels)—while removing exact duplicates resulting from the atomization.

Reporting For reporting results, we employ Datavzrd to create interactive tabular reports for recall and precision, as well as individual false positive and false negative variants. Datavzrd enables us to just provide the required data as TSV or CSV files combined with a configuration file that defines the rendering of each column. For the latter, one can choose from automatic link-outs, heatmap plots, tick plots, bar plots, or custom complex Vega-Lite plots (which can also be used to define alternative visualizations for an entire table view). For the former, we report a table containing for each callset and each read depth category (low, medium, high) precision and recall (while ignoring whether the genotype was predicted correctly), the underlying counts of true positives (TP), false positives (FP), and false negatives (FN), as well as the fraction of wrongly predicted genotypes. It is important to note that it cannot be excluded that the same variant in the truthset is predicted multiple times by a callset, e.g. as part of several complex replacements (see “Variant atomization” above). We therefore report two TP counts ${\mathrm{TP}}_{\text{query}}$ (number of TPs in the callset with the same matching variant from the truth potentially counted multiple times) and ${\mathrm{TP}}_{\text{truth}}$ (number of variants in the truth set that occur in the callset, each variant counted once, regardless how often it occurs in the callset), with ${\mathrm{TP}}_{\text{query}}\ge {\mathrm{TP}}_{\text{truth}}$. Following the established definitions, precision is then calculated as

An example can be seen in Figure 2. For reporting of individual FP and FN variants, we provide a Datavzrd table view for each that has one row per variant and a column for each callset. In order to visualize systematic patterns arising from the properties of callsets (e.g., using the same variant detection or mapping method), any kind of property can be annotated as a so-called “label” when registering a callset for evaluation with the pipeline. The labels are displayed using a categorical color coding in the header of the table views. Moreover, we perform a ${\mathrm{Chi}}^2$test for the association of the FP or FN pattern of each variant against the different labels in order to detect systematic effects. The variant/label combinations for which this test yields a significant result are then displayed in a separate table view for each type of label. These allow, for example, to spot variants that only occur when callsets use a particular variant caller. Thereby, significance is determined by controlling the false discovery rate over the p-values of the ${\mathrm{Chi}}^2$test using the Benjamini-Yekuteli procedure, as the variants could be both positively (e.g., being on the same haplotype) or negatively (e.g., being on different haplotypes) correlated. In order to combine the results with data provenance information we include the Datavzrd views into a Snakemake report [14], which automatically provides a menu structure for navigation between views, association with used parameters, code, and software versions as well as runtime statistics.

Continuous public evaluation

A central goal of the project was not to conduct a single benchmark and just publish the results, but rather provide a resource for continuous repeated and always up-to-date benchmarking, that is moreover open to any kind of contribution (callsets and code improvements, among others) from outside collaborators. In order to achieve this, we have developed the following approach (see Figure 1 for an illustration). We deployed the benchmarking workflow [15] as a module [16] into another Snakemake workflow that in addition has the ability to download callsets from Zenodo, using Snakemake’s Zenodo integration [17]. Then, we deployed this workflow into the GitHub repository [18] and configured GitHub Actions [19] to continuously rerun the workflow upon every commit on the main branch or any pull request [20]. In order to ensure that the workflow runs sufficiently fast (GitHub Actions offers only limited runtime and resources per job), we have precomputed benchmark dataset-specific central intermediate results (read depth and confidence derived stratification regions) that are computationally intensive to obtain, and deployed them along with the workflow code into the GitHub repository.

Figure 1. Continuous evaluation and reporting workflow.

Upon pull requests or pushes, a GitHub Actions workflow is triggered. This downloads data, runs the Snakemake-based evaluation pipeline, creates the Snakemake report and uploads it as an artifact. If the workflow is triggered on the main branch, its finalization triggers a second Github Actions workflow that builds and deploys the homepage at https://ncbench.github.io.

Upon each completion of the evaluation pipeline, a Snakemake report [21] is generated. In case of pull requests (e.g., contributing a feature or a new callset), the report is uploaded as a GitHub artifact [22], for inspection by the pull request author and the reviewer. In the case of the main branch, we utilize GitHub Actions to trigger the execution of a secondary GitHub Action pipeline in a repository that hosts the NCBench homepage [23]. This pipeline fetches the latest report artifact associated with the main branch and deploys it to the homepage. This way, the most recent results are automatically accessible on the homepage.