Flexible read-aware genotype imputation from sequence using biobank sized reference panels

Domain: biology
Published: 2025-12-13
Editorial outcome: minor_revision

Abstract
Inexpensive and accurate genotyping methods are essential to modern genomics and health risk prediction. Here we introduce QUILT2, a scalable and read-aware imputation method that can efficiently use biobank scale haplotype reference panels. This allows for fast and accurate imputation using short reads, as well as long reads (e.g. Oxford Nanopore Technologies (ONT) 1X, r2=0.937 at common SNPs), linked-reads and ancient DNA. In addition, QUILT2 contains a methodological innovation that is designed to enable imputation of the maternal and fetal genome using cell free non-invasive prenatal testing (NIPT) data. Using a UK Biobank reference panel and simulated NIPT data, we see accurate imputation of the mother (0.25X, r2=0.966, common SNPs) and modest imputation of the fetus (0.25X, r2=0.465, fetal fraction of 10%) at low coverage, with fetal imputation accuracy rising with coverage (4.0X, fetal r2=0.894). We show using simulated data that this could enable both GWAS and PRS for the mother and fetus, which could create clinical opportunities, and if phenotypes can be collected alongside clinical NIPT, the potential for large GWAS.

Introduction
Low-coverage whole genome sequencing (lc-WGS) was originally demonstrated to be a cost-effective alternative to DNA genotyping arrays for genotyping of modern genomes in 20121. However, early use of lc-WGS used statistical methods for phasing and imputation originally designed for microarrays2. Subsequently, efficient dedicated methods for imputation from lc-WGS were designed, which could operate both with or without reference panels, which were generally shown to be more powerful than approaches using arrays3–6. Additionally, imputation with reference panels has been demonstrated as a reliable approach for ancient DNA (aDNA) studies7. Leveraging large reference panels, aDNA imputation can be further improved6. Furthermore, human genetic studies have begun using long-read sequencing more and more for investigating disease-associated variants, and a strategy of mixing high and low coverages is favored for population-scale study8.

Previously, we developed QUILT for rapid genotype imputation of lc-WGS data with a reference panel4. Around the same time, Rubinacci et al. GLIMPSE, another method for genotype imputation4,5. Both QUILT and GLIMPSE share a core two-stage imputation framework, where per-sample imputation leverages a small set of conditioning reference haplotypes. For haplotype matching, QUILT employs the Li and Stephens model9, which is more robust, whereas GLIMPSE relies on the positional Burrows-Wheeler transform10 (PBWT), which offers computational efficiency. QUILT and GLIMPSE further differ in how they treat input data, with GLIMPSE using per-SNP genotype likelihoods, while QUILT uses and phases the sequencing reads directly, which is preferable for linked reads or long-read data. Overall, QUILT is at least as accurate as GLIMPSE and more robust, however GLIMPSE tends to be faster5. Since the publication of these methods, reference panels have expanded significantly in both the number of haplotypes as well as variant density, for example many studies now use the UK Biobank whole-genome sequencing (UKB WGS)11. To accommodate these growing panel sizes, the authors of GLIMPSE developed GLIMPSE26. GLIMPSE2 focuses on improving speed and facilitating biobank sized reference panels by incorporating more sparse data structures for both speed and memory efficiency, as well as the introduction of a dedicated version of PBWT, termed sparse PBWT.

For human GWAS, studies have estimated sample sizes in the millions are required to robustly detect signals that explain most of the heritability underlying common traits and diseases12. This problem is only exacerbated when we consider the need for GWAS in many different populations, to enable population specific effect size estimation, and to help realize the promise of precision medicine for all13. One option is to utilize DNA and phenotypes from pregnant women. Electronic health records (EHR) usually store rich phenotypes for pregnant women, which could include quantitative traits like height, weight, hyperglycaemia14, as well as socioeconomic traits like educational attainment, as well as disorders with age of onset before maternal age or during pregnancy, such as neuropsychiatric conditions like schizophrenia, as well as gestational diabetes mellitus15 and intrahepatic cholestasis of pregnancy16. A natural way to collect DNA from pregnant women is using cell free DNA from non-invasive prenatal testing (NIPT), which is a sensitive and specific screening technology which tests for the common fetal aneuploidies of trisomies 13, 18 and 21, and that, according to the American Society of Obstetricians and Gynecologists (ACOG, 2020), should be offered for all pregnancies regardless of risk. In some countries, NIPT is routinely offered to all pregnant women, with recent uptake rates at around half the population or more, with an uptake rate of 79% and 46% for Belgium and Netherlands respectively17. In 2018, it was suggested that NIPT had already been performed on ten million pregnant Chinese women18. Further, NIPT has already been used as a successful means of conducting GWAS, with NIPT based GWAS discovering novel associations for intrahepatic cholestasis of pregnancy16, gestational diabetes mellitus15, metabolites19, glycemic traits20 and neonatal phenotypes21. Separately, recent studies proposed high-resolution noninvasive prenatal screening to investigate the entire fetal exome from circulating cfDNA by whole exome sequencing (WES)22,23, which usually offers genome-wide off-target reads for non-coding regions for 1×−2× coverage24. The increased benefit of such a test over screening only for aneuploidies should increase the rate of adoption of such tests worldwide. However, state-of-the-art imputation methods for lc-WGS are designed for diploid samples, and there is no dedicated method designed for NIPT data taking both the fetus DNA and mother DNA into account. Given different time points of sequencing cfDNA in the blood of pregnant women, the fetal fraction (FF) varies between 4% and 50%23, and at the most commonly assayed gestational ages of 10 to 20 weeks, is about 10–15%25. Without dedicated imputation methods for NIPT data, imputation accuracy would be negatively correlated with fetal fraction in cfDNA26. As such, the underlying data suggests that imputation from NIPT can be improved, and the scale of NIPT suggests it would be very valuable27,28. Furthermore, imputing fetal genomes before childbirth may facilitate health risk management through polygenic risk scores (PRS) for perinatal traits in clinical practice.

In this study, we present QUILT2, a novel scalable method for rapid phasing and imputation from lc-WGS reads and cfDNA using very large haplotype reference panels. QUILT2 contains three key innovations, two technical and one methodological, compared to QUILT (or QUILT1). First, we introduce a memory efficient version of the positional Burrows-Wheeler transform (PBWT)10, which we call the multi-symbol PBWT (msPBWT). QUILT2 uses msPBWT in the imputation process to find haplotypes in the haplotype reference panel that share long matches to imputed haplotypes with constant computational complexity, and with a very low memory footprint. Second, we introduce a two stage imputation process, where we first sample read labels and find an optimal subset of the haplotype reference panel using information at common SNPs, and then use these to initialize a final imputation at all SNPs. This both speeds up imputation as well as decreases the amount of RAM required. Finally, we introduce a novel methodological innovation to QUILT2, which we call the QUILT2-nipt method, which assumes the observed sequencing reads are from the three haplotypes present in NIPT data, and which imputes both the maternal and fetal genomes using DNA from pregnant women. In what follows, we evaluated the accuracy, performance and scalability of QUILT2 for various data types on diploid samples. In addition, we evaluated the accuracy of imputing the maternal and fetal genomes from simulated NIPT data. Finally, we show how this data can enable both GWAS and PRS, in both the mother and the fetus.

Results
Model overview
QUILT2 is based upon and improves our previous model QUILT1 for diploid genotype imputation, which uses an iterative approach, as follows. First, Gibbs sampling is performed to generate a posterior estimate of the read labels (an estimate of what haplotype each read comes from) given the sequencing reads and the haplotypes in a small haplotype reference panel using the Li and Stephens algorithm9. Second, the small haplotype reference panel is updated based on the current estimate of the haploid dosages, again using Li and Stephens. Due to the second stage in haplotype selection, QUILT1 has linear time complexity in the number of haplotypes and variants in the full reference panel. For RAM and speed reasons, in QUILT1, we disable transitions except between every 32nd pair of SNPs, and store haplotypes as 32 bit integers (or a further representation of this). We refer to 32 consecutive SNPs as a grid.

With QUILT2, we first addressed the scaling of QUILT1. The use of PBWT is common for haplotype identification in genotype phasing and imputation5,6,29–31, including with GLIMPSE and GLIMPSE2, and derivatives exist that work on non-binary symbols, such as multi-allelic and Syllable PBWT32,33. We similarly developed our own analogue of PBWT, which we call the multi-symbol PBWT (msPBWT), which naturally operates on the symbols in the grids used in QUILT1 and QUILT2 (Supplementary Fig. 1). With msPBWT, and using pre-computed indices of the reference panel, we use an algorithm to identify haplotypes in the full reference panel that have long identical matches against a target imputed haplotype with computational time independent to the panel size. Next, we note with WGS derived reference panels, the vast majority of SNPs are very rare, and as such, will slow down methods that are linear in the number of SNPs. In addition, these rare SNPs are of marginal importance in determining which reference haplotype to copy from when there are centimorgan or longer haplotype matches between target and reference haplotype. As such, several methods for genotype phasing and imputation differentiate between common and rare SNPs, for both speed and RAM reasons, including GLIMPSE2 and other methods29,31,34. For QUILT2, we introduce such a two step process, where we iteratively perform imputation using only common SNPs and the reads that intersect them, to get an optimal subset of the full reference panel (Supplementary Fig. 2). We then load all reads at all sites, initialize their phase according to the current best estimate of the underlying haplotypes, and perform a final round of Gibbs sampling, to get final read labels, and from them haplotype phase and genotype dosage. Finally, methodologically, QUILT2 introduces a novel method for jointly imputing the mother and fetus utilizing the lc-WGS cfDNA data from NIPT (Fig. 1). This novel method includes changes to the mathematics as compared to the original diploid model, for instance the probability of the reads given the read labels and the haplotypes, as well as the probability of the read labels given the fetal fraction, to accommodate the three haplotypes present in different frequencies in cfDNA. Further, the change in the prior on read labels, which goes from uniform to unequal and dictated by the fetal fraction, was incorporated throughout the modeling, which necessitated substantial changes to the heuristics to minimize the Gibbs sampling getting stuck in local extrema. Further details, including forms PBWT, common variant initialized imputation, NIPT, as well heuristics, are described in Methods and Supplementary Note.

Diploid imputation performance
In the benchmarking of diploid imputation performance, we predominantly analyzed European samples in the 1000 Genomes. To reduce the computational burden, we analysed only chromosome 20, and unless noted, focus on the analysis with three different reference panels: the 1000 Genome Project35 (KGP, N (haps)=5,008, M (SNPs) = 2,115,074), the Haplotype Reference Consortium36 (HRC, N=54,330, M=884,932), and the UK Biobank 200K WGS genomes11,37 (UKB, N=400,022, M=14,075,021). We first investigated the imputation performance of QUILT2 and GLIMPSE2 for diploid genomes using various data types with multiple sequencing coverages. Throughout the results, unless noted, we stratified accuracy of European samples by allele frequencies from the separate Genome Aggregation Database (gnomAD v3.1.2 Non-Finish European). As exemplars, we used very rare to refer to SNPs with MAF of 0.01-0.02%, rare to refer to SNPs with MAF of 0.1–0.2%, and common to refer to SNPs with MAF of 10-20%.

We first considered imputation using short read Illumina data from the CEU population of the 1000 Genomes. As shown in Fig. 2a and Supplementary Fig. 3, QUILT2 was more accurate than GLIMPSE2 at sequencing coverage ≤0.5× regardless of the reference panel, while GLIMPSE2 was more accurate than QUILT2 for rare variants at coverage ≥2.0×. For higher coverage samples ≥1.0×, at common SNPs, accuracy was similar between the panels and between the methods, while at rarer SNP, far more rare SNPs were imputed, and more accurately imputed with the larger panels (Supplementary Fig. 3, Supplementary Tables 1, 2). For imputation of modern samples with smaller panels (KGP and HRC), we found that QUILT2 is less robust and accurate than QUILT1 and GLIMPSE2 for rare variants with default settings (Supplementary Fig. 3). For QUILT2 with the big UKB panel, we also evaluated the effect of the msPBWT haplotype selection versus the Li and Stephens approach of QUILT1, and note the msPBWT approach generates slightly more accurate results (Supplementary Fig. 4). In terms of phasing, the larger UKB panel substantially improved QUILT2 phasing accuracy compared to the smaller panels (Supplementary Fig. 5, 6), and mirroring the imputation results, phasing with QUILT2 was more accurate than GLIMPSE2 particularly for lower coverages.

We also assessed the imputation performance on aDNA from Afanasievo culture (~4.6 kya) with three reference panels. We followed the Mapache38 pipeline used in the GLIMPSE2 paper, which implements the practice by Mota et al. 20237. Interestingly, QUILT2 outperformed GLIMPSE2 across all frequency variants using KGP and HRC panel, while both performed similarly using the UKB panel (Fig. 2b, Supplementary Fig. 7). Similar to results for modern DNA, the imputation of aDNA showed similar accuracy for common variants across reference panels, while the UKB panel again improved the performance of rare variants.

Lastly, we evaluated the imputation accuracy using long sequencing reads from Oxford Nanopore Technologies (ONT) by calculating both the r2 at variant level and F1-score at sample level. To investigate the gain of using the UKB panel for long-read imputation, we included QUILT1 with HRC for comparison. We found that QUILT1-HRC performed better than QUILT2-UKB for coverage <1.0x but with fewer variants imputed (Fig. 2c). Additionally, we evaluated QUILT2 and GLIMPSE2, the inexpensive and scalable methods that can use the UKB panel. Across all coverages and reference panels, QUILT2 was substantially more accurate compared to GLIMPSE2 for long reads (Fig. 2c, Supplementary Table 3), e.g., common r2=0.937 vs r2=0.695 at 1× coverage.

Scalability and computational cost
We next evaluated the computational cost and scalability of QUILT2 along with GLIMPSE2 for both speed and memory usage, varying both the size of the haplotype reference panel, as well as the number of target samples (with 1x coverage) to run. As shown in Fig. 3a, QUILT2 has approximately constant runtime with respect to haplotype reference panel size, and the per-sample runtime is not affected appreciably by the input sample size. By contrast, the per-sample runtime of GLIMPSE2 is affected by the number of samples run, and only shows approximately constant computational complexity with respect to haplotype reference panel size for larger number of input samples (N>100) but not for smaller size, as has been previously noted by the authors of GLIMPSE26. In the extreme case, when only imputing one sample at a time with around 1 million haplotypes in the panel, QUILT2 would be an order of magnitude faster than GLIMPSE2. Further, QUILT2 with both msPBWT and the two-stage imputation enabled was ~3× faster and used ~4× less RAM, than QUILT2 with only msPBWT (Supplementary Fig. 8). Comparing QUILT2 and GLIMPSE2 for RAM, we see that both approaches show linear increases in RAM with reference panel size (Fig. 3a). Given around 1 million reference haplotypes and an average chunk size of around 5 Mb and 500 Kb buffer to impute on chromosome 20, the maximum RAM for both is less than 10 GBs. RAM usage is approximately constant for GLIMPSE2 with sample size, whereas it increases in QUILT2 (Fig. 3). In addition, QUILT2 also runs faster for lower coverage data (Supplementary Figs. 8, 9). Since the UKB panel is only accessible on the UK Biobank research platform (RAP), we also assessed the cost of imputing whole genomes on the RAP. Since both QUILT2 and GLIMPSE2 do not borrow information across samples, and as the cheapest way to run jobs on the RAP is using low or normal priority jobs which can be interrupted6, it is advantageous to impute large datasets using small batches, and then to combine results afterwards. In their previous work, the authors of GLIMPSE2 showed that optimal per-sample runtimes were achieved at approximately 100 samples6. Comparing per-sample costs for QUILT2 (using N=32) and GLIMPSE2 (using N=100), we see both methods are inexpensive, though QUILT2 (£0.237, £0.299, £0.371 and £0.475 for 0.1×, 0.5×, 1.0× and 2.0×, respectively) is more expensive than GLIMPSE2 (approximately £0.09 for all coverages). We note these are idealized estimates, and real costs will depend on prices at runtime, server congestion, etc, and furthermore, that default parameter settings for QUILT2 favour accuracy over speed.

Imputation performance on simulated NIPT data and one real NIPT sample
We next investigated the QUILT2-nipt method for maternal and fetal genome imputation. We simulated 30 NIPT samples by mixing the sequencing reads of the mother and child in the 30 CEU trios from KGP (Methods). Here, we focus on NIPT imputation using the UKB panel across different sequencing coverages and various fetal fractions (FF). Firstly, we compared the QUILT2-nipt method with the existing diploid method (i.e. QUILT2-diploid) for maternal genotype imputation. As shown in Fig. 4a and Supplementary Table 4, given FF=0.1, there is a clear gain of accuracy in QUILT2-nipt over QUILT2-diploid (0.25×, very rare r2=0.724 vs 0.695; 1.0×, very rare r2=0.831 vs 0.799; 2.0×, very rare r2=0.851 vs 0.839). At higher sequencing coverage, for higher fetal fractions, QUILT2-diploid performed increasingly poorly, while QUILT2-nipt still maintained high accuracy for common variants as expected (e.g. 4.0×, FF=0.2, common r2=0.994, 0.853 for QUILT2-nipt and QUILT2-diploid, respectively). In addition to simulations, we also confirmed the effectiveness and improvements of QUILT2-nipt over QUILT2-diploid for maternal genotype imputation with one real NIPT sample (Supplementary Fig. 10).

Next, we investigated the performance of QUILT2-nipt on fetal genotype imputation where there are no other approaches for comparison. As shown in Fig. 4b and Supplementary Table 5, both increasing fetal fraction, as well as sequencing coverage, contribute to improved accuracy, at both rare and common SNPs. For instance, at 0.25× coverage, accuracies at common SNPs improve with higher FF (r2=0.465, 0.657, 0.783 for FF=0.1, 0.2, 0.3, respectively), and they improve further at higher coverage with 4.0x (r2=0.894, 0.944 and 0.955, same FFs). However, we do not observe a gain in accuracy when the sequencing coverage exceeds 4.0×, assuming FF=0.3 (Supplementary Table 5). We note that r2=0.955 appears to be the upper limit for fetal genotype imputation using the UKB panel, versus 0.92 and 0.90 for the HRC and KGP panels, respectively (Supplementary Fig. 11). This contrasts with the performance for the mother, where the upper limit remains 1.0 regardless of the reference panel. This improved performance for the larger panel is driven by the improved phasing capabilities offered by the larger reference panel (Supplementary Fig. 5), enabled by msPBWT within QUILT2, which is especially important for phasing the reads from the paternally transmitted haplotype. In terms of the runtime, where mother and fetus are imputed together, we note that the QUILT2-nipt method is about 30% slower than the QUILT2-diploid method and imputing only common variants can reduce both the runtime and memory substantially (Supplementary Fig. 9).

Power in association testing and prediction with simulated NIPT data
Finally, we explored the usability of imputed results from QUILT2-nipt for both GWAS and PRS. To create realistic NIPT sequencing samples, we used UK Biobank samples, for which we have real phenotypic data, and used parent offspring pairs to represent the mother and the fetus (“Methods”). We used 11,028 parent-offspring pairs and simulated mean sequencing coverages of 0.25×, 1.0×, and 2.0×, with fetal fractions being simulated with mean 10% and following a normal distribution, with range from 5% to 15% (Supplementary Fig. 12). We conducted association tests on chromosome 1 for 25 quantitative traits using a linear mixed model on these 11,028 individuals (Methods). We only tested quantitative traits but expect binary traits to show broadly similar results.

We compared the GWAS association findings using QUILT2-nipt imputed maternal and fetal genotype dosages to a baseline that used imputed genotype dosages from arrays. With the QUILT2-nipt imputed data, we found that at common SNPs the r2 for maternal genotypes were 0.973, 0.994, 1.000 at 0.25×, 1.0×, 2.0× coverage respectively, and the common r2 for the fetal genotypes were 0.480, 0.75 and 0.82 (Fig. 5a). GWAS summary statistics based on different imputed datasets did not show signs of inflation of test statistics as noted by Q-Q plots (Supplementary Fig. 13). Next, we assessed the power of the different coverages and approaches, by measuring how many of 615 independent GWAS signals on chromosome 1 found using the entire UK Biobank participants could be identified (p<0.05) using the imputed datasets of 11,028 simulated NIPT samples (Methods). Using the array-based imputation, 190 associations were found, while for QUILT2-nipt maternal imputation, 190 (100% estimated relative power), 188 (99%) and 181 (95%) were found at 2.0×, 1.0×, and 0.25×, respectively. For QUILT2-nipt fetal imputation, 163 (86%), 159 (84%) and 106 (56%) of loci were validated for 2.0×, 1.0×, and 0.25× data respectively (Fig. 5b). Notably, even though the common r2 for fetal genotypes was only 0.82 at 2.0× with 10% FF, we still obtained a relative power of 86% in the GWAS. This is further validated in an error-prone setting where the estimates of fetal fraction are measured with uncertainty and thus the prior for QUILT2-nipt method is mis-specified (Supplementary Fig. 14).

In addition to GWAS, the fetal imputed genotypes from QUILT2-nipt can be used to generate fetal polygenic risk scores, which can have potential clinical value for perinatal traits. We therefore evaluated the performance of PRS using QUILT2-nipt imputed fetal genotype dosages for 1000 individuals (Methods). We examined 21 traits with incremental r2>0.1 (except for the waist-hip ratio) in unrelated white British Europeans39, which is the difference in r2 between the full prediction model (including PRS) versus the model with covariates alone (age, sex, age×sex, assessment center, genotyping array, 10 PCs). Using the reported PRS weights of each SNP, we calculated the PRS for each individual given different imputed datasets, that is the array-imputed, as well as NIPT imputed from 2.0×, 1.0×, and 0.25× coverage. We compared the incremental r2 of the NIPT imputed datasets to the array imputed datasets, referring to this comparison as the relative incremental r2 (Methods). As shown in Fig. 5c and Supplementary Fig. 15, with effective sequencing coverage for the mother at least 0.2× (e.g. 0.25×, 10% FF), we observed a relative incremental r2 >0.967. For the fetal PRS, while at low effective sequencing coverage, relative incremental r2 is modest (e.g. 0.25×, 10% FF, giving 0.025× coverage, r2=0.494), at high effective sequencing coverage (e.g. 2.0×, 20% FF, giving 0.4×), we observed a high relative incremental r2 >0.935. As expected, the PRS performance is aligned with the imputation accuracy.

Methods
QUILT2
Details about the QUILT model have been previously published4. In brief, QUILT is a Gibbs sampler that samples read labels based on the observed data (i.e. sequencing reads) and parameters of the model (e.g. haplotype reference panel), and from this, allows for phasing and genotype imputation. Here, we provide information about the mathematical changes to QUILT2 that form the basis of the QUILT2-nipt method. We also give a brief introduction to msPBWT, as well as the iterative process that uses common SNPs then all SNPs. More details about msPBWT, as well as the heuristics used in Gibbs sampling for NIPT are given in the Supplementary Note.

In the NIPT mode, we used Gibbs sampling to generate draws from H~P(H|O,λ,FF), where H is a vector of read (haplotype) membership labels, O is our reads (i.e. observed bases at SNPs), FF is the fetal fraction, which is required as an input and can be obtained from NIPT report59,60, and λ is the parameters of our model, including the recombination rate and reference haplotypes. In NIPT, sequencing reads come from three sources, and let us label them as 1=maternal transmitted, 2=maternal untransmitted and 3=paternal transmitted. As such, for each read v, we have the prior probability of haplotype membership as

P(H_v=h_v|FF) = 0.5, h_v=1
0.5 * (1-FF), h_v=2
FF, h_v=3

Let o be the realized observations for the random variable O (sequencing read bases and base qualities). Consider in the Gibbs sampling that we want to sample a new value for read indexed by v with read label h_v, conditional on all other read labels. Let H_v be this random variable at this point in the Gibbs sampler and let H_{-v} be a random variable representing the remaining read labels. We therefore need to calculate:

P(H_v=h_v|H_{-v}=h_{-v}, O=o, λ, FF) = P(O=o, H_v=h_v, H_{-v}=h_{-v}|λ, FF) / sum_{i=1}^3 P(O=o, H_v=i, H_{-v}=h_{-v}|λ, FF)

for h_v ∈ {1,2,3} and sample h_v using this probability. To do this we use P(O=o, H=h|λ, FF)=P(O=o|H=h, λ) P(H=h|λ, FF) which uses the prior probability of all reads, which is the product of the prior for one read described above. This prior probability does not cancel like it did for QUILT1 (diploid) where P(H=h|λ) was constant regardless of h_v. For P(O=o|H=h,λ) we note that given the vector of read labels, we split the reads into three sets reflecting their haplotypic origin, call that O_i, and can generate probabilities in the normal ways for HMMs, in an obvious extension to what is presented in the QUILT1 paper. For G_t^i the haploid (genotype) for haplotype i for some SNP t, we can calculate the posterior probability that this is the alternate allele (encoded by a 1) as follows

P(G_t^i=1|O, H, λ) = sum_{k=1}^K θ_{t,k} P(Q_t^i=k|O_i=o_i, λ)

where, θ_{t,k} is the probability that reference haplotype k carries the alternate allele at SNP t, Q_t^i is the reference haplotype carried at grid g for haplotype i, noting that P(Q_t^i=k|O_i=o_i, λ) is a standard output from an HMM. From these probabilities, a haplotype dosage can be taken as the sum across several Gibbs samplings, and from this, genotype dosages can be taken for the mother by summing up the value for i=1 and i=2, while for the fetus we choose i=1 and i=3.

We have introduced the msPBWT algorithm in QUILT2 to find haplotype matches for target haplotypes against the full reference panel, with computational complexity independent of the size of the full reference panel. An overview of msPBWT is given in Supplementary Fig. 1. First, the haplotype reference panel, containing binary information about whether haplotypes carry reference or alternate alleles, is encoded, and then ranked. This ranked matrix contains columns with 255 symbols representing the most commonly seen 32 SNP haplotype in this window (grid), and a final symbol (0) representing a haplotype not amongst the 255 most commonly seen, for which information (the haplotypes carried and their positions) are stored in helper matrices. Next, two indices are built: A, which allows for indices to be looked up at any grid to put the original haplotype reference panel in reverse sorted prefix order; and U, which allows for determination between grids about how the ordering of indices changes. Both A and U can be used with a two part algorithm that allows for rapid determination of which haplotypes in the haplotype reference panel contain long exact matches to the target haplotype. First, a vector f is found, representing at each grid, where the target haplotype would be inserted into the haplotype reference panel to place the haplotype reference panel in reverse sorted prefix order. Second, scanning up and down is performed, to maintain a list of the closest L haplotypes to the target haplotypes. Finally, filtering is done with respect to a minimum length of match, and when more matches are identified than requested haplotypes, further filtering is done for the haplotypes identified for the two haplotypes, which optimizes both the length of matches, as well as their uniqueness along the target region. Full details of msPBWT are given in the Supplementary Note.

Finally, we note that with high coverage WGS based reference panels, the vast majority of variants are very rare, and thus are not informative early in the iterative process, as they offer little discriminatory information about which reads come from which haplotypes. As such, we introduce a further two-step process, where we first impute using a set of common variants and the reads that intersect them, to get an initial estimate of the carried haplotypes. We then load all reads at all sites, initialize their phase according to the current best estimate of the underlying haplotypes, and perform a final round of read-sampling (Supplementary Fig. 2). Together, these changes made QUILT2 substantially faster and more memory efficient than QUILT1. Finally, we upgraded a heuristic, to perform efficient block Gibbs sampling between all sequential pairs of grids, to identify and solve phase switch errors with read labelling.

Datasets
We created various test and truth datasets to benchmark all methods. We used three reference panels for testing genotype imputation and phasing accuracy, which are the 200K WGS derived UK biobank reference panel, the HRC and the 1000 Genomes Project. Additionally, for benchmarking the computational performance, we used the very big reference panel of N=976,630 haplotypes (488,315 individuals), that is the UK biobank data imputed by the Genomics England haplotype reference panel61, which we referred to as the UKB-GEL. We used liftover to convert the HRC reference panel from the GRCh37 to the GRCh38 build using the Genome Analysis Toolkit62 (GATK) Picard Liftover VCF v2.22.2.

CEU. We downloaded the high coverage (>30x) CRAM files of 30 CEU trios in the 1000 Genome project, from which we called the true genotypes at sites in the UKB-GEL panel using bcftools63 v1.18 with option ‘call -A -m -f GQ -T ref.sites’. We filtered the true genotypes at sites with read depth lower than 10. Leveraging trios information, phasing was done first assuming Mendelian inheritance, excluding triple-heterozygous sites using bespoke R v.4.2.2 code; then, the excluded sites were phased with this scaffold using shapeit4 v4.2.264.

ONT. We downloaded high-coverage ONT alignment files of GRCh38 build (https://s3.amazonaws.com/1000g-ont/index.html?prefix=ALIGNMENT_AND_ASSEMBLY_DATA/FIRST_100/) for 11 samples in the EUR super-population from the first 100 ONT samples in the 1000 Genomes Project65. When evaluating the imputation accuracy with ONT reads, we used the called genotypes from high coverage short Illumina reads as truth, which was done using the same bcftools command.

Ancient DNA. We downloaded four high-coverage ancient DNA samples from the Afanasievo family (mother I3388, father I3950, son I3949, son I6714), who lived ∼4.6 thousand years ago (ka) in the Altai Mountains of Russia, with average coverages of 10.8×, 25.8×, 21.2×, and 25.3×, respectively, which have reliable genotype calls and phasing66. Since both the public BAM and VCF files are of GRCh37 coordinates, we first liftovered the VCF with true genotypes to GRCh38 and converted the downloaded BAM files back to raw FASTQ using ‘samtools fastq’ using samtools v1.18. Then we used the mapache38 snakemake workflow, which implements the best practice for aDNA imputation7, to map the FASTQ files against the human reference genome GRCh38 to generate the BAM files for imputation benchmarking.

NIPT cfDNA. For real NIPT cfDNA data, we found one sample that is publicly available67, which offers both the sequencing data from cfDNA in maternal plasma and saliva at high coverage (20x and 30x, respectively). We downloaded the raw FASTQ file from https://my.pgp-hms.org/profile/hu058D3E. Following best practice67, we called the genotypes of the mother using high coverage saliva DNA with GATK, and we generated the BAM file of GRCh38 coordinates from cfDNA of the pregnant women. Additionally, we simulated thousands of NIPT cfDNA with various coverages using the real sequencing reads from trios or duos individuals. The details are described in the next section.

NIPT cfDNA simulations from real reads
We simulated 3 sets of NIPT cfDNA reads for different analyses using the real sequencing reads. First, to evaluate the QUILT2-nipt and QUILT2-diploid method on maternal genotype imputation, we simulated NIPT cfDNA with various fetal fractions by mixing the maternal reads and offspring reads in the 30 CEU trios from the 1000 Genome Project. We only simulated BAM files for chromosome 20 for benchmarking purposes. Second, to evaluate the GWAS performance on NIPT samples and obtain as many samples as possible, we first identified 1060 trios and 4114 duos in the UK Biobank given the released kinship statistics37. Specifically, we defined the first degree relatives as the pairs with kinship value >0.1767 and IBS0 value <0.0012. To further identify the trios and duos, we defined that the difference between the parents and offspring in age must be greater than 15. Therefore, given a list of trios and duos, we simulated 11,028 NIPT samples by treating the offspring as the parent regardless of sex and vice versa. We only simulated the BAM files of 11,028 samples for chromosome 1 due to the limited budget and computational resources. Third, to evaluate the PRS performance that requires genotypes of the whole genome, we sampled only 1000 NIPT samples from the above identified duos and trios, which includes 500 unique parent-offspring pairs. For those 1000 NIPT samples, we simulated the BAM files for the whole genome.

Imputation benchmarking workflow
To access the performance of QUILT2 and GLIMPSE2 (v2.0.0) on various datasets, we developed a reproducible and configurable snakemake workflow for imputation benchmarking (https://github.com/Zilong-Li/lcWGS-imputation-workflow), which includes downsampling BAM/CRAM files to low coverage; creating reference panels by removing the target samples from the input panel; parallelly running by chunks and ligating results; automated reproducible pipeline with support to High-Performance Computing clusters.

Assessing imputation accuracy
To assess the imputation and phasing accuracy, we developed an efficient function vcfcomp in the vcfppR package v4.6.068, which can calculate various concordance metrics between two VCF/BCF files at either samples level or variants level. In the output VCF files, homozygous reference, heterozygous and homozygous alternative genotypes are coded as integer value 0, 1 and 2 respectively, while the genotype dosage are float values between 0 and 2. There are three concordance metrics defined here. First, we defined r2 as the squared Pearson correlation between the imputed dosages (test) and high-coverage genotype calls (truth). Second, we defined the non-reference concordance (NRC) rate as NRC=1 - (e0+e1+e2) / (e0+e1+e2+m1+m2), where e0,e1,e2 are the counts of mismatches for genotype 0,1 and 2 respectively, while m1 and m2 are the counts of matches between the imputed and truth for genotype 1 and 2. Third, we defined the F1-score as F1=2*TP/(2*TP+FP+FN), where TP is the counts of matches between the imputed and truth for genotype 1 and 2, FP is the counts of mismatches for genotype 1 and 2, and FN is the counts of mismatches for genotype 0 and 1.

To assess genotype imputation accuracy at variant level, we first grouped SNP by its frequency intervals and aggregated the concordance metrics (e.g. r2) for all samples across all SNPs in that frequency interval. On the other hand, to assess the genotype imputation accuracy at the sample level, we only aggregated the concordance metrics (e.g. NRC, F1) for genotypes of a single sample at a given allele frequency interval. The allele frequencies were taken from either the gnomAD v.3.1.2 database (for modern European samples, https://gnomad.broadinstitute.org/downloads#v3) or from the reference panel itself (for ancient DNA).

Assessing phasing accuracy
To assess phasing accuracy, we used phasing switch error (PSE) rates as follows, which is implemented in vcfcomp with option stats=pse. First, consider that we have both true haplotypes (truth) and the imputed haplotype (test) at sites where the truth genotypes are heterozygous. Next, define as discordant any test sites that are also not heterozygous. On the remaining sites, define a phase switch error when either the true haplotypes record a change where the haplotype carries the alternate allele between adjacent heterozygous sites when the test haplotypes do not or vice versa. We removed from consideration sites that were flipped, that is, yielding consecutive phase switch errors. The PSE rate is the number of phase switch errors divided by the total number of pairs of consecutive heterozygous sites examined and can be combined across discrete imputed windows.

Association testing and validations
Since the data consists of the related individuals from families, we used the linear mixed model implemented in gemma69 v0.98.5 for association tests accounting for the genetic relatedness. First, we computed the genetic relatedness matrix (GRM) using ‘gemma -gk 1’ based on the genotypes from the array with the default setting for SNP QC (-maf 0.01 -miss 0.05). Then we ran association tests on multiple imputed datasets in BIMBAM format (for dosages) using ‘gemma -lmm 3’ with the same GRM and sex, age and 10 PCs as covariates. The PCs were calculated based on the imputed array data for the target samples using PCAone70 v0.4.4. To prepare the BIMBAM format for gemma, we used bcftools query -f ‘%CHROM:%POS, %REF, %ALT[, %MDS]
’ to extract the dosages of the mother, and used bcftools query -f ‘%CHROM:%POS, %REF, %ALT[, %FDS]
’ to extract the dosages of the fetus from the VCF output of QUILT2.

We used the publicly available summary statistics of GWAS on more than 330 thousand individuals to find a list of the significant association hits in Europeans, which was downloaded from https://www.nealelab.is/uk-biobank (GWAS round 2). Then, to define the independent association signals that we need to validate using our target samples, we used the extract_instruments function with default parameters from the TwoSampleMR71 v.0.5.10 package to perform the clumping and curate a list of independent association signals for each trait.

PRS and incremental r2
To calculate the PRS, we downloaded the polygenic score of each trait for Europeans39 from the PGS catalog https://www.pgscatalog.org/publication/PGP000332/. We calculated the PRS for each individual using the imputed genotype dosages and effect size of each trait, which is given by

PRS_i = sum_{j=1}^M β^j G_{i,j}

Where PRS_i is the PRS for each individual, β^j is the effect size associated with each SNP, G_{i,j} is the imputed genotype dosage (float value between 0 and 2) for each SNP and each individual. Since the array-imputed data used genome builds GRCh37 while the lc-WGS imputed data used GRCh38, we used the harmonized score files that shared the same variants for both genome builds, which resulted in 1,095,616 out 1,109,311 reported SNPs (98.8% matches). We used the nextflow workflow (https://github.com/PGScatalog/pgsc_calc) to automate our analyses72. Then, we estimated the prediction accuracy using the incremental R2 (IR2), which is the difference in r2 for two regression models for each trait, with one having the PRS and one not having it, with each adjusted for age, sex and 10 PCs.

We ran the regression in R using the lm function. Therefore, we can obtain the IR2 for each trait and each dataset. To estimate the accuracy of PRS using the lc-WGS imputed data compared with the array-imputed data, we defined the relative incremental R2 (RIR2) as follows,

RIR2 = IR2_lc / IR2_array

Where IR2_lc and IR2_array are the IR2 for each trait calculated with the lc-WGS imputed data and the array-imputed data respectively. Further, to estimate the standard error (se) of RIR2, we ran the t.test in R on a vector of RIR2 (21 traits).

Benchmarking and computational cost
To benchmarking the performance of QUILT2 with respect to the size of reference panel, we used the one of the largest reference panel, that is the UK biobank data imputed by the Genomics England haplotype reference panel, which contains N=976,630 haplotypes (488,315 individuals) and M=7,179,683 SNPs for chromosome 20. Then we randomly sampled N=2e4, 2e5, 4e5, 6e5, 8e5 haplotypes to create 6 more smaller panels. To compare QUILT2 against GLIMPSE2, we split the imputation task into multiple chunks using the same chunk size as outputted by the GLIMPSE2_chunk program. Further to fairly measure the average RAM and time for each sample, we ran both software using only one thread, although both support multithreading. To estimate the realistic cost using the WGS-derived UKB-200K panel on the RAP, we maximized the usage of CPUs threads available as possible as we can. The benchmarking was performed on virtual machines of type mem2_ssd1_v2_x4 with 4 cores and 16 GBs memory on the RAP using normal priority. The cost was based on the instance rate card v2 (https://platform.dnanexus.com/resources/UKB_Rate_Card-Current.pdf).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents QUILT2, a genomic imputation tool, demonstrating enhanced scalability, computational efficiency, and accuracy over existing methods. It highlights QUILT2’s performance across diverse scenarios, including short and long sequencing reads, ancient DNA, and non-invasive prenatal testing (NIPT), where it imputes both maternal and fetal genotypes.\nI really appreciate the quality of the work; it’s very well presented in the paper, and the accompanying online resources (like the Snakemake workflow (!), tutorials, and comprehensive documentation) are excellent. Thanks!\nMajor Comments:\nThe authors mention two technical advancements (reduced computational complexity and two-stage imputation). How do these changes to the Positional Burrows-Wheeler Transform (-> resulting in msPBWT) influence the trade-off between speed and accuracy compared to PBWT? Also, what exactly contributes to QUILT2’s improved accuracy, considering that GLIMPSE also uses PBWT (to the best of my knowledge, not stated in the manuscript)?\nIs QUILT2 intended to fully replace QUILT1 for all reference panels / use cases, or are there cases where QUILT1 remains preferable? Maybe this should be stated in the manuscript.\nWhy is the imputation of long reads generally less accurate, as shown in Figure 2? Is this somewhere mentioned?\nThe authors claim that QUILT2 significantly outperforms GLIMPSE2 in scenarios where a single sample is imputed with a large reference panel. Can they elaborate on what makes QUILT2 so efficient in this context?\nWhen using QUILT2 on the UK Biobank Research Analysis Platform (RAP), are additional pre-processing steps needed for the reference panel? (e.g. QUILT2_prepare_reference.R ) Specifically, does the panel need to be in VCF format, and if so, is this format already available on RAP (I guess not)? Are there post-processing steps for QUILT available (e.g. ligate)? I think users would highly appreciate a tutorial how to create reference panels (e.g. locally or on UKB) run the imputation process + post-processing steps. How does the chunk size affect the quality? This could be interesting when conducting a WGS imputation for example on a local Slurm cluster.\nFor the new QUILT2-nipt2 approach, is the fetal fraction required as an input parameter? Is this somewhere stated?\nAncient DNA is briefly mentioned in the abstract but lacks focus in the introduction or in the title. Could the authors provide more context earlier in the manuscript to align the emphasis on ancient DNA with its later discussion? Or to state it differently, is the use case of ancient DNA central to the manuscript?\nGiven the reliance on simulated NIPT data, how closely do these simulations mirror real-world scenarios? Are there potential biases that might affect the conclusions?\nAre there specific genetic architectures or data types where QUILT2 might be less effective compared to other imputation methods?\nThe results suggest that QUILT2-nipt could significantly enhance GWAS and PRS performance using imputed fetal genotypes. Could the authors elaborate on the implications for clinical applications, particularly the accuracy requirements in prenatal risk assessments? Are there recommended r² thresholds for polygenic score calculations based on fetal fraction? Given the ongoing debate around the clinical utility of PGS and the uncertainty of fetal genotypes (depending on seq depth + FF, see also below), this topic could result in diverse opinions.\nMinor:\nThe abstract describes the fetal genome imputation results using QUILT2 at 0.25× sequencing coverage with a 10% fetal fraction as “accurate,” yet the r² is only 0.465. This statement seems potentially misleading—could the authors clarify this assessment?\nLine 154: Please confirm the minor allele frequency (MAF) levels used to define common SNPs.\nThere are several typos in the manuscript that need to be fixed.", "reviewer_2": "Summary of the key results, Originality & significance\nLi et al. reported an improved method for genotype imputation in data produce form short- and long-read sequencing data, named QUILT2. The methodology is primarily based on their previous work (QUILT) that makes use of sequencing reference data to impute undetermined genotypes of shallow-sequencing data (PMID: 34083788). The authors stated two technical and a methodological improvement over their previously published works (PMIDs: 27376236, 34083788). Technical improvements include: (i) a memory efficient PBWT with the purpose of finding haplotypes in the reference panel that shares long matches with constant computational complexity and “very low” memory footprint; (ii) a two-stage imputation approach that first samples read labels and makes use of common SNPs to find an optimal subset of haplotype reference panel. Methodological improvement includes assumes that three haplotypes are present in the NIPT data.\nThis reviewer did not find these improvements conceptually novel relevant to their previous works (PMIDs: 27376236, 34083788) and the existing work GMPSE2 (PMID: 37386250). In additions, although the authors claim applicability of their method (three-haplotype option) for non-invasive prenatal testing (NIPT) by discerning fetal and maternal genomes and imputing the fetal genome, in my view, the assumptions are far from real cfDNA component in blood stream of pregnant women (see below). And that it has been performed on simulated data from whole-genome sequencing data of trios. For the latter, the r2 values are specifically too low to claim that the methodology works on low-pass sequenced NIPT data.\nMajor comments:\n- The QUILT2-nipt method assumes the observed sequencing reads are distributed among three haplotypes (maternal transmitted, maternal untransmitted, and paternal transmitted). However, the complexity of real-world data, including recombination in maternal haplotypes and errors in fetal fraction estimation, may introduce inaccuracies that are not addressed. How do you address recombination sites between the maternal homologs?\n- Given my comment above, it is incorrect to assume that NIP data has “triploid nature”, it is rather tetraploid (prental chromosomes recombine).\n- R2 is very low for real NIPT data. Therefore application of QUILT2 can lead to wrong SNP calls. Could you clarify which parts of genomes cause this low r2 value?\n- The method has higher accuracy for common variants but its applications with (very) rare variants is challenging and substandard. While this is expected, the reliance on filtering out rare variants early in the iterative imputation process could overlook clinically relevant but less frequent variations. Specially in the context of NIPD. Please clarify how you would deal with this challenge?\n- The authors discuss the potential clinical utility of imputed fetal genomes for traits such as gestational hypertension and preeclampsia but does not provide sufficient empirical evidence or validation studies to support these claims. Without these and my concerns above, this reviewer finds such phrases more over-statements.\n- The authors briefly state the potential of QUILT2 to identify structural variations and copy number variants but did not provide how or a benchmarking for these applications. Please clarify.\n- For the diploid imputation performance, the authors indicate that they “used Chr 20 for convinience”. Why did you make such choice? Is this then correct that all the data presented in this manuscript are using Chr 20 only?\n- The authors devised a 32-bit grid strategy and allowed transion after 32nd pair of SNPs. It is not clear to this reviewer why 32 and how the transitions have happened?\n- Description of QUILT2 for (ONT) long-read sequencing sounds preliminary. Please clarfiy how you make use of long-read sequencing data and kb-Mb size haplotypes. How it differs from shor-read sequencing? How it performs for structural rearrangments? How do you potentially detect errornous SNVs of the ONT data?\nMinor comments:\n- In the beginning of the Results section “…QUILT2 is based upon, and improves, our previous model QUILT1 for diploid genotype imputation, which operates, which uses an iterative approach…”. Either use “which operates” or “which uses”\n- Please make figures that illustrate QUILT2, clearer. Specially for the grid search and using common SNPs as anchors.\n- The authors state that “In some countries, NIPT is routinely offered to all pregnant women, for instance in the Netherlands, with an average uptake rate of 46%”. This is not up to date. This uptake rate is from the TRIDENT study that was finalized in March 2023. Since April NIPT in the Netherland is a screening program with very high uptake rate. Please correct.\n- Please insert page and line numbers to make the review process easier and more efficient.", "reviewer_3": "The article presents QUILT2, a method for genotype imputation from low-coverage whole-genome sequencing (lc-WGS), building on the QUILT1 framework. QUILT2 integrates a memory-efficient implementation of the positional Burrows-Wheeler transform (msPBWT) and a two-stage imputation process, aiming to improve the speed and scalability of imputation. Additionally, it introduces a novel approach for imputing both maternal and fetal genomes from non-invasive prenatal testing (NIPT) data. While the method shows promise in some applications, especially for NIPT data, the manuscript raises several concerns.\nBiased benchmark\nOverall, the benchmark puts emphasis on QUILT2 outperforming GLIMPSE2, while it seems to avoid showing scenarios where GLIMPSE2 may be preferred. For example, it is unclear if GLIMPSE2 demonstrates better performance at higher coverage levels, as the benchmarking ends at 2x coverage, where differences seem to favor GLIMPSE2. This omission limits a comprehensive understanding of QUILT2's strengths and weaknesses. Additionally, computational time is reported at 0.25x coverage, where GLIMPSE2 is shown to be about four times faster, but it is unclear if this difference is maintained at higher coverage. A more comprehensive and fair evaluation of methods across coverage levels and contexts would enhance the manuscript's scientific validity.\nLack of citations to previous work\nThe authors seem to overstate the novelty of their approach by not clearly citing relevant prior work underlying QUILT2. Several key computational features of the QUILT2 method have been previously introduced by other methods, yet there is a noticeable lack of acknowledgment of these contributions.\nThe msPBWT shares a lot of similarities with the multi-allelic PBWT (Naseri et al, 2019) and/or Syllable-PBWT (Wang et al, 2023).\nPBWT for haplotype selection has been used in numerous methods already: EAGLE2, SHAPEIT4, SHAPEIT5, IMPUTE5, GLIMPSE1, and GLIMPSE2.\nTreating rare and common variants differently to speed up computation has been used in BEAGLE5, IMPUTE5, SHAPEIT5, and GLIMPSE2.\nAlso, not discussing methodological differences between GLIMPSE2 and QUILT2 is a significant oversight given the similarities between the two methods. The first mention of GLIMPSE2 occurs in the method comparison and does not provide any context for this comparison.\nPresentation of results\nI believe the whole Figures 2 and 3 should be carefully rethought and redesigned for clarity and more importantly for transparency.\n- Figure 2a: The comparison of methods for coverage > 0.25x is simply not visible.\n- Figure 2c: The choice of 0.25x coverage is unclear. It seems that the difference between methods do not really differ to what is already shown in 2a. What does the comparison look like for 1x coverage?\n- Figure 2d: Accuracy is presented differently for this analysis compared to 2a-c. That makes results harder to read.\n- Figure 3a-b: These two figures need to show computing time for datasets comprising thousands of samples and even more, cases where compute time and memory usage are a real issue.\n- Figure 3c: Why #samples=8 here and why stop at 1x while most plots show results up to 2x?\nFinally, it seems unfair to choose the reference panels for which QUILT2 performed the best in some panels of fig. 2: 2a shows results for UKB, while 2d for 1000G. The authors should either show results for all reference panels as the main figure or pick one and stick to it.\nLack of reference for major common use cases\nMajor low-coverage WGS initiatives, like Blended Genome-Exome (BGE), focus on higher sequencing depths (2-4x), yet the manuscript does not address QUILT2’s application at these depths for large-scale cohorts. I think all benchmarks should at least include 4x coverage configuration beside the coverages already presented.\nThe timing benchmark seems also to miss the point: compute time for single-sample imputation is not anymore an issue. Scalability for hundreds of thousands of samples as in typical GWAS is actually the critical aspect. This is not addressed in the benchmark, nor in the discussion. Readers will want to know how running times and RAM usage scale on large datasets comprising many thousands of samples.\nAncient DNA benchmark\nThe authors do not describe how imputation has been run on ancient DNA. Established pipelines involving GLIMPSE are already used in real-world ancient DNA applications, with GLIMPSE2 authors explicitly advising against its simplistic genotype likelihood model for complex scenarios like aDNA. If the authors choose to retain the ancient DNA benchmark, they should compare their tool with the recommended pipeline outlined by Mota et al., Nat Comm. 2023.\nLong reads\nThe application of QUILT to long reads is a key advantage, but this was demonstrated in prior work. It is unclear why similar results are highlighted in a main figure here or what this benchmark aims to achieve. If the focus is on long-read performance, the comparison should emphasize improvements over QUILT1, rather than contrasting with GLIMPSE2, which is not designed for long reads.\nNIPT analysis\nWe find this part of the manuscript useful to the community, nicely carried over and interesting. This is a significant advance in the field.\nMinor comments\n- The imputation performance of QUILT2 at indel sites is not clear. Could the accuracy for SNPs and indels be separated?\n- The statement \"For both human and animal studies, QUILT has been shown to be the most accurate and robust approach when a reference panel is available\" is overly assertive. Both the GLIMPSE2 study and this manuscript suggest minimal performance differences between methods.\n- The claim \"QUILT2 should retain greater accuracy over GLIMPSE2 in various non-human settings, as the effective population size, and hence SNP density, is often much higher than in humans\" requires further clarification. Where is the data showing this?\n- The statement \"We also showed the extent to which large haplotype reference panels can improve both phasing and imputation performance\" does not represent a major contribution of this work. This has been demonstrated previously when the reference panel was introduced (citation missing), and this study largely reaffirms those findings.\n- References for some reference panels (HRC, UKB, etc …) are missing.", "reviewer_4": "The authors introduce QUILT2, a method that builds on the original QUILT1 for accurate genotype imputation using reference panels. QUILT1 has a linear computational complexity that scales with reference panel sizes and number of SNPs, both of which are becoming increasingly large, thereby slowing the method runtime. QUILT2 addresses this with a scalable, memory-efficient design for low-coverage WGS imputation. Notably, it also enables imputation of fetal versus maternal genomes from cell-free DNA in non-invasive prenatal tests—an underexplored application. The paper effectively demonstrates the method's robustness and efficiency, though a few comments are provided below.\nMajor comments:\n1. Regarding framing, NIPT is a useful application in its own right, but the framing that it will substantially improve the ability of GWAS to reach millions of participants is somewhat of a funny one, particularly as most phenotypes studied by GWAS are adult-onset diseases or adult quantitative traits. Maybe reconsider NIPT framing as contributing substantially to GWAS in the introduction? Relatedly in the discussion, it is not clear how having common variants from a fetus would change clinical practice in cases of maternal preeclampsia or gestational hypertension, so perhaps clarifying or removing would be helpful.\n2. Relatedly, is using fetal genomes in a GWAS / PRS feasible at current accuracy rates? Variants at 1% MAF are only a little over 60% accurate even at 4x coverage and 0.3 FF. Could it introduce false positives? Most of the significant hits are only nominally significant both for mother and fetus (Fig 5b). Similarly, with the current state of PRS as predictors for traits, does matching array data performance with cfDNA really change the feasibility of its use in the clinic? This is an interesting and valuable analysis, but some text in the discussion addressing this caveat could be merited. Could include the study “Novel insights into the genetic architecture of pregnancy glycemic traits from 14,744 Chinese maternities” by Z. Zhu et. al. (2024, Cell Genomics) in the discussion relating to using NIPT data in GWAS.\nMinor comments:\n1. The manuscript could benefit from a careful read through to streamline and improve readability as there are some very complex sentences that are hard to parse. Some examples:\na. The first sentence of results lacks clarity: QUILT2 is based upon, and improves, our previous model QUILT1 for diploid genotype imputation, which operates, which uses an iterative approach, as follows. Suggest to remove “which operates” from the sentence.\nb. There are two i.e.’s in this 55-word sentence. It would be more readable if broken up or streamlined: “However, state-of-the-art imputation methods for lc-WGS are designed for diploid samples, i.e. where the sequencing reads come from two haplotypes 3–6, and there is no dedicated method designed for utilizing the triploid nature of NIPT data, i.e. where the sequencing reads come unequally from the maternal transmitted, maternal untransmitted, and paternal transmitted haplotypes.”\nc. “For diploid samples, we used the WGS derived UKB panel to show that QUILT2 is approximately as accurate as QUILT1, but substantially faster, and uses much less memory, and as such, should be much less expensive in cloud based applications.” Example streamlined suggestion: For diploid samples, the WGS-derived UKB panel showed that QUILT2 is as accurate as QUILT1 but significantly faster, more memory-efficient, and cost-effective for cloud-based applications.\nd. “In addition, while QUILT2 is slower than GLIMPSE2 when processing many samples, it is faster, by an order of magnitude on very large reference panels, when processing only a very small number of samples.” Suggestion: While QUILT2 is slower than GLIMPSE2 for many samples, it is an order of magnitude faster when using very large reference panels to process a small number of samples.\ne. Line 99-101: suggest that the sentences “Because of the second part of the iterative process, QUILT1 has linear time complexity in the number of haplotypes and variants in the full reference panel. For RAM and speed reasons…” be rewritten to “Given the two-part iterative process, QUILT1 has linear time complexity dependent on the number of haplotypes and variants in the full reference panel. To minimize RAM and runtime…”\n2. Including the number of SNPs as another factor in the runtime / memory usage could be interesting. For example, subset the SNPs in the lc-WGS data and see how it varies with the number of SNPs included. There is some analysis on this for the QUILT-nipt method vs QUILT2, though it only focuses on all variants vs only common, and does not mention the number of SNPs included at these thresholds.\n3. Given that low-coverage WGS imputation methods often differ from GWAS array imputation methods in that the former relies on genotype probabilities whereas the latter relies on genotype calls, is QUILT2 incompatible with GWAS array data? This point is worth clarifying to ensure use cases are clear to readers.\n4. Common variant definition of MAF=0.1-0.2 seems a little funny. Why not just MAF > 0.1?\n5. Given that QUILT2 imputed more accurately than GLIMPSE2 with UKB WGS but not phase 3 1kGP and HRC, is there better genotype error tolerance in GLIMPSE2? Notably, a 30x WGS version of 1kGP is available.\n6. In Figure 2a, it is hard to tell the coverage lines apart because points overlap. The authors might consider using a color gradient for increasing coverage, with QUILT2 vs GLIMPSE instead indicated by line type (e.g. dashed vs solid). Comparing accuracy across Figure 2a, 2b, and 2c is challenging as the scale is changing.\n7. Very interesting to see the ancient DNA imputation results. Though in Figure 2d (and corresponding supplemental figures) - is there a particular reason to switch from MAF as x-axis (for figures 2a-c) to coverage as x-axis and color as MAF in Fig2d? Could consider switching the axis to match the other plots, though not critical to do this.\n8. Non-reference concordance is a bit more of a forgiving metric than aggregate r^2. Do the authors care to comment on the relevance of this in the comparison of aDNA vs modern samples?\n9. In Figure 3, the legends are a bit confusing, since those shown in A differ from in B but are relevant to A-C. I would suggest extracting a common legend to the bottom of this figure, then reducing the y-axis for panel b, as it is hard to distinguish between lines. The description of memory usage in the text describing this figure could be clearer about the relative comparison between QUILT2 and GLIMPSE2 (note that GLIMPSE2 is not mentioned here): “We note that QUILT2 uses slightly higher RAM given 100 samples to run, and QUILT2 also runs faster with lower RAM for lower coverage data (Figure 3C, Supplementary Figure 5, 6).”\n10. The NIPT comparisons in Figure 4 are nice. The authors might want to comment on whether imputing only common fetal variants is the best strategy given the focus on embryonic screening for structural variants and potentially Mendelian variants, or whether lcWGS is an uncommon current application but with approaches like QUILT2 enabling more.\n11. Any thoughts on why there is a downward trend in aggregated R2 for common SNPs in Fig 4, both maternal and fetal genomes? More pronounced for diploid than NIPT and more pronounced with increasing coverage.\n12. Why were only quantitative traits used for the GWAS / PRS analyses? Do the authors expect a change in behavior for logistic mixed models? Don’t necessarily have to test binary traits, but just include a sentence explaining why quantitative should behave the same as binary traits.\n13. For the NIPT GWAS locus identification part, the authors use a loose definition of p < 0.05 to describe identification. Figure 5B shows that the vast majority are identified with nominal significance, so a concordance check on effect size (same direction at least?) may also be warranted. It might be helpful to reorder these bars so that the mother array followed by mother with decreasing coverage are shown, followed by the simulated fetus, as it’s otherwise hard to compare the relative change in performance. Blue and orange coloring in Figure 5B should be switched to not be confused with the blue/orange for fetus/mother color scheme in other plots. In terms of relative gains in discoverability expected across technologies, this publication might be worth discussing: Gaynor, S.M., Joseph, T., Bai, X., Zou, Y., Boutkov, B., Maxwell, E.K., Delaneau, O., Hofmeister, R.J., Krasheninina, O., Balasubramanian, S., et al. (2024). Yield of genetic association signals from genomes, exomes and imputation in the UK Biobank. Nat. Genet. 56, 2345–2351.\n14. Given the focus on NIPT, how many common clinically actionable HCMG genes or other common variants relevant to fetal Mendelian disease (e.g. CFTR deltaPhe508) are well-imputed? Is low-coverage WGS the best strategy here, or might a WEGS or blended strategy that enriches for higher coverage at interpretable coding variants be worth discussing?\n15. The discussion point about the clinical actionability of de novo CNVs e.g. 22q11.2 is interesting but also confusing, as imputation methods like QUILT2 use a reference panel to only impute existing variants in the population. This is framed as a potential future development opportunity for QUILT2 that seems highly speculative at present as a class of variation undetectable by any imputation method.\n16. Methods: “(cid:1) is the parameters of our model” - this parameter needs to be defined more clearly.\n17. Line 623: “Real big” reference panel - change to “very large” or “one of the largest”\n18. Not sure if it’s the pdf formatting, but for the notation in the supplemental note, there is no semi-colon.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "The manuscript presents QUILT2, a genomic imputation tool, demonstrating enhanced scalability, computational efficiency, and accuracy over existing methods. It highlights QUILT2’s performance across diverse scenarios, including short and long sequencing reads, ancient DNA, and non-invasive prenatal testing (NIPT), where it imputes both maternal and fetal genotypes.\nI really appreciate the quality of the work; it’s very well presented in the paper, and the accompanying online resources (like the Snakemake workflow (!), tutorials, and comprehensive documentation) are excellent. Thanks!\nThank you for taking the time to review our manuscript and for the positive comments.\nMajor Comments:\nThe authors mention two technical advancements (reduced computational complexity and two-stage imputation). How do these changes to the Positional Burrows-Wheeler Transform (-> resulting in msPBWT) influence the trade-off between speed and accuracy compared to PBWT? Also, what exactly contributes to QUILT2’s improved accuracy, considering that GLIMPSE also uses PBWT (to the best of my knowledge, not stated in the manuscript)?\nPBWT has constant computational complexity in haplotype reference panel size, and so does msPBWT, so in that sense, there is no trade off in speed. However, msPBWT is better in that it is much more memory efficient, and should be slightly faster, compared to PBWT, at the slight expense of accuracy.\nThe changes are motivated as msPBWT is a memory-efficient variant of PBWT that is designed to work around one of the key assumptions that facilitates rapid imputation in QUILT and QUILT2, where recombinations are disabled between except between every 32nd SNP. For densely imputed data, where recombinations occur infrequently compared to the number of SNPs, there should be nearly no loss in accuracy in using msPBWT vs PBWT, just gains in RAM and speed.\nMore generally, gains in accuracy in QUILT2 vs GLIMPSE2 are due to the underlying read-aware model in QUILT2 vs per-SNP independent genotype likelihood in GLIMPSE2, as well as potentially how the msPBWT / PBWT results are used by QUILT2 and GLIMPSE2 to select the new set of conditioning haplotypes. Gains in accuracy in QUILT2 (using msPBWT) vs QUILT1 (using a Li and Stephens model) are due to the much larger haplotype reference panels that QUILT2 is able to use.\nWe’ve made it more clearer now in the Introduction and Discussion by presenting the methodological differences between the methods.\nIs QUILT2 intended to fully replace QUILT1 for all reference panels / use cases, or are there cases where QUILT1 remains preferable? Maybe this should be stated in the manuscript.\nQUILT2 is not intended to fully replace QUILT1 as any msPBWT / PBWT based approach will be faster but less robust. As an example, consider random very rare genotyping errors in the haplotype reference panel. Using a Li and Stephens model, these will contribute very little to the copying probabilities when selecting new conditioning haplotypes, so can be ignored. However, PBWT style approaches require exact matches, so any single genotype error will break long exact haplotypes, hampering the ability of the model to use them in the next round of imputation.\nAs such, QUILT1 would be preferred in situations where there are smaller or less accurate haplotype reference panels, for instance in non-human studies. We now made this result more clearer in the text.\nIn the second paragraph of “Diploid imputation performance”, we highlighted:\nFor imputation of modern samples with smaller panels (KGP and HRC), we found that QUILT2 is less robust and accurate than QUILT1 and GLIMPSE2 for rare variants with default settings (Supplementary Figure 1).\nIn the Discussion, we noted:\nFor diploid samples, the WGS-derived UKB panel showed that QUILT2 is as accurate as QUILT1 but significantly faster, more memory-efficient, and cost-effective for cloud-based applications. However, QUILT2 is less robust than QUILT1 for smaller reference panels.\nWhy is the imputation of long reads generally less accurate, as shown in Figure 2? Is this somewhere mentioned?\nWe suspect the decrease in accuracy for the long read analysis shown in Figure 2 is predominantly driven by the higher error rate of the sequenced bases in the reads and the error in long read alignment. Secondarily, at low coverage, longer reads are at a disadvantage in terms of coverage of the two diploid haplotypes in normal imputation. As an extreme example, if the average read length is 100 kbp, then at 0.1X, it would be fairly common to have 0 sequencing reads from one haplotype over a 1 Mbp span.\nWe’ve previously mentioned this in the QUILT1 paper, saying\nThis explains why short-read (Illumina-based) sequencing outperforms long-read data (exemplified by ONT) because the slightly improved read label assignment possible from long reads is outweighed by the much lower per-base error rate of the short-read data.\nWe think it’s probably best not to mention this again in the text in this manuscript but are willing to do so if the reviewers think it wise\nThe authors claim that QUILT2 significantly outperforms GLIMPSE2 in scenarios where a single sample is imputed with a large reference panel. Can they elaborate on what makes QUILT2 so efficient in this context?\nTo be clear, this is in terms of speed not accuracy, as both methods impute samples independently from each other. As for speed, to the best of our understanding of how GLIMPSE2 works, this is caused by how both methods read in and process the reference haplotype data.\nWith QUILT2, haplotype reference panel processing is done beforehand, including making the necessary data structures for msPBWT, and the resulting data structures are stored as simple binary objects. Loading these up afterwards for any subsequent run is very fast and does not contribute meaningfully to run time.\nHowever, the process by which GLIMPSE2 loads and makes available the reference data for later PBWT based analyses is relatively slow and has linear computational complexity in the number of haplotypes in the haplotype reference panel. This fixed time cost contributes only minimally with large sample size, but contributes substantially at small sample sizes.\nHence why at low numbers of samples, QUILT2 can be substantially faster (more efficient) than GLIMPSE2.\nWhen using QUILT2 on the UK Biobank Research Analysis Platform (RAP), are additional pre-processing steps needed for the reference panel? (e.g. QUILT2_prepare_reference.R ) Specifically, does the panel need to be in VCF format, and if so, is this format already available on RAP (I guess not)? Are there post-processing steps for QUILT available (e.g. ligate)? I think users would highly appreciate a tutorial on how to create reference panels (e.g. locally or on UKB) run the imputation process + post-processing steps. How does the chunk size affect the quality? This could be interesting when conducting a WGS imputation for example on a local Slurm cluster.\nThank you for the comments. On the UKB RAP, there are reference VCF files available, which are generated and phased by BEAGLE5 and Shapeit5. The snakemake workflow we developed has the flexibility and the capability for pre-processing the reference panel (e.g. split the genome into chunks) and post-processing the results (ligating the chunks, measuring the accuracy and generating reports if a truth VCF is specified). So while we have not provided an explicit tutorial, the snakemake workflow should help users operate in the UKB RAP environment. Dedicated tutorials can be found elsewhere (e.g. from BEAGLE5). Note that in the benchmarking, we didn’t include the cost of additional pre-processing or post-processing steps if there are any, for any of the methods.\nIn terms of the chunk size, most imputation methods should be robust to the chunk size being too small, as long as the buffer is sufficient. Similarly, since most methods use approaches with conditioning haplotypes, the region should not be too large. In practice, windows of about 2-5 Mb seem appropriate for most modern methods, including QUILT2.\nFor the new QUILT2-nipt approach, is the fetal fraction required as an input parameter? Is this somewhere stated?\nYes, fetal fraction is required as an input parameter, and yes, this is stated on the website under the NIPT section, where it’s listed as a requirement “fflist” which is a file with a one-to-one correspondence with “bamlist” or “cramlist”. We’ve now added an explicit statement in the second paragraph in Methods “which is required as an input and can be obtained from NIPT report 52,53”, which includes references to how to calculate it including the paper “Bioinformatics Approaches for Fetal DNA Fraction Estimation in Noninvasive Prenatal Testing.” We chose to use an external value in QUILT2, rather than re-estimate it, for simplicity, and as it is easy to calculate using existing tools. We also showed QUILT2 is robust to misspecification of this parameter in the text.\nAncient DNA is briefly mentioned in the abstract but lacks focus in the introduction or in the title. Could the authors provide more context earlier in the manuscript to align the emphasis on ancient DNA with its later discussion? Or to state it differently, is the use case of ancient DNA central to the manuscript?\nThank you for the comments. We think showcasing the imputation performance of QUILT2 on ancient DNA advances a central theme of this manuscript that QUILT2 is able to accurately handle many different types of sequencing input and applications. In addition, the imputation of ancient DNA can also benefit from large reference panels, which was shown by the GLIMPSE2 paper and was also demonstrated by our results. We’ve added more text to this motivation in the Introduction\nAdditionally, imputation with reference panels has been demonstrated as a reliable approach for ancient DNA (aDNA) studies 7. Leveraging large reference panels, aDNA imputation can be further improved 6. Furthermore, human genetic studies have begun using long-read sequencing more and more for investigating disease-associated variants, and strategy of mixing high and low coverages is favored for population-scale study 8.\nGiven the reliance on simulated NIPT data, how closely do these simulations mirror real-world scenarios? Are there potential biases that might affect the conclusions?\nThank you for the comments. We chose this kind of simulation using real trio or duo haplotypes because the simulated NIPT data should then have the exact same recombination patterns between the fetus and the mothers genomes that are present in NIPT data. This means that we do not have to make assumptions about recombinations in meiosis, and so in that regard, we’re mirroring very closely the real-world case. However, at the molecular level, there is a difference between the cfDNA in the plasma of pregnant women and saliva DNA of non-pregnant adults (the one we used for simulation). It has been observed that fetal-derived cfDNA fragments tend to be shorter and enriched at different locations than the maternal cfDNA fragments\n(https://pmc.ncbi.nlm.nih.gov/articles/PMC7935715/). In this revision, we analyzed one real whole genome sequenced NIPT sample, which has 20X coverages of plasma cfDNA with estimated 6% FF and 30X coverages saliva DNA (Extended Figure 4), which confirms the effectiveness and improvements of QUILT2-nipt over QUILT2-diploid for imputing the mother.\nIn the first paragraph of “NIPT cfDNA imputation performance”, we added:\nIn addition to simulations, we also confirmed the effectiveness and improvements of QUILT2-nipt over QUILT2-diploid for maternal genotype imputation with one real NIPT sample (Extended Figure 4).\nAre there specific genetic architectures or data types where QUILT2 might be less effective compared to other imputation methods?\nThank you for the comments. Most imputation methods can either only work on a single data type, or are highly optimized for that data type. For instance QUILT2 and GLIMPSE2 for sequence based imputation, while other programs like Beagle (recent versions), Shapeit and IMPUTE are designed for genotyping microarray based imputation. As such, in designing a study, a researcher would first choose the genotyping method (sequence or arrays), and then their preferred imputation strategy.\nCertain study designs would certainly favour low coverage sequence versus arrays. From a cost and genotyping accuracy perspective, low coverage sequence should be favoured over genotyping microarrays for most human studies of common disease, as has been shown in our previously QUILT work, as well as by the GLIMPSE authors. This is driven by the increased accuracy of lc-WGS and imputation versus arrays at current price points, and that the genetic architecture suggests that most causal variants are comparatively common. There are some exceptions in which arrays would be favoured over lc-WGS. In particular, arrays directly type genetic variants. As such, in cases where there is a substantial contribution to the genetic architecture from variants which are hard to impute, in particular rare population specific variants which are poorly captured in large haplotype reference panels, or more generally, variants with little LD to other variants.\nThe results suggest that QUILT2-nipt could significantly enhance GWAS and PRS performance using imputed fetal genotypes. Could the authors elaborate on the implications for clinical applications, particularly the accuracy requirements in prenatal risk assessments? Are there recommended r² thresholds for polygenic score calculations based on fetal fraction? Given the ongoing debate around the clinical utility of PGS and the uncertainty of fetal genotypes (depending on seq depth + FF, see also below), this topic could result in diverse opinions.\nThank you for the comment. For GWAS, we think it is a simpler case, where we envision the data already having been collected. The genotypes would then be imputed with accuracy as expected based on sequencing depth, and the results could be used like a normal GWAS.\nFor clinical applications using PGS, it is more complicated. As shown in the paper, increasing sequencing depth improves fetal genotype imputation accuracy, and as such, the r2 for polygenic score calculation - i.e. the squared correlation between the polygenic score from fetal NIPT imputed genome vs a high quality WGS genome. Therefore it should be a matter of cost, not technical possibility, at least for PRS and common variants. As such, the question would then be one of value - whether the costs and benefits of the higher coverage required to ensure a certain reliability of prenatal PRS for the PRS to be useful, outweigh the cost and benefit of not having that PRS but still being able to perform aneuploidy detection. We note that methods which take the uncertainty of genotype imputation into PRS interpretation are emerging, e.g.\nhttps://pubmed.ncbi.nlm.nih.gov/37490908/.\nFor clinical diagnoses relying on rare variants, lcWGS based imputation is unlikely to be sufficiently accurate, and in such cases, it is more likely that deep fetal exome sequencing would be preferable.\nWe have changed the text in the Discussion\nFrom a clinical perspective, relevant traits that are important for health care, such as hyperglycemia, gestational diabetes and birth weight, are genetically influenced by both maternal and fetal genomes 14,43,46,48–50. Therefore, the imputation of the fetal genome before childbirth affords the opportunity to calculate PRS prenatally, which may impact future health risk management if the PRS is sufficiently predictive, though future studies will be required to fully evaluate this potential. Recent papers, which take the uncertainty of genotype imputation into PRS interpretation into account, are of particular relevance here51. Nonetheless, the imputation of specific variants, in particular rare variants, remains challenging. Therefore it is unlikely to be used for clinical diagnosis, and in that context, recent approaches using non-invasive deep fetal exome sequencing are more promising 22,23. We note that further extension to QUILT2-nipt should still operate on fetal exome sequencing, so in addition to accurate rare variants, such data would also allow for genotype imputation and PRS calculation.\nMinor:\nThe abstract describes the fetal genome imputation results using QUILT2 at 0.25× sequencing coverage with a 10% fetal fraction as “accurate,” yet the r² is only 0.465. This statement seems potentially misleading—could the authors clarify this assessment?\nThank you for this comment. We had meant “accurate” to refer to the mother, but this phrasing is unclear. We’ve changed the language to the following.\n“... we see accurate imputation of the mother (“r2 = 0.966) and modest imputation of the fetus (r2 = 0.465)”\nLine 154: Please confirm the minor allele frequency (MAF) levels used to define common SNPs.\nWe stratify SNPs into bins based on gnomAD allele frequency, and while we show all bins in figures, we focus on the following exemplar bins, as stated in the text. We note that this would give very similar results, and effectively identical qualitative results, had we defined common SNPs otherwise e.g. 5-50%.\nAs exemplars, we used ‘very rare’ to refer to SNPs with MAF of 0.01-0.02%, ‘rare’ to refer to SNPs with MAF of 0.1-0.2%, and ‘common’ to refer to SNPs with MAF of 10-20%.\nThere are several typos in the manuscript that need to be fixed.\nWe apologize for any typos. We have done our best to correct any typos we have found, but please let us know if you find more.\nReviewer #1 (Remarks on code availability):\nSee above. I think the quality of the Github repo is very very high.\nThank you!\nReviewer #2 (Remarks to the Author):\nSummary of the key results, Originality & significance\nLi et al. reported an improved method for genotype imputation in data produce form short- and long-read sequencing data, named QUILT2. The methodology is primarily based on their previous work (QUILT) that makes use of sequencing reference data to impute undetermined genotypes of shallow-sequencing data (PMID: 34083788). The authors stated two technical and a methodological improvement over their previously published works (PMIDs: 27376236, 34083788). Technical improvements include: (i) a memory efficient PBWT with the purpose of finding haplotypes in the reference panel that shares long matches with constant computational complexity and “very low” memory footprint; (ii) a two-stage imputation approach that first samples read labels and makes use of common SNPs to find an optimal subset of haplotype reference panel. Methodological improvement includes assumes that three haplotypes are present in the NIPT data.\nThis reviewer did not find these improvements conceptually novel relevant to their previous works (PMIDs: 27376236, 34083788) and the existing work GMPSE2 (PMID: 37386250). In additions, although the authors claim applicability of their method (three-haplotype option) for non-invasive prenatal testing (NIPT) by discerning fetal and maternal genomes and imputing the fetal genome, in my view, the assumptions are far from real cfDNA component in blood stream of pregnant women (see below). And that it has been performed on simulated data from whole-genome sequencing data of trios. For the latter, the r2 values are specifically too low to claim that the methodology works on low-pass sequenced NIPT data.\nThank you for taking the time to review our manuscript and for the constructive comments.\nMajor comments:\n- The QUILT2-nipt method assumes the observed sequencing reads are distributed among three haplotypes (maternal transmitted, maternal untransmitted, and paternal transmitted). However, the complexity of real-world data, including recombination in maternal haplotypes and errors in fetal fraction estimation, may introduce inaccuracies that are not addressed. How do you address recombination sites between the maternal homologs?\nWe’ll answer this jointly with the comment below. Also, we analyzed one real NIPT sample during this revision.\n- Given my comment above, it is incorrect to assume that NIP data has “triploid nature”, it is rather tetraploid (prental chromosomes recombine).\nYes this is an issue theoretically and it is one that might negatively impact imputation, but probably less than you think. To expand, you are correct that in the blood of the pregnant mother, if we consider the entire chromosome of origin of any fragment of DNA, than these chromosomes will have four haplotypes: the first haplotype of the mother, the second haplotype of the mother, a mosaic maternal haplotype (derived from crossovers of the maternal haplotypes) and a paternal haplotype. However, if we consider a single arbitrary short stretch on the chromosome, the mosaic maternal haplotype will exactly match one of the two maternal haplotypes, and as such there are only three distinct haplotypes, present at three frequencies. The problem therefore is the scale at which these two truths diverge, and what it means for imputation.\nHuman recombination is about 1 cM / Mbp, with one obligate crossover during meiosis. As such we would expect on the order of a few recombinations (crossovers) per chromosome, with a higher rate in the telomeres. With QUILT2, we typically impute in windows that are between about 1 to 5 Mbp. So as an example, on a 100 Mbp chromosome, if there were two maternal crossovers, and if we imputed in 2 Mbp windows with 500 Kbp buffers, we would expect this problem in the central region in about 2 out of 50 windows - in the other 48 windows, there is no error in this key assumption in QUILT2-nipt. In the 2 problematic windows, in the worst case scenario, that the crossover happens in the middle, then we’ll incorrectly impute one of the halfs. As such, for the imputed maternal chromosome, we would expect imputation problems in about 2 x 1 Mbp / 100 Mbp or about 2% of the chromosome. We note that this problem would mostly affect the imputed fetal genome (through difficulties in determining the transmitted maternal haplotypes), and would only very minorly affect the imputation of the maternal genome (with both maternal haplotypes).\nThis assumption could be fixed through an alternate modelling strategy that adds in a new hidden variable that explicitly models which of the two maternal haplotypes is the transmitted one, along the region to be imputed. This is a consideration for a future model, however this would be a substantial change to the model, and the potential benefit of this (improving fetal imputation in ~2-5% of the genome, or less with shorter windows) would need to be balanced against other potential changes.\nIn terms of errors in fetal fraction estimation, we showed in the manuscript in Figure 5a that the model is robust to errors in this parameter by giving QUILT2 biased fetal fraction. These errors are in line with, or above, what would be expected from commonly used methods to estimate fetal fraction.\n- R2 is very low for real NIPT data. Therefore application of QUILT2 can lead to wrong SNP calls. Could you clarify which parts of genomes cause this low r2 value?\nTo clarify, we impute/genotype known SNPs in a reference panel, and for many people, the preferred output is a dosage, which is the estimated number of alternate alleles. The r2 is very accurate for the mother, and it’s also accurate for the fetus with effective sequencing coverage > 0.4x, e.g. r2 >0.866 for 4x cfDNA with 10% Fetal Fraction, for common variants (Also, see our response to Reviewer 4 Minor Comment 14 about imputing common variants that causes Mendelian diseases ).\nRegarding the imputation error, we expect mostly a generally uniform imputation error rate across the genome, however some subtle differences are possible. In particular we expect that non-uniform imputation error across the genome is principally caused by variation in recombination rate, variation in read mapping quality, and variation in haplotype reference panel accuracy (including phasing accuracy), itself also likely related to the previous two points.\nFor the recombination rate, a higher recombination rate will affect imputation both for the reason given in the previous answer, and as a higher recombination rate will shorten the length of haplotypes copied identically between the reference panel and the focal haplotype to be imputed, making its detection using lc-WGS more difficult. As such in particular the telomeres should have lower accuracy, as should regions of the genome with high or unusual recombination rate and coalescent histories, like the MHC.\nFor read mapping quality, clearly erroneous read mapping will negatively impact imputation. However, with appropriate minimum mapping quality thresholds, this should be a minor concern. The main concern would likely be from reads from regions of the genome with recent duplications, where one of those duplications is not captured appropriately in the reference genome.\nFinally variation in reference panel accuracy, both singleton genotyping errors, and phase switch errors, will affect imputation accuracy, and will be non-uniformly distributed. We expect this to largely follow the same regions as the previous two issues.\n- The method has higher accuracy for common variants but its applications with (very) rare variants is challenging and substandard. While this is expected, the reliance on filtering out rare variants early in the iterative imputation process could overlook clinically relevant but less frequent variations. Specially in the context of NIPD. Please clarify how you would deal with this challenge?\nThank you for the comments. Two things here. First, with respect to the iterative process, while the rare variants are not used in the early rounds of imputation, they are used and output at the end. So an output imputation file will contain imputed genotypes at all SNPs in the haplotype reference panel, regardless of frequency.\nSecond, with regards to clinical applications of rare variants imputed from lc-WGS. Our results clearly suggest that we can not impute rare variants to sufficient accuracy for clinical purposes. If that is the aim, then alternative approaches, such as the deep prenatal whole exome sequencing technology recently developed by Brand, H. et al. (High-Resolution and Noninvasive Fetal Exome Screening. N Engl J Med 2023), should be considered. We note that further extension to QUILT2-nipt can work on such data (cfDNA NIPT whole exome sequencing), so in such a case, one could get rare variants directly from the WES, and common variants / PRS across the whole genome using QUILT. We’ve added to the discussion with these points as follows.\nNonetheless, the imputation of specific variants, in particular rare variants, remains challenging. Therefore it is unlikely to be used for clinical diagnosis, and in that context, recent approaches using non-invasive deep fetal exome sequencing are more promising 22,23. We note that further extension to QUILT2-nipt should still operate on fetal exome sequencing, so in addition to accurate rare variants, such data would also allow for genotype imputation and PRS calculation.\n- The authors discuss the potential clinical utility of imputed fetal genomes for traits such as gestational hypertension and preeclampsia but does not provide sufficient empirical evidence or validation studies to support these claims. Without these and my concerns above, this reviewer finds such phrases more over-statements.\nThank you for the comments. We agree that this is an over-statement for gestational hypertension and preeclampsia, where the contribution of the fetal genome is not supported by the paper (PMID: 36947680) we’re referencing. However, there are many papers showing how both the mother and fetal genomes contribute to the birth weights and fetus growth, which are also linked to maternal hyperglycemia and gestational diabetes (ref: PMID: 29463506, 31043758, 32841251, 31043758, 37798380, 38297123). Specifically, the findings by Chen et al. 2020 (https://pubmed.ncbi.nlm.nih.gov/32841251/ ) concluded that, “Alleles elevating blood pressure are associated with shorter gestational duration through a maternal effect and are associated with reduced fetal growth through a fetal genetic effect. Alleles that increase blood glucose in the mother are associated with increased birth weight, whereas risk alleles for type 2 diabetes in the fetus are associated with reduced birth weight.”\nTherefore, we discuss potentials for these traits instead and note that future studies will be required to fully evaluate these potentials with the following changes in text.\nIn the second paragraph of Discussion, the original sentence was\nFrom a clinical perspective, as we move toward leveraging PRS in health risk management, the imputation of the fetal genome before childbirth should afford the opportunity to change clinical practice for relevant traits, such as gestational hypertension and preeclampsia.\nThe new sentence is\nFrom a clinical perspective, relevant traits that are important for health care, such as hyperglycemia, gestational diabetes and birth weight, are genetically influenced by both maternal and fetal genomes 14,43,46,48–50. Therefore, the imputation of the fetal genome before childbirth affords the opportunity to calculate PRS prenatally, which may impact future health risk management if the PRS is sufficiently predictive, though future studies will be required to fully evaluate this potential.\n- The authors briefly state the potential of QUILT2 to identify structural variations and copy number variants but did not provide how or a benchmarking for these applications. Please clarify.\nAs is standard, we have included in the Discussion areas of potential future methodological improvement. This capability does not yet exist in QUILT2, and it is therefore not benchmarked. We chose to highlight this as the methodological changes would be difficult but reasonable, and the potential clinical impact high. Specifically, there are already clinical tests that offer sub-chromosomal CNV calls, e.g. Roche Harmony for 22q11.2 microdeletions, however with limited sensitivity at reasonable specificity. We think it is worth investigating whether an approach based on QUILT2 (which partitions reads) would be more accurate.\nFor further discussion on this, including discussion of our ongoing work in this area, please see our response to Reviewer 4, Minor Comment 15.\n- For the diploid imputation performance, the authors indicate that they “used Chr 20 for convinience”. Why did you make such choice? Is this then correct that all the data presented in this manuscript are using Chr 20 only?\nIt is generally standard in imputation benchmarking to use a single chromosome for extensive testing. This makes analysis faster when there are multiple methods and parameters to evaluate. It is also reasonable as while imputation accuracy varies along the genome (see above) it generally does not vary appreciably between chromosomes. Similar studies like the GLIMPSE/GLIMPSE2 paper do similar things when comparing methods and parameter choices.\nHowever, we do evaluate the performance for the whole genome in some analysis, e.g. we imputed the whole genome to conduct the GWAS/PRS. In addition, in this revision, we analyzed a real NIPT sample and evaluated the performance of QUIILT2-nipt vs QUILT2-diploid by imputing the whole genome as well.\n- The authors devised a 32-bit grid strategy and allowed transion after 32nd pair of SNPs. It is not clear to this reviewer why 32 and how the transitions have happened?\nThis technical detail is the same as for the previously published QUILT (QUILT1) manuscript. That being said, the 32 number is chosen for computational reasons. We restrict ourselves to bi-allelic SNPs so haplotypes can only have reference or alternate bases i.e. 0’s (reference) or 1’s (alternate), and as such, 32 consecutive SNPs in one reference haplotype can be encoded with a single 32 bit integer, which is highly efficient. Transitions (recombinations) are then forced to happen only every 32 SNPs. This is a deficiency, as recombinations might occur between those 32 SNPs, however, in the context of WGS reference panels, where SNPs are overwhelmingly rare, and recombinations infrequent, occurring many thousands or tens of thousands of SNPs apart, it has at most very modest effects on accuracy, while substantially speeding up imputation.\n- Description of QUILT2 for (ONT) long-read sequencing sounds preliminary. Please clarfiy how you make use of long-read sequencing data and kb-Mb size haplotypes. How it differs from shor-read sequencing? How it performs for structural rearrangments? How do you potentially detect errornous SNVs of the ONT data?\nTo reiterate first, we do not discover novel variants including structural variants. The technical details of how QUILT2 can handle long reads are broadly the same as in the QUILT1 paper, and as such, have been covered before. Here we mostly use the long read analysis to reaffirm this capability with QUILT2.\nIn addition, as is standard, we discussed how it can be useful for future method developments. For example, most SV methods are alignment based, and suffer as the true underlying phase of the reads is unknown - it would be easier to genotype them if we could phase them into their underlying haplotypes. In this field, people normally use the “whatshap haplotag” program to accomplish it, however this requires high coverage sequencing data. This could be done using QUILT2, with lower coverage data. More generally we could envision an iterative process of underlying haplotype proposal, alignment, and read phasing, using QUILT2 and other methods.\nFor further discussion on this, including discussion of our ongoing work in this area, please see our response to Reviewer 4, Minor Comment 15\nMinor comments:\n- In the beginning of the Results section “…QUILT2 is based upon, and improves, our previous model QUILT1 for diploid genotype imputation, which operates, which uses an iterative approach…”. Either use “which operates” or “which uses”\nThank you. We’ve corrected the language, i.e. “which uses” was removed.\n- Please make figures that illustrate QUILT2, clearer. Specially for the grid search and using common SNPs as anchors.\nWe appreciate your comment, however there are a lot of key ideas to get across in QUILT2, namely the two technical improvements (msPBWT & using common SNPs as anchors), as well as the methodological improvement (NIPT). It is difficult to get all three across in a single figure without being overwhelming. We chose to focus Figure 1 on the msPBWT and NIPT ideas, and the link between NIPT and diploid. We have two Extended online figures that cover msPBWT (Extended Online Figure 1) and the anchor SNP idea (Extended Online Figure 2) in high detail (these took a lot of time to make!). We hope they are sufficient for the average reader. We note that the caption for Figure 1 includes the following, to direct the reader to these figures\nMore details about the msPBWT algorithm are shown in Extended Online Figure 1. Details about the iterative approach using common and all SNPs, which is not shown in this schematic, is shown in Extended Online Figure 2.\n- The authors state that “In some countries, NIPT is routinely offered to all pregnant women, for instance in the Netherlands, with an average uptake rate of 46%”. This is not up to date. This uptake rate is from the TRIDENT study that was finalized in March 2023. Since April NIPT in the Netherland is a screening program with very high uptake rate. Please correct.\nThank you for the comments. We could not find a more up to date reference for this figure, in either English or in Dutch. We have additionally added a reference to a number for Belgium. The new sentence is\nIn some countries, NIPT is routinely offered to all pregnant women, with recent uptake rates at around half the population or more, with an uptake rate of 79% and 46% for Belgium (79%) and Netherlands (46%) respectively 17.\nWe would happily update this further if provided alternative references to use.\n- Please insert page and line numbers to make the review process easier and more efficient.\nThank you for this nice suggestion. We’ve now added the page and line numbers.\nReviewer #2 (Remarks on code availability):\nThere are different GitHub repositories refered in \"Code availability\" section of the manuscript. Please make a clear scheme how they are related to each other.\nWe deleted some redundant repositories, and reiterated the functionality of each repository. QUILT2 is the main focus of the manuscript and is listed first. The second repository is the code used to run the benchmarking in the manuscript.\nThe codes require extensive documentations, indicating what function or piece of code does what. For instance, \"## known ped, weird format\" is not really informative.\nQUILT2 has a README with instruction and usage information and examples of how to use it. QUILT2 has a defined user-interface after installation which is heavily documented. We do not expect the users to use internal QUILT2 functions so we have not put the same care into documenting them.\nReviewer #3 (Remarks to the Author):\nThe article presents QUILT2, a method for genotype imputation from low-coverage whole-genome sequencing (lc-WGS), building on the QUILT1 framework. QUILT2 integrates a memory-efficient implementation of the positional Burrows-Wheeler transform (msPBWT) and a two-stage imputation process, aiming to improve the speed and scalability of imputation. Additionally, it introduces a novel approach for imputing both maternal and fetal genomes from non-invasive prenatal testing (NIPT) data. While the method shows promise in some applications, especially for NIPT data, the manuscript raises several concerns.\nThank you for taking the time to review our manuscript and for the many constructive comments.\nBiased benchmark\nOverall, the benchmark puts emphasis on QUILT2 outperforming GLIMPSE2, while it seems to avoid showing scenarios where GLIMPSE2 may be preferred. For example, it is unclear if GLIMPSE2 demonstrates better performance at higher coverage levels, as the benchmarking ends at 2x coverage, where differences seem to favor GLIMPSE2. This omission limits a comprehensive understanding of QUILT2's strengths and weaknesses. Additionally, computational time is reported at 0.25x coverage, where GLIMPSE2 is shown to be about four times faster, but it is unclear if this difference is maintained at higher coverage. A more comprehensive and fair evaluation of methods across coverage levels and contexts would enhance the manuscript's scientific validity.\nThank you for the comments. We appreciate this and have reworked the benchmarks and their presentation to try and be more fair to both methods. We have added 4.0x coverage into the benchmarking.\nIn addition, in the text, we have acknowledged the better performance of GLIMPSE2 including for rare variants at higher coverage. The original sentence was:\nQUILT2 was more accurate than GLIMPSE2 at sequencing coverage ≤0.5×, regardless of the reference panel, and particularly at 0.1×-0.25×.\nThe new sentence is:\nQUILT2 was more accurate than GLIMPSE2 at sequencing coverage ≤0.5× regardless of the reference panel, while GLIMPSE2 was more accurate than QUILT2 for rare variants at coverage ≥2.0×\nWe have also tried to better emphasize differences in reference panels. The original sentence was:\nNo method was uniformly superior for these rare SNPs at higher coverage samples, with QUILT2 being slightly more accurate than GLIMPSE2 for the larger panel (UKB) but not for the smaller panels (KGP and HRC)\nThe new sentence is:\nFor imputation of modern samples with smaller panels (KGP and HRC), we found that QUILT2 is less robust and accurate than QUILT1 and GLIMPSE2 for rare variants with default settings (Supplementary Figure 1)\nIn addition, regarding the computational time, we apologize that we forgot to write the coverage of samples used in Figure 3A, which is 1.0x instead of 0.25x as the reviewer mentioned. We now updated Figure 3 to make this more clear, as well as including the cost on RAP for coverage 2.0x as the reviewer requested in another comment. About the broader point, about computational run time as a function of sequencing depth, this is now covered through Figure 3B, where we show cost, driven by runtime, as a function of sequencing depth. QUILT2 is faster for lower coverage samples than higher coverage samples, as some expensive calculations in the HMM can be skipped if windows of 32 SNPs have no reads in them, and more generally, as there are fewer reads to iterate over in the HMM.\nLack of citations to previous work\nThe authors seem to overstate the novelty of their approach by not clearly citing relevant prior work underlying QUILT2. Several key computational features of the QUILT2 method have been previously introduced by other methods, yet there is a noticeable lack of acknowledgment of these contributions.\nThe msPBWT shares a lot of similarities with the multi-allelic PBWT (Naseri et al, 2019) and/or Syllable-PBWT (Wang et al, 2023).\nPBWT for haplotype selection has been used in numerous methods already: EAGLE2, SHAPEIT4, SHAPEIT5, IMPUTE5, GLIMPSE1, and GLIMPSE2.\nTreating rare and common variants differently to speed up computation has been used in BEAGLE5, IMPUTE5, SHAPEIT5, and GLIMPSE2.\nAlso, not discussing methodological differences between GLIMPSE2 and QUILT2 is a significant oversight given the similarities between the two methods. The first mention of GLIMPSE2 occurs in the method comparison and does not provide any context for this comparison.\nThank you for the comment. You are right that we have not done a sufficient job referencing the literature and previous methods, for which we apologize.\nWe note that Syllable-PBWT, b-PBWT (original “binary” PBWT from Richard Durbin) and IMPUTE4 were referenced in the section of “Haplotype matching with msPBWT” in the Supplementary Note, however the reference page was somehow not included when we compiled the Supplementary Note previously, a problem which has now been fixed.\nWe now acknowledge the relevant works. To begin, we have substantially re-worked the second paragraph in the introduction, to offer a brief introduction to QUILT as well as GLIMPSE, as well as what motivated GLIMPSE2.\nPreviously, we developed QUILT for rapid genotype imputation of lc-WGS data with a reference panel 4. Around the same time, Rubinacci et al. introduced GLIMPSE, another method for genotype imputation 4,5 . Both QUILT and GLIMPSE share a core two-stage imputation framework, where per-sample imputation leverages a small set of conditioning reference haplotypes. For haplotype matching, QUILT employs the Li and Stephens model 9, which is more robust, whereas GLIMPSE relies on the positional Burrows-Wheeler transform 10 (PBWT), which offers computational efficiency. QUILT and GLIMPSE further differ in how they treat input data, with GLIMPSE using per-SNP genotype likelihoods, while QUILT uses and phases the sequencing reads directly, which is preferable for linked reads or long-read data. Overall, QUILT is at least as accurate as GLIMPSE and more robust, however GLIMPSE tends to be faster 5. Since the publication of these methods, reference panels have expanded significantly in both the number of haplotypes as well as variant density, for example many studies now use the UK Biobank whole-genome sequencing (UKB WGS) 11. To accommodate these growing panel sizes, the authors of GLIMPSE developed GLIMPSE2 6. GLIMPSE2 focuses on improving speed and facilitating biobank sized reference panels by incorporating more sparse data structures for both speed and memory efficiency, as well as the introduction of a dedicated version of PBWT, termed sparse PBWT.\nIn addition we better reference the history of the key ideas of QUILT2 in the results, first with msPBWT\nThe use of PBWT is common for haplotype identification in genotype phasing and imputation 5,6,29–31, including with GLIMPSE and GLIMPSE2, and derivatives exist that work on non-binary symbols, such as multi-allelic and Syllable PBWT 32,33.\nSecond, for the use of rare and common variants\nAs such, several methods for genotype phasing and imputation differentiate between common and rare SNPs, for both speed and RAM reasons, including GLIMPSE2 and other methods 29,31,34.\nPresentation of results\nI believe the whole Figures 2 and 3 should be carefully rethought and redesigned for clarity and more importantly for transparency.\nWe appreciate the constructive comment. We have re-designed Figures 2 and 3.\n- Figure 2a: The comparison of methods for coverage > 0.25x is simply not visible.\nWe have adjusted the y-axis so that the difference is more visible. We have also added in the 4x benchmarking.\n- Figure 2c: The choice of 0.25x coverage is unclear. It seems that the difference between methods do not really differ to what is already shown in 2a. What does the comparison look like for 1x coverage?\n0.25x was chosen because the data was only collected at slightly above 0.25x coverage. However, since the haplotagged data is not in common use, and to facilitate overall greater readability of this figure, we have removed this analysis from Figure 2.\n- Figure 2d: Accuracy is presented differently for this analysis compared to 2a-c. That makes results harder to read.\nWe made that choice because we followed the design of GLIMPSE2 paper, i.e use non-reference concordance as y-axis and use sequencing coverage for x-axis. We’ve now made the x-axis the same as the other figures, but we stick to the non reference concordance rate for the y-axis because it’s more informative than the r2 of genotype dosages for ancient DNA analyses.\n- Figure 3a-b: These two figures need to show computing time for datasets comprising thousands of samples and even more, cases where compute time and memory usage are a real issue.\nWe’ll answer this jointly with the same comment below in “Lack of major common use cases”.\n- Figure 3c: Why #samples=8 here and why stop at 1x while most plots show results up to 2x?\nWe now added 2x. We will give a more complete answer to why we compare N = 8 for QUILT2 vs N = 100 for GLIMPSE2 for cost reasons in the section below “Lack of reference for major common use cases”.\nFinally, it seems unfair to choose the reference panels for which QUILT2 performed the best in some panels of fig. 2: 2a shows results for UKB, while 2d for 1000G. The authors should either show results for all reference panels as the main figure or pick one and stick to it.\nThank you for the comment. We now focus primarily on the UKB panel.\nLack of reference for major common use cases\nMajor low-coverage WGS initiatives, like Blended Genome-Exome (BGE), focus on higher sequencing depths (2-4x), yet the manuscript does not address QUILT2’s application at these depths for large-scale cohorts. I think all benchmarks should at least include 4x coverage configuration beside the coverages already presented.\nThank you for the comment. We’ve now added 4x in the benchmarking and re-designed the Figure 2.\nThe timing benchmark seems also to miss the point: compute time for single-sample imputation is not anymore an issue. Scalability for hundreds of thousands of samples as in typical GWAS is actually the critical aspect. This is not addressed in the benchmark, nor in the discussion. Readers will want to know how running times and RAM usage scale on large datasets comprising many thousands of samples.\nWe appreciate this suggestion, and we partly agree. However a key consideration of current analyses is that large reference panels increasingly live on the cloud (e.g. the UKB WGS), and can’t be run locally. The cheapest way to run jobs on the cloud is to use low priority, preemptible jobs. On the UKB RAP, high priority jobs are about three times more expensive. As such it makes sense to run jobs with fewer samples so they run quickly. In the GLIMPSE2 paper, in their supplement, we note\nIn practice, we found that jobs that run on large machines (such as “mem3_ssd1_v2_x32”) are interrupted and executed “on-demand” after ~35 minutes of computation.\nWe saw similar results recently as shown in the below screenshot\nThis screenshot shows that 3 out of 11 jobs of GLIMPSE2 were interrupted when run with N=100. For those preempted jobs, we had to rerun several times to be able to finish them normally, with a corresponding increase in cost. Note that all costs reported in the paper are idealized, i.e. cost of the least expensive machine times the amount of time used, as done in the GLIMPSE2 paper.\nBoth QUILT2 and GLIMPSE2 do not use information from other samples when imputing, and as such, can be completely parallelized per-sample. In the supplement to the GLIMPSE2 paper, they show they achieve optimal per-sample efficiency with N ~= 100 samples, so it would make sense to target that sample size in any given job. QUILT2 by contrast has no such per-sample efficiency costs so would presumably be run in very small jobs e.g. 8 samples like we showed before (a multiple of the number of available 4 cores) as we’ve previously highlighted, as downstream tasks like merging VCFs are very inexpensive when done correctly. As such, the issue of running large datasets with many thousands of samples is really just an issue of how to most efficiently run per-sample, which is why we’ve focused as we have on comparatively small sample sizes like shown in Figure 3.\nIdeally, if the cloud platform is not overloaded and all jobs with 1h runtime can be expected to finish in normal priority, it makes sense to run QUILT2 with N =32 samples. Here are the costs for different numbers of samples as an input batch (as always, the idealized cost, assuming there is no pre-emption).\nQUILT2(N=8) GLIMPSE2(N=8) QUILT2(N=32) GLIMPSE2(N=100)\n0.1x 0.2603240 0.3202330 0.2373724 0.09143739\n0.5x 0.3284220 0.3111820 0.2994665 0.08885302\n1.0x 0.4072600 0.3138550 0.3713537 0.08961626\n2.0x 0.5210532 0.3321573 0.4751142 0.09484219\nAs such, we’ve now amended the Figure and the text in the benchmark to reflect these numbers (N = 32 for QUILT2, N = 100 for GLIMPSE2), as these are reasonable sample sizes that someone might choose for arbitrarily large GWAS.\nSince both QUILT2 and GLIMPSE2 do not borrow information across samples, and as the cheapest way to run jobs on the RAP is using low or normal priority jobs which can be interrupted 6, it is advantageous to impute large datasets using small batches, and then to combine results afterwards.\nIn their previous work, the authors of GLIMPSE2 showed that optimal per-sample runtimes were achieved at approximately 100 samples 6. Comparing per-sample costs for QUILT2 (using N = 32) and GLIMPSE2 (using N = 100) we see both methods are inexpensive though QUILT2 (£0.237, £0.299, £0.371 and £0.475 for 0.1×, 0.5×, 1.0× and 2.0× respectively) is more expensive than GLIMPSE2 (approximately £0.09 for all coverages). We note these are idealized estimates, and real costs will depend on prices at run time, server congestion, etc, and furthermore, that default parameter settings for QUILT2 favour accuracy over speed.\nWe further make the following mention in the Discussion\nIn terms of speed, QUILT2 is slower and less efficient than GLIMPSE2 when run at scale, however it is faster when using very large reference panels to process a small number of samples.\nAncient DNA benchmark\nThe authors do not describe how imputation has been run on ancient DNA. Established pipelines involving GLIMPSE are already used in real-world ancient DNA applications, with GLIMPSE2 authors explicitly advising against its simplistic genotype likelihood model for complex scenarios like aDNA. If the authors choose to retain the ancient DNA benchmark, they should compare their tool with the recommended pipeline outlined by Mota et al., Nat Comm. 2023.\nThank you for the comment. In the initial submission we did use the recommended pipeline from Mota et al Nat Comm 2023. We’ve made this more clear by mentioning this in the Results section where we say\nWe also assessed the imputation performance on aDNA from Afanasievo culture (~4.6kya) with three reference panels. We followed the Mapache 38 pipeline used in the GLIMPSE2 paper, which implements the practice by Mota et al. 2023 7.\nWe use the Mota et al pipeline to generate the BAM files used by QUILT2 and GLIMPSE2 directly. We note in the Mota et al 2023 paper that they show that corrected methods for genotype likelihood calculation are only marginally more accurate than using a simplistic model anyway.\nWe also lightly re-wrote the methods. The original text said\nThen we used the mapache 23 snakemake workflow, to map the FASTQ files against the human reference genome GRCh38 to generate the BAM files for imputation benchmarking\nWhile the new sentence reads\nThen we used the mapache 23 snakemake workflow, which implements the best practice for aDNA imputation 7, to map the FASTQ files against the human reference genome GRCh38 to generate the BAM files for imputation benchmarking.\nLong reads\nThe application of QUILT to long reads is a key advantage, but this was demonstrated in prior work. It is unclear why similar results are highlighted in a main figure here or what this benchmark aims to achieve. If the focus is on long-read performance, the comparison should emphasize improvements over QUILT1, rather than contrasting with GLIMPSE2, which is not designed for long reads.\nThank you for the comment. Indeed there are two factors at play: QUILT1 (small panel) vs QUILT2 (large panel); and QUILT2 (long read favourable method) to GLIMPSE2 (less favourable to long reads). It feels worthwhile to evaluate both changes simultaneously, which we now do, in Figure 2C. While this generally highlights the performance of QUILT2 using UKB, QUILT1-HRC is more accurate than QUILT2-UKB at 0.5x coverage. Nonetheless as far more SNPs are imputed using the UKB panel, this would undoubtedly still prove a draw for most users.\nThe re-worked Results paragraph is as follows\nLastly, we evaluated the imputation accuracy using long sequencing reads from Oxford Nanopore Technologies (ONT) by calculating both the r2 at variant level and F1-score at sample level. To investigate the gain of using the UKB panel for long-read imputation, we included QUILT1 with HRC for comparison. We found that QUILT1-HRC performed better than QUILT2-UKB for coverage <1.0x but with fewer variants imputed (Figure 2c). Additionally, we evaluated QUILT2 and GLIMPSE2, the inexpensive and scalable methods that can use the UKB panel. Across all coverages and reference panels, QUILT2 was substantially more accurate compared to GLIMPSE2 for long reads (Figure 2c, Supplementary Table 3), e.g common r2=0.937 vs r2=0.695 at 1× coverage.\nNIPT analysis\nWe find this part of the manuscript useful to the community, nicely carried over and interesting. This is a significant advance in the field.\nThank you for this positive comment. We note that as part of this review process, we have also added in the analysis of a real NIPT sample in Extended Figure 4.\nMinor comments\n- The imputation performance of QUILT2 at indel sites is not clear. Could the accuracy for SNPs and indels be separated?\nThank you for the comments. QUILT2/QUILT1 don’t support indel imputation. We have not benchmarked the performance for indels.\n- The statement \"For both human and animal studies, QUILT has been shown to be the most accurate and robust approach when a reference panel is available\" is overly assertive. Both the GLIMPSE2 study and this manuscript suggest minimal performance differences between methods.\nThank you for the comment. We have toned down this statement, replacing it with\nOverall, QUILT is at least as accurate as GLIMPSE and more robust, however GLIMPSE tends to be faster.\n- The claim \"QUILT2 should retain greater accuracy over GLIMPSE2 in various non-human settings, as the effective population size, and hence SNP density, is often much higher than in humans\" requires further clarification. Where is the data showing this?\nWe have changed this to “QUILT2 should be more accurate than GLIMPSE2 in instances where SNP heterozygosity is high, as would be expected in species with large effective population sizes”. If the reviewer would like, we could further amend or remove it.\nWe do not have data specifically to verify this, however believe it is a reasonable extrapolation. We showed in the paper that QUILT2 retains increased accuracy over GLIMPSE2 for long read sequencing in humans. The difference between using long and short reads is that with long reads, the expected number of SNPs per read is increased. This is akin to using short reads with much higher SNP density. For instance in Anopheles, heterozygosity is ~1.5%, vs 0.1% in humans (from the paper Genetic diversity of the African malaria 1 vector Anopheles gambiae).\n- The statement \"We also showed the extent to which large haplotype reference panels can improve both phasing and imputation performance\" does not represent a major contribution of this work. This has been demonstrated previously when the reference panel was introduced (citation missing), and this study largely reaffirms those findings.\nWe had meant this in the context of the paper. While others have shown this, it was important to show in our work, as, for instance, the accuracy of the NIPT imputation, which is more difficult due to the lower coverage and three haplotypes, depends on it. Nonetheless, we have removed this sentence.\n- References for some reference panels (HRC, UKB, etc …) are missing.\nThank you for the comments. We’ve added these references in the main text at their first mention in the results section on diploid genotype imputation.\nReviewer #3 (Remarks on code availability):\nThe code is both accessible and functional. Additionally, the documentation is comprehensive and thorough.\nThank you.\nReviewer #4 (Remarks to the Author):\nThe authors introduce QUILT2, a method that builds on the original QUILT1 for accurate genotype imputation using reference panels. QUILT1 has a linear computational complexity that scales with reference panel sizes and number of SNPs, both of which are becoming increasingly large, thereby slowing the method runtime. QUILT2 addresses this with a scalable, memory-efficient design for low-coverage WGS imputation. Notably, it also enables imputation of fetal versus maternal genomes from cell-free DNA in non-invasive prenatal tests—an underexplored application. The paper effectively demonstrates the method's robustness and efficiency, though a few comments are provided below.\nThank you for taking the time to review our manuscript and for the positive comments.\nMajor comments:\n1. Regarding framing, NIPT is a useful application in its own right, but the framing that it will substantially improve the ability of GWAS to reach millions of participants is somewhat of a funny one, particularly as most phenotypes studied by GWAS are adult-onset diseases or adult quantitative traits. Maybe reconsider NIPT framing as contributing substantially to GWAS in the introduction? Relatedly in the discussion, it is not clear how having common variants from a fetus would change clinical practice in cases of maternal preeclampsia or gestational hypertension, so perhaps clarifying or removing would be helpful.\nThank you for the comments. About the potential for NIPT to contribute substantially to GWAS, we are imagining the GWAS being done on the mother, using their imputed genome, with QUILT-nipt leading to more accurate imputation, and hence more powerful GWAS. We envision the collection of either survey data, or more ideally linking to EHR, to capture adult quantitative traits and binary traits. This could include (through EHR) things like lab measures, as well as simple traits like height and (pre-pregnancy) weight, as well as things like educational attainment, which are important traits in the field of genetics broadly. For binary disease traits, yes there are some traits, for example heart attacks and breast cancer, that are well studied, whose onset is typically after female maternal age, so would be poorly captured directly, though would be amenable to gwas-by-proxy. However there are many traits of interest with much earlier onset including Type 1 Diabetes, neuropsychiatric disorders (e.g. depression and schizophrenia), auto-immune (IBD, RA), epilepsy, etc, that would make GWAS conducted on phenotypes collected from women of maternal age very valuable.\nTo better rationalize GWAS in the Introduction we have added the following about the value of the phenotypes that can be collected in pregnant women\nOne option is to utilize DNA and phenotypes from pregnant women. Electronic health records (EHR) usually store rich phenotypes for pregnant women, which could include quantitative traits like height, weight, hyperglycaemia 14, as well as socioeconomic traits like educational attainment, as well as disorders with age of onset before maternal age or during pregnancy, such as neuropsychiatric conditions like schizophrenia, as well as gestational diabetes mellitus 15 and intrahepatic cholestasis of pregnancy 16. A natural way to collect DNA from pregnant women is using cell free DNA from non-invasive prenatal testing (NIPT), which is a sensitive and specific screening test which tests for the common fetal aneuploidies of trisomies 13, 18 and 21, and that, according to the American Society of Obstetricians and Gynecologists (ACOG, 2020), should be offered for all pregnancies regardless of risk.\nAnd further rationalized the potential for GWAS in NIPT by specifically highlighting existing GWAS using NIPT\nIn 2018, it was suggested that NIPT had already been performed on ten million pregnant Chinese women 18. Further, NIPT has already been used as a successful means of conducting GWAS, with NIPT based GWAS discovering novel associations for intrahepatic cholestasis of pregnancy 16 , gestational diabetes mellitus 15, metabolites 19, glycemic traits 20 and neonatal phenotypes 21.\nAbout the discussion, we have re-written the text surrounding potential clinical applications, and discuss this in response to other comments. We have put the re-written text below.\nFrom a clinical perspective, relevant traits that are important for health care, such as hyperglycemia, gestational diabetes and birth weight, are genetically influenced by both maternal and fetal genomes 39,42,44–47. Therefore, the imputation of the fetal genome before childbirth may afford the opportunity to change clinical practice as we move toward leveraging PRS in health risk management. Recent papers, which take the uncertainty of genotype imputation into PRS interpretation into account, are of particular relevance here48. Nonetheless, the imputation of specific variants, in particular rare variants, remains challenging. Therefore it is unlikely to be used for clinical diagnosis, and in that context, recent approaches using non-invasive deep fetal exome sequencing are more promising 15,16. We note that QUILT2-nipt should still operate on fetal exome sequencing, so in addition to accurate rare variants, such data would also allow for genotype imputation and PRS calculation.\n2. Relatedly, is using fetal genomes in a GWAS / PRS feasible at current accuracy rates? Variants at 1% MAF are only a little over 60% accurate even at 4x coverage and 0.3 FF. Could it introduce false positives? Most of the significant hits are only nominally significant both for mother and fetus (Fig 5b). Similarly, with the current state of PRS as predictors for traits, does matching array data performance with cfDNA really change the feasibility of its use in the clinic? This is an interesting and valuable analysis, but some text in the discussion addressing this caveat could be merited. Could include the study “Novel insights into the genetic architecture of pregnancy glycemic traits from 14,744 Chinese maternities” by Z. Zhu et. al. (2024, Cell Genomics) in the discussion relating to using NIPT data in GWAS.\nThank you for the comment. We’ll answer minor points first before addressing the larger issue. First the specific quoted accuracy, you’re highlighting 0.01%, for variants at 1%, it is closer to 0.80 accuracy. Second, with respect to Figure 5b, this figure shows the relative equivalency of different NIPT coverages versus arrays. For a real PRS, people would use standard PRS methods like clumping and thresholding or LDpred, which may or may not use variants with marginal p-values below genome-wide significance (most do). Third, semi-related, methods are emerging that take the uncertainty of genotypes into account, e.g. https://pubmed.ncbi.nlm.nih.gov/37490908/, which would increase the benefit of using imputation genotypes.\nFor the larger issue, about whether matching array data performance with cfDNA really changes the feasibility of its use in the clinic. Today, the answer is likely no. However, as GWAS sample sizes increase (hopefully facilitated by this study!), and as we get more accurate methodology for building PRS, so too should the accuracy of array data (and cfDNA) based PRS. We think this study, as the first dedicated method for genotype imputation and hence PRS calculation from lc-WGS cfDNA NIPT, plays an important role in showcasing the utility of this data type for making PRS available for eventual potential clinical applications. We have modified the text of the Discussion and specifically highlighted in this response that future studies will be required to properly evaluate the potential of prenatal PRS.\nWe’ve re-worked the discussion part as follows:\nFrom a clinical perspective, relevant traits that are important for health care, such as hyperglycemia, gestational diabetes and birth weight, are genetically influenced by both maternal and fetal genomes 14,43,46,48–50. Therefore, the imputation of the fetal genome before childbirth affords the opportunity to calculate PRS prenatally, which may impact future health risk management if the PRS is sufficiently predictive, though future studies will be required to fully evaluate this potential. Recent papers, which take the uncertainty of genotype imputation into PRS interpretation into account, are of particular relevance here51. Nonetheless, the imputation of specific variants, in particular rare variants, remains challenging. Therefore it is unlikely to be used for clinical diagnosis, and in that context, recent approaches using non-invasive deep fetal exome sequencing are more promising 22,23. We note that further extension to QUILT2-nipt should still operate on fetal exome sequencing, so in addition to accurate rare variants, such data would also allow for genotype imputation and PRS calculation.\nMinor comments:\n1. The manuscript would benefit from a careful read through to streamline and improve readability as there are some very complex sentences that are hard to parse. Some examples:\nThank you for the suggestions. We’ve now corrected many typos and updated the text\na. The first sentence of results lacks clarity: QUILT2 is based upon, and improves, our previous model QUILT1 for diploid genotype imputation, which operates, which uses an iterative approach, as follows. Suggest to remove “which operates” from the sentence.\nThank you. We removed “which operates”\nb. There are two i.e.’s in this 55-word sentence. It would be more readable if broken up or streamlined: “However, state-of-the-art imputation methods for lc-WGS are designed for diploid samples, i.e. where the sequencing reads come from two haplotypes 3–6, and there is no dedicated method designed for utilizing the triploid nature of NIPT data, i.e. where the sequencing reads come unequally from the maternal transmitted, maternal untransmitted, and paternal transmitted haplotypes.”\nThank you. We have removed “i.e”s and updated the sentence.\nHowever, state-of-the-art imputation methods for lc-WGS are designed for diploid samples, and there is no dedicated method designed for NIPT data taking both the fetus DNA and mother DNA into account.\nc. “For diploid samples, we used the WGS derived UKB panel to show that QUILT2 is approximately as accurate as QUILT1, but substantially faster, and uses much less memory, and as such, should be much less expensive in cloud based applications.” Example streamlined suggestion: For diploid samples, the WGS-derived UKB panel showed that QUILT2 is as accurate as QUILT1 but significantly faster, more memory-efficient, and cost-effective for cloud-based applications.\nWe adopted your version. Thank you!\nd. “In addition, while QUILT2 is slower than GLIMPSE2 when processing many samples, it is faster, by an order of magnitude on very large reference panels, when processing only a very small number of samples.” Suggestion: While QUILT2 is slower than GLIMPSE2 for many samples, it is an order of magnitude faster when using very large reference panels to process a small number of samples.\nWe adopted your version. Thank you!\ne. Line 99-101: suggest that the sentences “Because of the second part of the iterative process, QUILT1 has linear time complexity in the number of haplotypes and variants in the full reference panel. For RAM and speed reasons…” be rewritten to “Given the two-part iterative process, QUILT1 has linear time complexity dependent on the number of haplotypes and variants in the full reference panel. To minimize RAM and runtime…”\nYour suggested alternate sentence has a different meaning to the original sentence. It’s not the iterative process itself, but rather specifically the second part, that yields linear computational complexity in the number of reference haplotypes (this is the key bit, both the first and second part yield linear time complexity in the number of SNPs). We have left this sentence unchanged, as we think it’s less of a problem compared to the others you highlighted.\n2. Including the number of SNPs as another factor in the runtime / memory usage could be interesting. For example, subset the SNPs in the lc-WGS data and see how it varies with the number of SNPs included. There is some analysis on this for the QUILT-nipt method vs QUILT2, though it only focuses on all variants vs only common, and does not mention the number of SNPs included at these thresholds.\nThank you for the comment. We have partially done this, in that in Supplementary Figure 5, we evaluated the runtime and memory by adjusting a threshold which determines at which frequency SNPs are considered common vs rare in QUILT2. More generally, usually, the number of SNPs in a panel is a more minor concern than the number of haplotypes. In Figure 3, where we show the scaling properties of the methods, the number of SNPs in the panels varies as well, as monomorphic SNPs are removed from smaller panels. So in that sense, we do evaluate how the number of SNPs effects run time, and we don’t think another dedicated analysis would be particularly fruitful for readers of the manuscript.\n3. Given that low-coverage WGS imputation methods often differ from GWAS array imputation methods in that the former relies on genotype probabilities whereas the latter relies on genotype calls, is QUILT2 incompatible with GWAS array data? This point is worth clarifying to ensure use cases are clear to readers.\nQUILT2 as written uses sequencing reads directly. It could be modified to use array genotype information, but there are many dedicated methods already to do this. We note that QUILT2 operating on sequencing reads is a distinct advantage of the method, as we’ve shown in the paper, e.g. with the long read analysis, in contrast with GLIMPSE2, which only uses input information in a per-SNP manner.\nWe have amended the Introduction in general, and have included this sentence below. Implied but not stated is that QUILT cannot use array input. We can make this more explicit if the reviewer would like us to.\nQUILT and GLIMPSE further differ in how they treat input data, with GLIMPSE using per-SNP genotype likelihoods, while QUILT uses and phases the sequencing reads directly, which is preferable for long-read data.\n4. Common variant definition of MAF=0.1-0.2 seems a little funny. Why not just MAF > 0.1?\nWe stratified variants into bins based on minor allele frequency, as is standard in methods papers like these, as variant frequency is strongly correlated with imputation. This is most important for variants at low frequency. We were perhaps overzealous in our bin choices for common variants. Nonetheless more bins at higher frequencies makes for slightly nicer figures, and reporting this way ensures consistency between figures and text. If the reviewer remains keen we can change this e.g. to 5 or 10 to 50%, but it is unlikely to affect a qualitative interpretation of the results. For reference, we used 20-50% for the QUILT1 paper.\n5. Given that QUILT2 imputed more accurately than GLIMPSE2 with UKB WGS but not phase 3 1kGP and HRC, is there better genotype error tolerance in GLIMPSE2? Notably, a 30x WGS version of 1kGP is available.\nThank you for the comment. Yes we do use the higher coverage 30X 1kGP panel. It is tricky to know for sure what’s driving the differences, which are minimal. You’re right that it might have to do with better genotype error tolerance and how this relates to the PBWT based heuristics.\n6. In Figure 2a, it is hard to tell the coverage lines apart because points overlap. The authors might consider using a color gradient for increasing coverage, with QUILT2 vs GLIMPSE instead indicated by line type (e.g. dashed vs solid). Comparing accuracy across Figure 2a, 2b, and 2c is challenging as the scale is changing.\nThank you for the comments. We now re-designed the plot and added 4x in the benchmarking as suggested by another reviewer. We hope this will facilitate interpretation.\n7. Very interesting to see the ancient DNA imputation results. Though in Figure 2d (and corresponding supplemental figures) - is there a particular reason to switch from MAF as x-axis (for figures 2a-c) to coverage as x-axis and color as MAF in Fig2d? Could consider switching the axis to match the other plots, though not critical to do this.\nThank you for the comment. We agree this was confusing. We have re-designed it now to be more in line with the other figures.\n8. Non-reference concordance is a bit more of a forgiving metric than aggregate r^2. Do the authors care to comment on the relevance of this in the comparison of aDNA vs modern samples?\nWe agree, but here we are following the style in the GLIMPSE2 paper and other papers on ancient DNA.\n9. In Figure 3, the legends are a bit confusing, since those shown in A differ from in B but are relevant to A-C. I would suggest extracting a common legend to the bottom of this figure, then reducing the y-axis for panel b, as it is hard to distinguish between lines. The description of memory usage in the text describing this figure could be clearer about the relative comparison between QUILT2 and GLIMPSE2 (note that GLIMPSE2 is not mentioned here): “We note that QUILT2 uses slightly higher RAM given 100 samples to run, and QUILT2 also runs faster with lower RAM for lower coverage data (Figure 3C, Supplementary Figure 5, 6).”\nWe appreciate this suggestion. We have re-designed the figure now. In addition, we have changed the text to be clearer about the relative comparison of RAM usage with sample size for GLIMPSE2 and QUILT2.\nRAM usage is approximately constant for GLIMPSE2 with sample size, whereas it increases in QUILT2 (Figure 3). In addition, QUILT2 also runs faster for lower coverage data (Supplementary Figure 5, 6). We do note that both methods have reasonably low RAM usage, which we don’t estimate to be a hindrance for most users.\n10. The NIPT comparisons in Figure 4 are nice. The authors might want to comment on whether imputing only common fetal variants is the best strategy given the focus on embryonic screening for structural variants and potentially Mendelian variants, or whether lcWGS is an uncommon current application but with approaches like QUILT2 enabling more.\nThank you for the comment. Sequencing cfDNA and NIPT has been proven to be the most accurate method for aneuploidy detection, and it is increasingly cheap and easy to perform. Specifically NIPT using genome-wide lcWGS, vs targeted sequencing, has been adopted at scale with many millions of samples already performed in both the Netherlands and China. We hope that, with this paper, people will also see the benefits of imputing the mother and the fetus for both research and clinical opportunities. We also hope that, as we’ve shown, accuracy improves as sequencing coverage improves, and in such cases, we expect to see further improvement in the ability to genotype Mendelian (rare) variants. Further, for structural and copy number variants, we’ve mentioned this as an area for further improvement in the Discussion. Finally for rare variants, imputing these, at low coverage, is very difficult, and likely insufficient for clinical purposes, hence why in the Discussion we highlight alternative approaches like non invasive fetal exome sequencing. We think we’ve covered these issues in the Discussion in the second and third paragraphs though if the reviewer feels we could do more to address them we are happy to add more text in.\n11. Any thoughts on why there is a downward trend in aggregated R2 for common SNPs in Fig 4, both maternal and fetal genomes? More pronounced for diploid than NIPT and more pronounced with increasing coverage.\nThank you for the comments, which are very good points. The aggregated R2 starts dropping towards high frequency variants can be explained by several factors. First, for high frequency variants, it is more likely that the paternal alleles/reads are different to the maternal allele/reads, which will be more pronounced given high depth sequencing. Then, when the paternal alleles are different to the maternal alleles, the ability to distinguish them in the model explains how the results will look like. With the diploid model, where the paternal origin is wrongly modeled as either the maternal transmitted or maternal untransmitted reads, it just can not distinguish them and hence will be confused. As an example, given the truth allele of maternal untransmitted, maternal transmitted and paternal origin is 1, 1, and 0 respectively (hence the mother is 1|1 homozygous alternatives), if the fold difference between maternal untransmitted (1) and paternal origin (0) is very low, say 2-fold (means fetal fraction is very high), then there is 1 / 3 chance to call the mother 1|0, which is wrong. In contrast, if the fetal fraction is very low, say fold change in the coverage of maternal untransmitted vs paternal origin is 9-fold, then there is only 1/10 chance calling the mother genotype wrong. This downward trend in the diploid model was also observed by other studies, e.g the recent one https://www.sciencedirect.com/science/article/pii/S2211124724011501#sec3.\nNow, for the nipt mode, relatedly, even though we model the paternal origin reads as a third distinct haplotype, the higher fetal fraction can still blur the distinction between the paternal haplotype and the maternal untransmitted haplotype, making it easier for the model to confuse them e.g. at ff = 10%, we expect 45% maternal untransmitted vs 5% paternal origin reads, a nine-fold difference, while at ff = 30%, we expect 35% mat-u vs 15% pat, around 2-fold difference. As an example, given the truth genotype of mother, father and fetus is 0|1, 1|1 and 1|1 respectively, if the fold change between maternal untransmitted (0) and paternal origin (1) is only 2 compared to 9 (when fetal fraction is low), then it challenges the model more in distinguishing them, hence it might be a higher chance to call the fetus as 0|1 instead of 1|1.\n12. Why were only quantitative traits used for the GWAS / PRS analyses? Do the authors expect a change in behavior for logistic mixed models? Don’t necessarily have to test binary traits, but just include a sentence explaining why quantitative should behave the same as binary traits.\nWe expect differences in power between array vs lcWGS for quantitative traits to be similar for binary ones. We chose to use only quantitative traits so that all traits could be analysed in the same way. We have added a sentence to the results to make it clearer that we expect binary traits to show broadly the same results\nWe only tested quantitative traits but expect binary traits to show broadly similar results.\n13. For the NIPT GWAS locus identification part, the authors use a loose definition of p < 0.05 to describe identification. Figure 5B shows that the vast majority are identified with nominal significance, so a concordance check on effect size (same direction at least?) may also be warranted. It might be helpful to reorder these bars so that the mother array followed by mother with decreasing coverage are shown, followed by the simulated fetus, as it’s otherwise hard to compare the relative change in performance. Blue and orange coloring in Figure 5B should be switched to not be confused with the blue/orange for fetus/mother color scheme in other plots. In terms of relative gains in discoverability expected across technologies, this publication might be worth discussing: Gaynor, S.M., Joseph, T., Bai, X., Zou, Y., Boutkov, B., Maxwell, E.K., Delaneau, O., Hofmeister, R.J., Krasheninina, O., Balasubramanian, S., et al. (2024). Yield of genetic association signals from genomes, exomes and imputation in the UK Biobank. Nat. Genet. 56, 2345–2351.\nThank you for the comments. First, thank you for suggesting a better way to present the results in Figure 5. We’ve now re-designed it as your suggestion. Second, yes, a p-value threshold of 0.05 is lenient, however the main purpose of this is as a means to show that results are generally consistent when doing using arrays, or lc-WGS using NIPT. As an alternative to checking concordance of effect size direction, below we check the Pearson correlation between effect sizes versus arrays, and we see they are incredibly similar, ranging from 0.998, to 0.935. As such, any similarity in the number of loci shown in Figure 5B is undoubtedly from the high similarity in underlying genotypes.\n14. Given the focus on NIPT, how many common clinically actionable HCMG genes or other common variants relevant to fetal Mendelian disease (e.g. CFTR deltaPhe508) are well-imputed? Is low-coverage WGS the best strategy here, or might a WEGS or blended strategy that enriches for higher coverage at interpretable coding variants be worth discussing?\nThank you for the comments. In the paper, we didn’t look into the specific common variants that are relevant to fetal Medelian disease. Assuming for example the CFTR deltaPhe508 variant shows similar imputation accuracy as for other variants at similar frequency (1.5% in Europeans by gnomAD), then given 1x, 2x, 4x and 10x NIPT lcWGS with 10% FF, the imputation accuracy r2 for the fetal genotype would be around 0.713, 0.809, 0.866, and 0.885 respectively. This may be sufficient to request confirmatory sequencing, or in the context of pre-existing knowledge of being a carrier of the variant, but in itself does not meet reasonable diagnostic thresholds. As such, we’ve modified our discussion to emphasize this, and the value of recent development of non-invasive fetal exome sequencing.\nNonetheless, the imputation of specific variants, in particular rare variants, remains challenging. Therefore it is unlikely to be used for clinical diagnosis, and in that context, recent approaches using non-invasive deep fetal exome sequencing are more promising 22,23. We note that further extension to QUILT2-nipt should still operate on fetal exome sequencing, so in addition to accurate rare variants, such data would also allow for genotype imputation and PRS calculation.\n15. The discussion point about the clinical actionability of de novo CNVs e.g. 22q11.2 is interesting but also confusing, as imputation methods like QUILT2 use a reference panel to only impute existing variants in the population. This is framed as a potential future development opportunity for QUILT2 that seems highly speculative at present as a class of variation undetectable by any imputation method.\nPeople have previously called CNVs including de novo ones using array data e.g. PennCNV. People have also looked at calling more complicated chromosomal aberrations using array data from blood e.g. work by Po-Ru Loh and colleagues using the UKB. That being said, QUILT2 is one of the very few imputation methods that works directly on sequencing reads, rather than genotype probabilities / genotype likelihoods. With QUILT2, Gibbs sampling on read labels yields read label estimates that directly assign reads to haplotypes. As such, it is also a read phasing method. We have a new paper in preparation describing improved SV genotyping accuracy by using the per-read phasing information already output from QUILT2, as this is additional evidence, beyond what current methods use, to generate a likelihood for whether an SV call is correct or not.\nIn the future, a simple modification to the QUILT2 model could be developed to consider a “copy number”. For recurrent microdeletions and duplications, with known start and end points, this would not be difficult, as the currently used prior probability on read labels (based on fetal fraction) would be modified to include copy number, given known start and end points. One could then for instance compare the likelihood of a normal model (with normal copy number) to other models (e.g. paternal origin duplication between LCR A and LCR D in 22q11), and “call” the event if a sufficient difference in likelihood were present. For general de novo CNVs, more complicated future development that could adopt varying copy number (e.g. as a hidden variable itself) would permit a general purpose solution, however that would be more complicated to implement. Even further in the future there could be scope to implement a read phasing and local realignment model, like used for variant calling in Platypus and other more modern variant callers e.g. the GATK HaplotypeCaller.\nAs such we think it’s a reasonable speculation to motivate interest in this result and future work based on this method or using these ideas for future methods.\n16. Methods: “𝜆 is the parameters of our model” - this parameter needs to be defined more clearly.\nWe have changed this to\n𝜆 is the parameters of our model, including the recombination rate and reference haplotypes\nWe also mention later another parameter “FF is the fetal fraction, which is required as an input and can be obtained from NIPT report,”\n17. Line 623: “Real big” reference panel - change to “very large” or “one of the largest”\nThank you, we have corrected the language.\n18. Not sure if it’s the pdf formatting, but for the notation in the supplemental note, there is no semi-colon.\nThank you for this. We wrote “semi-colon” but we meant colon. We have corrected this in the supplement.\nReviewer #4 (Remarks on code availability):\nThe README provides detailed instructions, including conda environment installation, test run information, and tutorial pages for the various use cases of QUILT2. A number of issues have been opened on the github, which is normal with any bioinformatics software. Importantly, the authors seem responsive to those ongoing issues. We did not try to install the software, but based on the github repository, the usability seems high.\nThank you for your kind comments.\nReviewer #5 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nR1\nI would like to add a further comment regarding large-scale GWAS analysis. I understand the rationale for splitting jobs across samples, which provides a runtime advantage for QUILT2 when dealing with small sample sizes. In your response to Reviewer 3 regarding large scale GWAS analysis, you note: “As such, the issue of running large datasets with many thousands of samples is really just an issue of how to most efficiently run per-sample, which is why we’ve focused as we have on comparatively small sample sizes like shown in Figure 3.”\nIf I understand correctly, the authors suggest splitting large datasets both into chunks and by sample. My concern is that this approach introduces substantial overhead in starting and managing multiple machines, as it generates a large number of tasks per imputation job, potentially offsetting any runtime advantage. While fast per-sample processing can be beneficial, higher parallelisation does not make it automatically faster in a cloud environment (though it could be advantageous on a local Slurm cluster). I am curious whether the authors could provide an overall estimate of cloud-based imputation performance (e.g., using a larger sample size for both methods -> runtime, costs). By comparison, current state-of-the-art imputation servers using microarray data typically use chunking only, which works efficiently in the cloud. Even when using spot instances, which may be interrupted more frequently, the cost advantage is substantial (I think up to 90% cheaper than on-demand instances). I would also appreciate clarification on whether the Snakemake pipeline supports per-sample chunking, as this is not as straightforward for end-users as region-based chunking (due to splitting and merging in correct order). Finally, I am unfamiliar with the distinction between “high” and “low” priority jobs. In the UK Biobank instance rating card, the main difference is between spot and on-demand instances (0.2336 vs 1.2832 GBP/hour for the instance type used), not sure if these terms can be used interchangeably\nre: estimate of cloud-based imputation performance: We have presented what we think is a reasonable strategy for imputing samples, and this costing includes the overhead cost of starting VMs. In the UKB environment, we suggest running a reasonably similar number of samples (N = 32) as the GLIMPSE2 authors (N = 100) samples. This costing and setup reflects what would be most important, the total cost, when running a large number of samples. It is certainly true that in different computer architectures, including cloud vs academic, as well as different sample size setups, e.g. running very few samples, a different split between sample chunking and region chunking might make more sense. But overall we think we’ve presented a comparison that accurately reflects the main issues that users would experience and is fair to both QUILT2 as well as GLIMPSE2.\nre: Snakemake: the Snakemake pipeline does not currently perform per-sample chunking. However this is generally easier to do compared to per-region chunking. Since the imputed SNPs will be the same for all samples, merging VCFs afterwards is easy using tools like the GATK MergeVCFs or BCFTOOLS merge.\nRe “high” vs “low”: When submitting a job, the VM is assigned depending on the specified priority, which takes three values: high, normal or low. High priority jobs will use on demand instances, while low priority jobs will use spot instances. A normal priority job running on a spot VM can be interrupted and after several restarts it will turn into a high priority job running on an on-demand VM. More on this\nhttps://dnanexus.gitbook.io/uk-biobank-rap/administrator/managing-usage-and-storage-costs\nand https://dnanexus.gitbook.io/uk-biobank-rap/working-on-the-research-analysis-platform/managing-jobs/managing-job-priority.\nR2\nThanks to the authors for addressing my comments.\nGiven their reply and that they only showed that their approach partly may work in an NIPT context (i.e. only one sample). If the editor decides to publish this manuscript, I recommend to remove NIPT (fetal cell-free DNA) in the abstract and the title and instead limit their claims as potential application of QUILT2 in an NIPT context. The low imputation capability and testing QUILT2 on only one actual NIPT sample, is far from reality.\nIn my opinion, QUILT2 is far from \"enabling imputation of the maternal and fetal genome using cell free non-invasive prenatal testing (NIPT) data.\"\nGiven the editorial need to make the title more succinct, as well as a need to remove punctuation from the title, we have removed the specific reference to cfDNA NIPT imputation in the title. However we have left it in the abstract, as this is a major feature of this work. QUILT2 is the only method that performs imputation on cfDNA NIPT data, and allows for imputation of the mother and the fetus both, and also improves the accuracy of imputation for both, with the accuracy of the maternal imputation being quite high. As NIPT sequencing is performed millions of times a year, we think it is important to highlight this part of QUILT2 in the abstract."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have provided clear and thoughtful responses to my comments (Reviewer #1), addressing both major and minor concerns. There are no further comments from my side.\nNevertheless, I was not fully convinced by their reply to Reviewer 3’s comment on “Lack of reference for major common use cases,” specifically regarding the timing benchmark. Please find below my reply to this point:\n----------------------------------------------------------------------\nREPLY TO ANSWER “Lack of reference for major command use cases” regarding the Timing benchmark:\nI would like to add a further comment regarding large-scale GWAS analysis. I understand the rationale for splitting jobs across samples, which provides a runtime advantage for QUILT2 when dealing with small sample sizes. In your response to Reviewer 3 regarding large scale GWAS analysis, you note: “As such, the issue of running large datasets with many thousands of samples is really just an issue of how to most efficiently run per-sample, which is why we’ve focused as we have on comparatively small sample sizes like shown in Figure 3.”\nIf I understand correctly, the authors suggest splitting large datasets both into chunks and by sample. My concern is that this approach introduces substantial overhead in starting and managing multiple machines, as it generates a large number of tasks per imputation job, potentially offsetting any runtime advantage. While fast per-sample processing can be beneficial, higher parallelisation does not make it automatically faster in a cloud environment (though it could be advantageous on a local Slurm cluster). I am curious whether the authors could provide an overall estimate of cloud-based imputation performance (e.g., using a larger sample size for both methods -> runtime, costs). By comparison, current state-of-the-art imputation servers using microarray data typically use chunking only, which works efficiently in the cloud. Even when using spot instances, which may be interrupted more frequently, the cost advantage is substantial (I think up to 90% cheaper than on-demand instances). I would also appreciate clarification on whether the Snakemake pipeline supports per-sample chunking, as this is not as straightforward for end-users as region-based chunking (due to splitting and merging in correct order). Finally, I am unfamiliar with the distinction between “high” and “low” priority jobs. In the UK Biobank instance rating card, the main difference is between spot and on-demand instances (0.2336 vs 1.2832 GBP/hour for the instance type used), not sure if these terms can be used interchangeably.", "reviewer_2": "Thanks to the authors for adressing my comments.\nGiven their reply and that they only showed that their approach partly may work in an NIPT context (i.e. only one sample). If the editor decides to publish this manuscript, I recommend to remove NIPT (fetal cell-free DNA) in the abstract and the title and instead limit their claims as potential application of QUILT2 in an NIPT context. The low imputation capability and testing QUILT2 on only one actual NIPT sample, is far from reality.\nIn my opinion, QUILT2 is far from \"enabling imputation of the maternal and fetal genome using cell free non-invasive prenatal testing (NIPT) data.\"", "reviewer_4": "The authors have been responsive to my (and seemingly other) reviewer comments, and I support the advances of QUILT2 moving towards publication.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": ""}]

Decision letter