Multi-trait genome-wide association studies (GWAS) use multi-variate statistical methods to identify associations between genetic variants and multiple correlated traits simultaneously, and have higher statistical power than independent univariate analysis of traits. Reverse regression, where genotypes of genetic variants are regressed on multiple traits simultaneously, has emerged as a promising approach to perform multi-trait GWAS in high-dimensional settings where the number of traits exceeds the number of samples. We extended this approach and analyzed different machine learning methods (ridge regression, random forests and support vector machines) for reverse regression in multi-trait GWAS, using genotypes, gene expression data and groundtruth transcriptional regulatory networks from the DREAM5 SysGen Challenge and from a cross between two yeast strains to evaluate methods. We found that genotype prediction performance, in terms of root mean squared error (RMSE), allowed to distinguish between genomic regions with high and low transcriptional activity. Moreover, model feature coefficients correlated with the strength of association between variants and individual traits, and were predictive of true trans-eQTL target genes, with complementary findings across methods. 1 Scientific Background Genome-wide association studies (GWAS) aim to find statistical associations between genetic variants and traits of interest using data from a large number of individuals [1, 2]. When multiple correlated traits are studied simultaneously, joint, multi-trait approaches can be more advantageous than studying the traits individually, due to increased power from taking into account cross-trait covariances and reduced multiple-testing burden by performing a single test for association to a set of traits [3, 4]. The most commonly used multi-trait GWAS approaches are based on a multivariate analysis of variance (MANOVA) or canonical correlation analysis (CCA) [4]. However, these are applicable only to cases where the number of traits is small. When analysing the effects of genetic variants on molecular traits (gene or protein expression levels, metabolite concentrations) or imaging features, we have to deal with a large number, often an order of magnitude or more greater than the sample size, of correlated traits simultaneously. For such cases, the standard procedure is still to conduct univariate linear regression or ANOVA tests for each genetic variant against each trait separately. While efficient algorithms exist to undertake this task, [5, 6, 7] the massive multiple-testing problem results in a significant loss of statistical power. An alternative approach to multi-trait GWAS has been to reverse the functional relation between genotypes and traits, and fit a multivariate regression model that predicts genotypes from multiple traits simultaneously, instead of the usual approach to regress traits on genotypes. The first study to do this explicitly used logistic regression and showed a significant increase in power compared to univariate methods [8]. Although the method as presented in [8] is still only valid when the number of traits is small, extending multivariate regression methods to high-dimensional settings is straightforward. Thus a recent study used L2-regularized linear regression of single nucleotide polymorphisms (SNPs) on gene expression traits to identify trans-acting expression quantitative trait loci (trans-eQTLs), and showed that this approach aggregates evidence from many small trans-effects while being unaffected by strong expression correlations [9]. ∗Corresponding author ar X iv :2 11 1. 00 10 8v 1 [ qbi o. G N ] 2 9 O ct 2 02 1 Despite these advances, several limitations and open questions remain unanswered in high-dimensional GWAS. First, linear models effectively search for the linear combination of traits that is most strongly associated to the genetic variant, but there is no a priori biological reason why only linear combinations should be considered. Second, while regularization allows to deal with high-dimensional traits, it does not address the problem of variable selection. For instance, in the case of gene expression, we expect that trans-eQTLs are potentially associated with many, but not all genes. Indeed, in [9] a secondary set of univariate tests is carried out to select genes associated to trans-eQTLs identified by the initial multi-variate regression. Finally, a systematic biological validation and comparison of the available methods is lacking. Here we address these questions by considering a wider range of machine learning methods (in particular, random forests regression (RFR) and support vector regression (SVR)) for reverse genotype prediction from gene expression traits. Hypothesizing that true trans-eQTL associations are mediated by transcription regulatory networks, we use simulated data from the DREAM5 Systems Genetics Challenge, and real data from 1,012 segregants of a cross between two budding yeast strains [10] together with the YEASTRACT database of known transcriptional interactions [11], to validate and compare these methods against L2-regularized linear regression (ridge regression (RR)). 2 Materials and Methods 2.1 Reverse genotype prediction For genotype prediction using machine learning models, the expression values were treated as explanatory variables whereas the genotype value of a variant was treated as a response variable. The prediction performance was measured by computing the root mean squared error (RMSE) between the predicted and the actual genotype value of variants. 2.2 Trans-eQTL target prediction Trans-eQTL target prediction was done using weights assigned to the features by the machine learning methods. We computed the area under the receiver operating characteristic (AUROC) curve to measure prediction performance. The weights (feature importance in case of RFR, and coefficients for RR and SVR) were compared against the true targets in the ground truth for each variant.