Dear editor We write to raise conceptual and methodological concerns regarding the article by Lu et al. (‘Impact of control selection strategies on GWAS results: a study of prostate cancer in the UK Biobank’) [1]. A central claim of the paper is that nearest-neighbour propensity score matching represents an effective, or even superior, strategy for control selection in genome-wide association studies (GWAS). We contend that this conclusion is unsupported and the core of the issue is a fundamental mismatch between the methodology employed and the aims of GWAS. Nearest-neighbour propensity score matching is a tool developed for causal inference, where the objective is to estimate treatment effects while mitigating confounding, particularly when matching on categorical phenotypes, such as sex or socioeconomic factors. GWAS, in contrast, is an association-based framework that estimates marginal Single Nucleotide Polymorphism (SNP)–phenotype relationships using regression or mixed models that can explicitly adjust for confounding variables such as age, technical factors, and genetic principal components (PC). Within this context, matching is neither required for unbiased estimation nor theoretically motivated as a primary design strategy. Conditioning on a variable (genotype PC) that directly influences both the exposure (genotype) and outcome (disease risk) creates a collider bias. This results in control individuals who are genetically similar to cases, thus attenuating or eliminating precisely the subtle allele-frequency differences that GWAS aims to detect. Modern GWAS tools such as REGENIE use linear mixed models, which can account for all degrees of relatedness and population structure across all the samples simultaneously. Thus, using propensity score matching on the same covariates further amplifies the overfitting effect, resulting in groups that can be so similar that the actual causal variants no longer show a statistical difference. This risk is not merely theoretical. Recent work has shown that even the standard practice of including PCs as covariates in the association analysis can induce collider bias in certain settings [2]. More aggressive conditioning on genome-wide similarity is therefore likely to exacerbate rather than mitigate such bias. These limitations extend to rare-variant and gene-burden analyses, where the manuscript suggests matching may be particularly advantageous. Because rare variants are strongly structured by the same fine-scale ancestry captured by genome-wide SNPs, matching implicitly conditions on the expected rare-variant spectrum and attenuates genuine signals precisely where power is already limited. Furthermore, by artificially restricting the dataset, this directly impacts the ability to identify rare variants with small effects. In summary, nearest-neighbour matching using genotype PCs, is conceptually misaligned with the goals of GWAS. By utilizing a causal-inference tool unsuited to association testing and introducing a collider variable, the approach risks removing true biological signal, introducing bias, and sacrificing power, both for polygenic as well as rare-variant architectures. We urge the authors to clarify or correct these points to prevent misinterpretation by the field. Propensity score matching is conceptually inappropriate for GWAS. Matching reduces power and using genetic principal components can induce collider bias. Risks from bias are empirically supported and may be worsened by aggressive matching. None declared.

To the Editor, The Isaac pipeline described in the Bioinformatics article ‘Isaac: ultra-fast whole-genome secondary analysis on Illumina sequencing platform’ (Raczy et al., 2013) has been used in cancer sequencing studies (Burns et al., 2018; Quigley et al., 2018) and in the ongoing UK 100 000 Genomes Project (100KGP) (Turnbull et al., 2018). Whilst Isaac has been benchmarked with respect to variant calling (Raczy et al., 2013), there has been less extensive evaluation of its suitability for other analyses routine to cancer genomics. Estimating the fraction of cancer cells with individual somatic mutations is central to cancer genome studies, including characterization of clonal architecture (Dentro et al., 2017). Estimation of these cancer cell fractions (CCFs) is however contingent on unbiased assessment of the fraction of reads supporting variant allele frequencies (VAFs). We demonstrate that VAFs computed by Isaac are biased by the preferential soft clipping of reads supporting non-reference alleles, with deleterious consequences on downstream analyses reliant on unbiased CCF estimation. Reads supporting heterozygous single-nucleotide polymorphism (SNP) reference and alternate alleles can be expected to occur with equal probability when sequencing normal tissue. Due to limited sequencing depth, the exact number of reads supporting reference and alternate alleles will not be equal in many instances, even when the aligner is unbiased. It can be expected however that the median VAF of a large number of heterozygous SNPs will be 0.5. We assessed heterozygous SNP VAF distributions in whole genome sequencing (WGS) data from the germline of 25 multiple myeloma (MM) tumor-normal pairs aligned to GRCh38Decoy assembly using Isaac v03.16.02.19. Germline variants were called using Starling v2.4.7 (Raczy et al., 2013) and VAFs were calculated directly from alignment files using alleleCount (Van Loo et al., 2010). Median VAFs per sample ranged from 0.478 to 0.479 (Fig. 1A), with this consistent skew indicating that Isaac can exhibit bias toward the reference allele. Evidence of reference bias from Isaac. (A) Heterozygous SNP VAF distributions in sequencing data from 25 normal samples. Dashed line represents expected median VAF of 0.5. Whiskers extend 1.5 times inter-quartile range and values outside of this range are not shown. (B) Proportion of reads covering SNP positions supporting the reference (blue line) and alternate (red line) alleles with that read position soft clipped. (C) SNV VAFs from 25 tumor-normal pairs. (D) CCFs of clonal mutation clusters identified by DPClust. Blue and grey dashed lines denote the median putative clonal mutation cluster CCF and a CCF of 1, respectively. (E) Ccube sample purity estimates. Distribution differences assessed using Wilcoxon Signed-Rank test Evidence of reference bias from Isaac. (A) Heterozygous SNP VAF distributions in sequencing data from 25 normal samples. Dashed line represents expected median VAF of 0.5. Whiskers extend 1.5 times inter-quartile range and values outside of this range are not shown. (B) Proportion of reads covering SNP positions supporting the reference (blue line) and alternate (red line) alleles with that read position soft clipped. (C) SNV VAFs from 25 tumor-normal pairs. (D) CCFs of clonal mutation clusters identified by DPClust. Blue and grey dashed lines denote the median putative clonal mutation cluster CCF and a CCF of 1, respectively. (E) Ccube sample purity estimates. Distribution differences assessed using Wilcoxon Signed-Rank test Isaac has a parameter (–clip-semi-aligned) that invokes the soft clipping of reads at each end until a stretch of five consecutive bases are matched with the reference sequence (here, we term this ‘alt-clipping’ to distinguish from soft clipping performed for other reasons). This parameter is present in all Isaac versions after v01.13.06.20 and was used to align the 25 normal samples. To test whether alt-clipping is responsible for the reference bias exhibited by Isaac, we re-aligned the 25 normal samples with Isaac without alt-clipping. When alt-clipping was not performed, the median VAF of heterozygous SNPs in each sample equaled 0.500 (Fig. 1A), thereby showing that alt-clipping introduces reference bias. Alt-clipping results in the clipping of the majority of reads supporting the alternate allele where the variant position is within five bases of either read end (Fig. 1B). Fewer reads supporting the reference allele are soft clipped and VAFs therefore become biased towards the reference allele. If the preferential soft clipping of reads supporting the alternate allele is responsible for the reference bias, then we would expect the effect to be negated if the ends of all reads were soft-clipped by five bases. To further validate the effect of alt-clipping, we therefore soft-clipped five bases at each end of all reads in the alt-clipped alignments (here, we term this ‘balanced-clipping’). The median VAF of heterozygous SNPs in balanced-clipped alignments equaled 0.500 in each sample (Fig. 1A), demonstrating that the preferential clipping of reads supporting the alternate allele introduced by alt-clipping causes reference bias. To assess the effect of alt-clipping on the analysis of cancer genomes, we aligned the tumor MM WGS data using Isaac with and without alt-clipping. Somatic single nucleotide variants (SNVs) were called using Strelka v2.4.7 (Kim et al., 2018). Unlike SNPs in normal samples, we do not know the true VAF of SNVs in tumor samples, as they can be affected by copy number aberration, normal sample contamination and clonal heterogeneity. SNV VAFs from alt-clipped alignments were however lower than SNV VAFs from the same samples from alignments generated without alt-clipping (P < 2.2 × 10−16; Fig. 1C), indicating that alt-clipping also affects somatic SNV VAFs. To test the effect of alt-clipping on subclonal reconstruction, we ran Battenberg and DPClust, which uses a Dirichlet process to model subclonal fractions (Nik-Zainal et al., 2012). If SNV VAFs do not exhibit allelic bias, we would expect DPClust to identify clusters of mutations with CCFs centered on 1, representing clonal mutations. When DPClust was run using alignments generated with alt-clipping, the median CCF of putative clonal mutation clusters (defined as the cluster with a CCF closest to 1) was 0.959, compared to 0.983 when run using alignments generated without alt-clipping (Fig. 1D). Finally, we assessed the effect of alt-clipping-induced reference bias on tumor sample purity estimation. Sample purities estimated using Ccube (Yuan et al., 2018) were smaller when computed using alt-clipped alignments than non-alt-clipped alignments (P = 1.3 × 10−3; Fig. 1E), demonstrating that alt-clipping also affects purity estimation. Reference bias introduced by Isaac through alt-clipping can affect downstream processes, potentially making conclusions unreliable for many types of cancer analysis. If unbiased VAFs are required, Isaac should be run with soft clipping of semi-aligned reads disabled, or an alternative aligner such as BWA (Li and Durbin, 2009) should be used. Although realignment can be performed where clipped alignments have been previously produced, this may be cost or time-prohibited. For example, projects such as 100KGP have already sequenced and aligned >10 000 tumor-normal genome pairs. In such cases, equally clipping all reads would enable downstream analyses reliant on unbiased VAFs without the need for sequencing data realignment. While the Isaac aligner version assessed in this study (v03.16.02.19) was released in April 2016, as of November 2019 it is still being used with reference-bias-introducing alt-clipping in 100KGP. Whether reference bias has affected previous studies using the Isaac aligner is difficult to predict. It is clearly essential that aligners, such as Isaac, be evaluated to ensure that the data they produce are not systematically biased. This work was supported by grants from Cancer Research UK [grant number C1298/A8362], Myeloma UK and a David Forbes Nixon Foundation Fellowship (to M.K.). Conflict of Interest: none declared.

MOTIVATION: Sufficiently powered case-control studies with next-generation sequence (NGS) data remain prohibitively expensive for many investigators. If feasible, a more efficient strategy would be to include publicly available sequenced controls. However, these studies can be confounded by differences in sequencing platform; alignment, single nucleotide polymorphism and variant calling algorithms; read depth; and selection thresholds. Assuming one can match cases and controls on the basis of ethnicity and other potential confounding factors, and one has access to the aligned reads in both groups, we investigate the effect of systematic differences in read depth and selection threshold when comparing allele frequencies between cases and controls. We propose a novel likelihood-based method, the robust variance score (RVS), that substitutes genotype calls by their expected values given observed sequence data. RESULTS: We show theoretically that the RVS eliminates read depth bias in the estimation of minor allele frequency. We also demonstrate that, using simulated and real NGS data, the RVS method controls Type I error and has comparable power to the 'gold standard' analysis with the true underlying genotypes for both common and rare variants. AVAILABILITY AND IMPLEMENTATION: An RVS R script and instructions can be found at strug.research.sickkids.ca, and at https://github.com/strug-lab/RVS. CONTACT: [email protected] SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Summary: SNPLINK is a Perl script that performs full genome linkage analysis of high-density single nucleotide polymorphism (SNP) marker sets. The presence of linkage disequilibrium (LD) between closely spaced SNP markers can falsely inflate linkage statistics. SNPLINK removes LD from the marker sets in an automated fashion before carrying out linkage analysis. SNPLINK can compute both parametric and non-parametric statistics, utilizing the freely available Allegro and Merlin software. Graphical outputs of whole genome multipoint linkage statistics are provided allowing comparison of results before and after the removal of LD. Availability: SNPLINK is freely available for non-commercial research institutions. For full details see www.icr.ac.uk/cancgen/molgen/MolPopGen_Bioinformatics.htm Contact: [email protected]