Benchmarking DNA Sequence Models for Causal Variant Prediction in Human Genetics

We appreciate your comments and suggestions. We have uploaded a revised version of our manuscript.

The numbers of both variants and traits are small.

The causal variants of Mendelian traits were collected from a single source (Smedley et al., 2016). The causal variants of complex traits were collected from also a single source (UKBiobank). It would be more meaningful if the authors could incorporate more complex traits and disorders with well-defined causal variants, such as those from FinnGen or other large BioBanks, to improve the benchmark’s generalizability.

The authors need to address how biases in the single dataset curation might affect benchmarking results and suggest methods for mitigating these biases, such as incorporating data from additional cohorts or stratified sampling.

Thank you for your feedback. We do not believe that the number of traits (196 in total) is considered small in the field. We now include expanded datasets with millions of negative variants, but we note that the number of putative causal variants remains very low. However, we still believe that benchmarking on these putative causal non-coding variants has provided numerous insights throughout the paper. We have done our best to incorporate as many trait-associated variants as possible. We are not aware of any larger publicly available resource of putative causal non-coding variants for Mendelian traits. We also note that the OMIM variants from Smedley et al. (2016) are themselves gathered from the global literature, not from a single source. We envision our data curation and benchmarking as a continued effort, and we hope to add more variants in the future. FinnGen variants, unfortunately, do not have LD scores publicly available yet. The UK BioBank has been a pioneer and contains the largest amount of public statistics available, and was therefore our dataset of choice to kickstart this project.

The way for the benchmark study to define the causal variants is too simple and might be controversial. Both OMIM pathogenic annotation and PIP>0.9 do not imply the variant is causal.

Thank you for your suggestion. We do clarify that these are putative causal variants:
“In this work, we present TraitGym, a curated dataset of genetic variants that are either known to be causal or are strong candidates across 113 Mendelian and 83 complex traits, along with carefully constructed control variants.”
and specifically in the task definition:
“Task definition. The task is to classify whether a variant is putatively causal for any trait or not.”

In some sections, we might speak of causal variants as a shorthand, as this is the norm in the literature. For example, Smedley et al. (2016) also use “causal” in their abstract to refer to the exact same OMIM variants. In the case of complex traits, concurrent work by expert statistical geneticists (Fabiha et al., 2016) also refers to PIP > 0.9 variants as “causal”.

Many baseline methods were missed. The authors may refer to a very recent variant effect prediction benchmark study (doi: 10.1186/s13059-024-03314-7). 24 methods for variant effect prediction were benchmarked. In another benchmark study in 2023 (doi: 10.15252/msb.202211474), 55 variant effect predictors were benchmarked. Many of those methods were ignored in this study.

We appreciate your suggestion, but we would like to note that the great benchmarks mentioned above concern variant effect prediction for coding variants only, while our study is precisely aiming to fill the gap in non-coding variant effect prediction.

The evaluation metric used by the benchmark study is too simple based on the binary classification problem. In most of the cases, the number of non-causal variants are far more than causal variants. The auPRC will be affected by the imbalance ratio. Other metrics, such as auROC, precision, recall, f1-score are also necessary to provide more comprehensive benchmark results.

Thank you for your suggestion. It is correct that in all our datasets the proportion of causal variants of interest is low. A recent review of recommendations for applying ML in genomics recommends the AUPRC in such cases (Whalen et al., 2022). We do not include precision, recall and F1-scores since these require setting specific thresholds while zero-shot scores from different models do not have any obvious cutoffs. We have added AUROC results, however. Please see our updated text:

“While the AUPRC is generally recommended for imbalanced datasets where we are mostly interested in the positive minority class (Whalen et al., 2022), we also report the area under the receiver operating characteristic (AUROC). The main difference we see is a slight relative improvement of Enformer and Borzoi zero-shot scores (Figure A.6).”

Thank you for your thorough review and insightful suggestions. We have revised our manuscript and uploaded it above.

The threshold for minor allele frequency (MAF) in pathogenic OMIM variants, set at 0.1%, appears arbitrary, particularly since a 1% threshold is typically standard for rare variants. Conducting ablation studies would clarify the impact of this threshold on model performance and help justify this design choice.

Were ablation studies conducted to assess the impact of the 0.1% MAF threshold for OMIM variants on model performance? A comparison with a 1% threshold would align with conventional standards and strengthen the study’s design choices.

Thank you for your suggestion. We have added the suggested ablation, which yields very similar results. We note that the 0.1% cutoff is often used in the literature, e.g. in a rare disease screening application (Adhikari et al., 2020). Please see our updated text below:

“When using a more relaxed MAF cutoff of 1%, only 19 additional positive variants are included, resulting in very similar results (Figure A.1)“

The control matching process would also benefit from closest gene-specific stratification. By matching positive and negative samples more closely at the closest gene level, gene-specific biases could be reduced, making the control set construction more rigorous.

Thank you for your suggestion. We have included additional versions of the dataset matched by gene. This analysis is certainly of biological interest but requires dropping many causal variants that cannot be properly matched, especially for Mendelian traits. Please see our updated text:

“Also, we explored matching negatives from the same gene rather than from the same chromosome, which required dropping many variants that could not be properly matched, but with similar overall conclusions (Figure A.2).“

“When matching negatives from the same gene rather than the same chromosome, performance is lower overall, but Borzoi obtains a small advantage (Figure A.5).”

Finally, the methodology of matching controls based on Euclidean distance across TSS distance, MAF, and LD score assumes these features are equally relevant, which may not hold true. Given that MAF and LD score are more population-dependent, while TSS distance is more functionally related, this approach risks introducing biases in the control set.

Thank you for your suggestion. We acknowledge that matching multiple features always involves a tradeoff. We empirically find that our approach, which includes robust scaling, results in the matched features (TSS distance, MAF, LD score) themselves having very low predictive power (Table A.4), as intended.

The leave-one-chromosome-out (LOCO) validation approach in the linear probing evaluations, while aiming to test model generalization, has notable limitations. Chromosomes differ significantly in size and gene density, leading to inconsistencies in split sizes and class balance for variant groups. Smaller chromosomes, such as chromosome 21, have fewer variants, which can skew model performance and impact the validity of linear probing. Furthermore, features like minor allele frequency (MAF) and linkage disequilibrium (LD) score vary by chromosome and could add further bias if not carefully accounted for during these experiments.

Class proportion is maintained since we match each positive with a fixed number of negatives from the same chromosome. To account for differences in sample sizes between different chromosomes, we perform a weighted average of AUPRC. MAF and LD score naturally vary, but there is no reason to believe that they would vary in a biased manner towards a specific chromosome (fine-mapping results already exclude sex chromosomes and the MHC locus). We also note that a very similar LOCO approach has been adopted by expert statistical geneticists in a recent preprint (Fabiha et al., 2024), which appeared while our ICLR submission was under review.

Thank you for your helpful comments and suggestions. We have uploaded a revised version of our manuscript.

Feature extraction across model classes could be enhanced. For example, incorporating the L2 distance and cosine distance between reference and alternate embeddings (rather than only output features) in both supervised and self-supervised models could provide additional insights into model performance.

Please expand the feature set for supervised and self-supervised models to include L2 distances or cosine similarity metrics for embeddings

We now compare different zero-shot scoring approaches for gLMs. Please see our updated text:

“Good results have also been obtained comparing the embeddings of the alternate and reference alleles (DallaTorre et al., 2023; Mendoza-Revilla et al., 2024). We evaluate these different scoring approaches for each model (Table A.8) and choose the best performing one when benchmarking against other models (Table A.9).”

The benchmark would benefit from including variant classes that directly impact gene expression.

Please integrate gene expression variants from fine-mapped GTEx data (SuSiE) to enhance the benchmark’s relevance for expression-related causal variants

The focus of our paper is on variants affecting high-level traits, which are known to have limited overlap with eQTLs (Mostafavi et al., 2023). We nevertheless have analyzed the overlap between complex trait causal variants and fine-mapped eQTL variants, and included our findings in the revised manuscript. Please see our updated text:

“eQTL colocalization. We found that 103 putative causal variants for complex traits (9%) overlap with fine-mapped GTEx eQTL variants (Lonsdale et al., 2013; Wang et al., 2021); we found no such overlap for Mendelian trait variants, as expected given their low allele frequencies. The low overlap of complex trait and eQTL variants is well known and Mostafavi et al. (2023) discuss several hypotheses for the cause. We found that eQTL-overlapping variants are much easier to predict than non-eQTL-overlapping variants, across all model types (Figure 8D). We also note that Borzoi achieves a wide margin compared to other models and little is gained from ensembling. We observed that eQTL-overlapping variants are enriched in exonic variants (Fisher’s exact p = 8 × 10^−8) and, among non-exonic variants, they have lower TSS distances (Mann Whitney p = 4 × 10^−4), all of which could explain their increased predictability.”

We appreciate your thoughtful review and suggestions. We have revised our manuscript and uploaded it above.

The technical contribution may be considered low with minimal innovation.

We would like to emphasize that our main contribution is the much-needed systematic benchmarking of both the latest functional-genomics-supervised and self-supervised models on the task of predicting causal variants in human genetics. Our results provide insights into the capabilities and limitations of different approaches, and we believe that this kind of benchmarking study will facilitate future technical work that will help to advance the field.

More discussion about gLM promoter and its difference from GPN-MSA and other self-supervised trained models would be helpful.

We have now expanded the description of the model. Please see our updated text:

“gLM-Promoter was trained on 512-bp sequences centered at TSSs of protein-coding genes from reference genomes of animal species. TSS coordinates were obtained from the gene annotations available at NCBI Datasets (O’Leary et al., 2024). Species available at NCBI Datasets were subsampled, among those with gene annotations, to keep at most one per family. This resulted in 434 reference genomes. gLM-Promoter’s training objective follows GPN: base-pair-level tokenization and masked language modeling of local windows of 512-bp with downweighting of repeat positions (soft-masked in the reference genome). gLM-Promoter’s architecture follows ByteNet (Kalchbrenner et al., 2017; Yang et al., 2024), consisting of blocks alternating dilated convolutions and feedforward layers. Hyperparameters are displayed in Table A.6. Training took approximately 2 weeks using 4 NVIDIA A100 40GB GPUs.”

We now also analyze its performance across variant consequences. Please see our updated text:

“We also inspected the performance of gLM-Promoter on different consequences, given that it was trained only on promoters (Figure A.8). gLM-Promoter’s zero-shot scores perform better on proximal non-exonic and 5’ UTR variants, which lie in the regions of the gene with the highest overlap with the model’s training data (512 bp around the TSS). Except for the aforementioned classes in Mendelian traits, linear probing outperforms zero-shot scores.”

A cross comparison between the curated genetic variants with those from ClinVar would strengthen this study.

We now include an investigation into expert-curated ClinVar variants, finding that they lack coverage of the regulatory genome. Please see our updated text:

“Furthermore, expert-reviewed pathogenic variants in ClinVar are highly skewed towards coding and splice region variants, containing only a single promoter variant and no intergenic variants (Table A.7).”

Why choosing to use common variants as controls in the study of Mendelian traits? The causal variants of Mendelian traits are typically rare (i.e., low MAF).

That is correct. The reason why we use MAF > 5% variants as controls is that this is a recommended standalone criterion (BA1) to classify a variant as benign for Mendelian disorders [1]. This negative set is common in the literature, see e.g. CADD [2]. It would indeed be interesting to use benign rare variants as controls, but, in practice, we are not aware of any robust standard approach to confidently call rare benign variants.

References:
[1] Richards, Sue, et al. "Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology." Genetics in Medicine 17.5 (2015): 405-423.
[2] Rentzsch, Philipp, et al. "CADD-Splice—improving genome-wide variant effect prediction using deep learning-derived splice scores." Genome Medicine 13 (2021): 1-12.