Aggv3 site QC¶
For all variants on autosomes and chrX, we applied a number of quality control (QC) metrics. All the metrics are reported in the INFO column in a separate siteQC VCF file corresponding to each shard and subshard. The first iteration of the siteQC VCFs do not include a FILTER column based on these QC metrics, but this will be updated in the future. We do not remove any variants from these files.
The biallelic VCF files contain a FILTER field provided by Illumina. The sites can fail on two conditions:
LowGTR- variants with genotyping rate <0.9-
LowMLSQ- variants with Machine Learning Site Quality < 0.1. The MLSQ takes into consideration genotyping rate and genotyping quality from all samples. This filtering is applied only to variants passing the LowGTR filter.We have observed that variants and sites have been flagged by the
LowMLSQfilter which are seemingly fine in terms of missingness, depth, and quality. This becomes apparent for some multiallelic sites where the most common variant may be flagged with this filter while the less common variant receives aPASSfilter even though the site metrics and quality are largely similar if not identical. This may be caused by the model which was trained on UK Biobank data which would have different characteristics than Genomics England data, hence this 'over-filtering' was not observed in the UK Biobank data. We are working with Illumina to better understand this, and to provide an alternative solution or strategy for filtering (April, 2026).
A file containing paths to the siteQC VCFs is available in BED format at s3://512426816668-gel-data-resources/dragen3.7.8/AggV3_resources/manifests/site_qc/2026-01-06/siteqc_shards.bed. For guidance on working with these VCFs, please refer to our code book.
Site QC of autosomes¶
To calculate the SiteQC metrics, we first decomposed the multiallelic VCF (Processing of Multiallelic VCFs), and calculated the following per-variant quality metrics. They can be found in the INFO field of site QC VCF files.
| Metric TAG | Description | Further description |
|---|---|---|
MEDIAN_DP |
Median DP calculated across all samples. | |
MEDIAN_GQ |
Median GQ calculated across samples with complete genotypes. | The median GQ was calculated from GTs in which fully missing genotypes were filtered out. Values are capped at 99. The filtering method was equivalent to: bcftools query $infile -e "GT=\\".\\"" [...] |
MISSINGNESS_RATE |
Fraction of samples with fully missing genotypes. | Missingness is calculated as: missing_genotypes/(complete_genotypes+missing_genotypes). Missing genotypes: bcftools query $infile -e "GT=\\".\\"" [...]. Complete genotypes: bcftools query $infile -e "GT!=\\".\\"" [...] |
AB_RATIO |
For each heterozygous genotype, a binomial test is performed on reads supporting the called alleles (ref/alt or alt1/alt2). AB_RATIO is the proportion of heterozygotes passing the test (p>0.01) out of all heterozygotes. | The AB ratio is a measure of the evidence supporting whether a heterozygous call is correct or not. This is achieved by testing the distribution of reads supporting ref and alt alleles for each genotype. We apply a binomial test with a strict threshold of p-value >0.01. The ratio is number of heterozygous calls passing this test divided by the total number of heterozygous calls for a variant. For heterozygotes containing a reference allele, LAD[0] and LAD[1] are used for the binomial test, where LAD[0] refers to the REF, and LAD[1] refers to the ALT allele. For heterozygotes where there is no reference allele (e.g., with a 1/3 genotype in the multiallelic VCF), we calculate the binomial test on LAD[1] and LAD[2]. This ensures that AB ratio is calculated using the correct numbers even if the sites are represented as 0/1 and 1/0 in the biallelic VCF (the LAD retains 3 values). |
AN |
Total number of alleles in called genotypes calculated across all samples. | These values are all calculated using the BCFtools plugin fill-tags. |
AC_Hom |
Allele counts in homozygous genotypes calculated across all samples. | As above. |
AC_Het |
Allele counts in heterozygous genotypes calculated across all samples. | As above. |
AC_Hemi |
Allele counts in hemizygous genotypes calculated across all samples. | As above. |
AF |
Allele frequency calculated across all samples. | As above. |
Site QC of chrX¶
Sex chromosome QC was handled in a similar manner to autosomal QC, however some metrics were calculated separately for XX and XY samples. The table below provides a summary of the sex-specific metrics and their names. Sex chromosome karyotype was determined on the basis of DRAGEN 3.7.8 ploidy estimation.
| Metric type | Metric | Metric TAG |
|---|---|---|
| Metrics calculated on all samples | AN, AC, AC_Hemi, AC_Het, AC_Hom | AN, AC, AC_Hemi, AC_Het, AC_Hom / unchanged |
| Metrics calculated on XX and XY samples separately | AN, AC, AC_Hemi, AC_Het, AC_Hom, MISSINGNESS_RATE, MEDIAN_DP, MEDIAN_GQ | e.g. AN_XX or MEDIAN_DP_XY |
| Metrics calculated only on XX samples | AB_RATIO | AB_RATIO_XX |
Site QC of chrY¶
Variant-level QC metrics are challenging for chromosome Y because standard QC methods (e.g., Hardy–Weinberg equilibrium, heterozygosity-based filters, allele balance expectations) assume diploidy and are not directly applicable to a largely haploid chromosome. In addition, the Y chromosome contains extensive repetitive and multi-copy regions that introduce systematic mapping ambiguity, which is not fully captured by generic site-level QC metrics such as depth or variant quality scores.
Reliable Y chromosome analysis typically requires ad hoc and analysis-specific QC, including restriction to high-confidence regions, removal of heterozygous calls, and context-dependent thresholds for coverage and mapping quality. Because appropriate variant filtering on chromosome Y depends strongly on the downstream application, we recommend performing variant QC in an analysis-specific manner. Providing chromosome Y QC metrics is not on our product development roadmap.
Hardy-Weinberg Equilibrium¶
Hardy-Weinberg Equilibrium (HWE) mid-p-value -log10 scores were calculated using PLINK2's --hardy command with the midp and log10 modifiers. The scores are reported in a negative log10 form to preserve the accurate representation of very small p-values. Calculations were performed separately within each inferred unrelated super-population from samples with XX or XY karyotypes. Only the NEG_LOG10_MIDP_HWE value from PLINK2's output is retained in the final dataset and renamed accordingly.
The results are in a tsv format with the following structure:
| CHROM | POS | REF | ALT | ID | neg_log10_midp_hwe_AFR | neg_log10_midp_hwe_AMR | neg_log10_midp_hwe_EAS | neg_log10_midp_hwe_EUR | neg_log10_midp_hwe_SAS |
|---|---|---|---|---|---|---|---|---|---|
| chr1 | 100000 | A | G | chr1:100000:A:G | 0.0996516 | 1.72044 | 2.58358 | 0.171292 | 0.171292 |
For non-PAR regions of chromosome X, HWE p-values were calculated using XX karyotypes only, while incorporating allele frequency information from XY karyotypes. This approach is implemented in PLINK2 and described by Graffelman and Weir (2016). Pseudoautosomal (PAR) regions of chromosome X are treated as autosomes for HWE calculations.
Dataset location¶
The HWE results are currently stored separately to the other siteQC results. A file containing paths to the HWE tsv files is available in BED format at s3://512426816668-gel-data-resources/dragen3.7.8/AggV3_resources/manifests/hardy_weinberg_equilibrium/2026-09-26/hardy_weinberg_equilibrium_shards.bed and in the shard lookup tool. For guidance on working with these files, please refer to our code book.
Mendelian Inconsistencies¶
Mendelian inconsistencies (MIs) were calculated using bcftools +mendelian2 with the --mode a option. We additionally calculated the per-trio counts of MIs using bcftools +mendelian2 with the --mode c option (see the plugin's documentation). For both modes we used the biallelic msVCFs as an input. We chose bcftools +mendelian2 over PLINK2 because we observed issues using PLINK2 where genotypes within trios were being mismatched in particular repetitive regions.
MI data were only calculated for autosomal and chrX subshards. For each of these we provide:
- Annotated site VCF file (with the number of MIs for each site in the new INFO/MERR field without sample-level or genotype information),
- Per-site TSV file with the same information (with columns "CHROM POS REF ALT ID MERR"),
- Per-trio TXT file with the number of MIs per trio.
Trio selection¶
We generated a curated set of ~16,700 family trios (Proband, Mother, Father) to use as an input to calculate MIs using the following selection criteria:
- AggV3 samples were filtered for "complete" referrals/families (i.e. 3 members where there is a Father/Mother/Proband).
- Each individual sample can only be present in a single trio.
- Where multiple trios were present in a family (e.g. Mother, Father, Proband, and Mother, Father, Sibling), only the trio containing the proband was kept.
- In NHS-GMS parents often have two associated samples included, and because of that the parents could technically be present in multiple trios under different sample IDs. However, due to our sample deduplication approach in AggV3 there will only be one sample for either parent resulting in only one of the trios being selected.
- Trios where any member is not XX or XY were removed.
- Cases of suspected uniparental disomy (UPD) were removed. This data has been derived from our healthcare pipelines.
- Trios where full kinship could not be established using our inferred relatedness matrix were removed (see Outliers section below).
Outliers¶
We initially screened the trio list for outliers by only applying the above criteria 1-3 where an outlier was defined as a trio having a count of MIs beyond ±SD*4 from the mean. During this analysis <10 trios were flagged with an excess MIs. Deeper analysis using kinship information from our Population Structure and Relatedness data revealed that kinship was not correctly assigned (that is potential parent-proband swap) for all of these trios. We therefore applied criterium 4 to ensure that this is systematically captured. From our screening using criteria 1-3 we also observed <10 trios as having a number of MIs between a mean ±SD*2 and mean ±SD*4. For these cases, we observed no obvious issues, nor from our internal healthcare findings which could explain any excess in MIs and have left these trios in the final list of trios. In summary, after applying criteria 1-4 no outliers have been removed.
Dataset location¶
The MIs dataset is currently stored separately to the other siteQC results. File containing paths to the Mendelian inconsistencies per-site VCF and TSV files, as well as per-trio TXT files, are available at s3://512426816668-gel-data-resources/dragen3.7.8/AggV3_resources/manifests/mendelian_inconsistencies/2026-09-24/mendelian_inconsistencies_* and in the shard lookup tool. For guidance on working with these files, please refer to our code book. The list of all trios used in the analysis is available at s3://512426816668-gel-data-resources/dragen3.7.8/AggV3_resources/mendelian_inconsistencies_resources/2026-09-24/aggv3_trio_24092026.fam.