OmiGA for ultra-efficient molecular quantitative trait loci mapping

Domain: neuroscience
Published: 2026-02-12
Editorial outcome: minor_revision

Abstract
Molecular quantitative trait loci (molQTL) mapping is one of the most popular approaches to systematically characterize functional impacts of genomic variants, leading to advanced understanding of the regulatory mechanisms underpinning complex traits and diseases. However, when applied to high-throughput molecular phenotypes, the existing molQTL mapping tools often implement simple linear models, overlooking complex inter-individual relatedness, leading to false positives and insufficient statistical power. Here, we introduce OmiGA, an ultra-efficient omics genetic analysis toolkit, for molQTL mapping based on linear mixed model in populations with complex relatedness. Both computational simulations and real data analyses demonstrate that OmiGA outperforms the existing popular tools regarding molQTL discovery power, fine mapping of causal variants, colocalization of molQTL and trait associations, and computational efficiency. In summary, we recommend OmiGA for molQTL mapping in populations with complex relatedness, for example, those in the Farm animal Genotype-Tissue Expression project and family-based molQTL studies in humans.

Introduction
The majority of trait-associated variants discovered by genome-wide association studies (GWAS) are non-coding and act through regulating intermediate molecular phenotypes (for example, gene expression and alternative splicing)1,2. Molecular quantitative trait loci (molQTL) mapping has been proposed to be a powerful approach to characterize the regulatory effects of genomic variants in natural populations, a step forward to deciphering the regulatory circuitry of complex traits and diseases3–7. The conventional cis-molQTL mapping often involves two-step analytical framework: 1) is to perform a nominal association test and obtain the marginal effects of genetic variants on molecular phenotypes (e.g., gene expression) within the cis-region (e.g., 1Mb around transcription start site, TSS) of the target gene; 2) is to conduct multiple testing correction to identify eGenes, i.e., genes with significant cis-molQTL6. When analyzing these high-throughput molecular phenotypes (e.g., millions of gene expression across tissue and cell types), most existing molQTL mapping tools (e.g., FastQTL, tensorQTL, MatrixEQTL) implement simple linear regression models (LM) to reduce computational burdens8–10, resulting in the incomplete control of inter-individual relatedness, particularly in populations with complex relatedness (Fig. 1a).

Relatedness and population structure are crucial confounding variables that cannot be neglected11,12. To mitigate the impacts of relatedness among samples, a straightforward strategy is to analyze unrelated individuals only. This leads to a loss of statistical power of molQTL discovery, particularly in farm animals, where most individuals are related due to the recent intensive breeding and limited effective population size. In the conventional GWAS, the linear mixed model (LMM) has been proposed as a standard approach to account for the complex relatedness among individuals13,14. However, such GWAS tools are not designed for high-throughput molecular phenotypes. This has limited the widespread application of LMM, as estimating variance components for millions of molecular phenotypes using the restricted maximum likelihood (REML) algorithm entails intensive computation demands15,16. Furthermore, existing molQTL mapping tools often ignore non-additive genetic regulatory effects, even though they are prevalent for molecular phenotypes17,18. Identifying variants with non-additive effects will allow us to better understand the underlying biology of complex traits19, particularly in farm animals17.

In this study, we developed an ultra-efficient and powerful Omics Genetic Analysis toolkit (OmiGA) based on the LMM for molecular QTL mapping, which supports GPU acceleration. By conducting extensive simulations, we demonstrated that OmiGA outperformed LM-based tensorQTL regarding the statistical discovery power of molQTL and computational efficiency, especially for its robustness in controlling false positives in the presence of population stratification and complex relatedness. More specifically, in cis-eQTL mapping, OmiGA outperformed the commonly used tensorQTL in both CPU and GPU version. It was also >26 and >13 times faster than the existing tools (GCTA20 and LDAK21) in estimating the additive and dominant heritability of molecular traits, respectively. To showcase the utility of OmiGA, we mapped the cis-eQTL within four actual datasets sourced from pigs and humans. The pig datasets represent the population with complex relatedness, while human datasets represent the population consisting of unrelated individuals. We re-mapped cis-eQTL for 34 tissues in the pilot phase of the PigGTEx project, leading to a 64.54% increase in eGene discovery. Therefore, we recommend OmiGA (https://omiga.bio) for molQTL mapping in populations with complex relatedness, such as studies in the Farm animal Genotype-Tissue Expression (FarmGTEx) project22.

Results
Overview of OmiGA
OmiGA is a toolkit designed for ultra-efficient molQTL mapping analysis, while accounting for population structure and relatedness, with technical details in Methods. The toolkit has two key modules (Fig. 1b, c): (1) LMM-based methods for molQTL mapping (e.g., cis-QTL, context-dependent cis-QTL, and trans-QTL) with additive and/or non-additive effects; and (2) additive and/or non-additive (cis, trans, and global) heritability partitioning for high-throughput molecular phenotypes. OmiGA requires three fundamental input data components: (1) complete or imputed genotype data; (2) a numerical matrix of molecular phenotypes with genomic annotations (e.g., TSS); and (3) covariates accounting for the batch effects and other hidden confounders (e.g., cell type composition in bulk tissue samples).

OmiGA performs molQTL analysis which comprises three steps: (1) The data rotation step, which involves constructing the genetic relationship matrix (GRM) using genotype data, performing once the GRM factorization using Singular Value Decomposition (SVD) (with computational complexity of O(n3)), and rotating the genotype, phenotype, and covariate using the factorization results with O((m+p+c)n2), where n represent the sample size, p is the number of molecular traits, c is the number of covariates, and m is the number of variants; (2) The nominal association test step, which involves estimating molQTL effects, \hat{\beta} and V(\hat{\beta}), using a REML-free accelerated iterative dispersion that based on the newly developed AIDUL algorithm (see Methods)22. The computational complexities of AIDUL algorithm are lower, at O(tpnc2+mpc2), compared to its original version (i.e., iterative dispersion update to fit LMM, IDUL23), which is O(tmpnc2), with the iteration number denoted as t. Additionally, it computes the association P-value for each variant across all molecular phenotypes using Wald test statistics F = \frac{\hat{\beta}_{snp}^2}{V(\hat{\beta}_{snp})} \sim F(1, n-c-2);
(3) The multiple-testing correction step, which involves two optional strategies that includes permutation-free aggregated Cauchy association test (ACAT)24 and P-value-free FDR control method, implemented in ClipperQTL25. Compared to thousands of permutations in tensorQTL, ClipperQTL requires fewer permutations to achieve comparable statistical power of molQTL discovery. The permutation test step of ClipperQTL involves shuffling phenotypes using two optional strategies to serve as response variables. This is followed by estimating approximate partial correlation coefficients (APCC) between the permuted phenotypes and genotypes, which assumes that the relatedness contribution is negligible for the permuted phenotypes (as detailed in Methods).

To detect molQTL with dominant effects, we have implemented four optional statistical models (Methods). The computational procedures are similar to those used for additive molQTL mapping above, except when employing models that concurrently incorporate both additive and dominant GRMs (i.e., LMM(d+A+D)). These models estimate additive genetic variance \sigma^2_a, dominant genetic variance \sigma^2_d and residual variance \sigma^2_e for each phenotype, necessitating iterative computation of the inverse of the variance-covariance matrix (V-1) tens of thousands of times (computational complexity of O(tpn3)). To address this computational challenge, we introduce the Pre-Allocated Eigen Decomposition (PAED) method to estimate the relevant variance components under the null model before conducting the nominal association test. This method reduces the computational complexity to O(\kappa n3+tpn2), where \kappa is a finite value (i.e., 1000) dependent on the computational precision. PAED constructs a ‘prior’ combined GRM (cGRM) based on the variance component estimates and pre-computes the Eigen Decomposition, avoiding direct computing V-1 during the REML iteration. Subsequently, we optimize the statistical model by replacing the two GRMs with the cGRM, simplifying the model structure and enhancing computational efficiency. Ultimately, we leverage the AIDUL algorithm to efficiently calculate nominal association statistics.

We developed a specialized heritability estimation module for high-throughput molecular phenotypes within OmiGA. Unlike existing tools not optimized for molecular phenotypes, which repetitively estimate variance components for each phenotype, necessitating O(tpn3) computations for V-1 in every REML iteration, our approach is more efficient. We first perform Eigen Decomposition of the GRM once, obtaining V-1 through matrix-vector multiplication for each phenotype, resulting in a total computational complexity of O(n3+tpn2). For datasets with large sample size (e.g., >1000), we estimate variance components using two weighted least square regressions employing the IDUL algorithm. In the context of heritability models such as LMM(A_c+A_t) and LMM(A_c+A_t) with one GRM constructed from variants in the cis-region, we introduce a cis low-rank approximation (cisLRA) method. This reduces the V-1 computation complexity from O(n3) to O(k3) in the REML iteration, where k denotes the rank of the approximate low-rank matrix and k \ll n. Furthermore, for heritability estimation involving both global additive and dominant GRMs (LMM(A_g+D_g)), we leverage the PAED method to enhance computational speed and efficiency.

Runtime and memory requirements
We evaluated the computational performance of OmiGA in comparison with existing tools that are specifically designed for molQTL analysis (i.e., LMM-based APEX and LM-based tensorQTL) in three datasets, including MAGE (Multi-ancestry Analysis of Gene Expression)26, GEUVADIS (Genetic European Variation in Disease)27 and GIAD042 (300 pigs)28. MAGE comprises 731 individuals from five distinct populations that encompasses a total of 19,539 expressed genes and 14,267,797 variants. GEUVADIS comprises 462 individuals from five distinct populations that encompasses a total of 17,703 expressed genes and 10,313,747 variants. GIAD042 comprises 300 individuals from three distinct populations that encompasses a total of 14,008 expressed genes and 15,495,927 variants (see Methods for details). We conducted comparisons of these tools in the same computing environment, which was allocated 120 GB of memory and 10 CPU threads (Intel Xeon CPU Max 9462) within a single server node. We then performed cis-eQTL mapping using OmiGA, APEX, and tensorQTL (CPU version). Since APEX does not support direct GRM calculation, we pre-computed the GRM using GCTA v1.94.1 and converted it to the sparse matrix format file that APEX accepts. Both OmiGA and tensorQTL utilized the same input files. We repeated five times to compare their average elapsed time and maximum memory usage.

We observed that OmiGA was up to ~3.6 times faster than tensorQTL while requiring similar memory consumption (Table 1). To ensure that APEX could run properly, we set GRM elements to 0 for individuals with an estimated genetic relationship of more than 4th degree (i.e., genetic relatedness < 0.044), as detailed in the previous study29. Although the runtime of OmiGA was 0.21h for GEUVADIS, which was slightly faster (~1.8 times) than that of APEX based on a sparse GRM, OmiGA used much fewer computational resources than APEX (Table 1). More specifically, OmiGA required ~14 GB of memory to complete the whole analysis, only with ~16% of the maximum memory usage of APEX (~90 GB). Note that APEX cannot handle the MAGE and GIAD042 dataset with >14 million variants, despite the allocation of 500 GB of memory. To evaluate the scalability of OmiGA, we simulated a large-scale dataset containing 30,000 individuals with 19,539 tesed genes and 10 million variants. Our results demonstrate that OmiGA exhibits strong scalability: 1) When using 10 CPU threads (Intel Xeon CPU Max 9462), the LMM-based implementation of OmiGA outperformed the LM-based tensorQTL in terms of runtime; 2) When processing a dataset of the same scaled dataset under the same statistical model (LM), the CPU version of OmiGA is faster than GPU version of tensorQTL in completing the entire analyses (Supplementary Fig. 1).

Regarding heritability estimation, we compared OmiGA with two commonly used tools, namely GCTA v1.94.1 and LDAK v6. Because the GCTA and LDAK cannot automatically estimate heritability for tens of thousands of molecular phenotypes, we develop in-house pipelines. For global additive heritability estimation (i.e., Ag), we used the global additive GRM constructed by GCTA with --make-grm option as input for all three software. OmiGA showed >26 and >434 times faster than GCTA and LDAK, respectively (Table 2). For dominant heritability estimation (i.e., Ag+Dg), we used the dominant GRM constructed by GCTA with --make-grm-d option as input. OmiGA showed >13 and >204 times faster than GCTA and LDAK, respectively (Table 2). OmiGA enables efficient estimation of trans- and cis-heritability (i.e., At+Ac), whereas GCTA and LDAK do not yet provide a convenient way to achieve this (Table 2). In addition, we observed that the estimates of variance components from OmiGA were nearly identical to those obtained from GCTA and LDAK, except for the estimation of dominant heritability using LDAK (Supplementary Figs. 2–4).

OmiGA shows powerful performance in the simulations
Given the comparison results of resource requirements for different software above, we only considered OmiGA and tensorQTL, representing LMM and LM methods, respectively, to compare their statistical performances via extensive computational simulations. We then calculated the genomic inflation factor (\lambda), false discovery rate (FDR), F1-score, and true positive rate (TPR) (Methods) to evaluate their performances. One of the main aims of the simulations was to explore the influences of family relatedness on molQTL mapping. To mimic the effect of population stratification and relatedness among individuals, we conducted simulations in two scenarios i.e., simulations with close relatedness individuals (SIM-CLOSE) and with unrelated individuals (SIM-UNREL) (see Methods). SIM-CLOSE used the GIAD042 dataset with close family relatedness, containing 300 pigs from three different populations (Fig. 2a). SIM-UNREL used the MAGE dataset containing 731 unrelated individuals from five different sub-populations (Fig. 2b). As expected, there were higher degrees of relatedness among individuals in GIAD042 than in MAGE (Fig. 2a, b). Additionally, GIAD042 exhibited a much stronger inter-chromosome LD than MAGE (Fig. 2c, d, Supplementary Figs. 5–7).

We simulated gene expression levels with causal variants randomly sampled from cis-regions (designated as eGene) or no causal variants (designated as non-eGene), under various genetic architectures. We quantified \lambda by performing genome-wide eQTL mapping, where \lambda was defined as the median chi-squared statistic divided by its expected value at the null simulations. We defined the FDR of cis-eQTL discovery as the proportion of non-eGenes with multiple-testing corrected P values less than a given threshold (e.g., 0.05). We defined the statistical discovery power of eQTL as the proportion of true positives (TPR) under a specified false discovery rate (FDR) (e.g., 5%). We mimicked the effect of population stratification by generating a mean phenotype difference between three pig populations and mimicked the effect of relatedness by specifying polygenic background effects at different levels (Methods).

We first compared the performance of OmiGA and tensorQTL to account for the population structure and relatedness. As expected, simulation results showed that the test statistics from LM were inflated across all scenarios with relatedness of different degrees (Fig. 2e, f, Supplementary Fig. 8a–d). The genomic inflation factors in the test statistics of LM becomes smaller when increasing the number of genotype PCs in the model. The inflation in the test statistics for SIM-CLOSE population was stronger than that for SIM-UNREL population. This was due to the stronger long-range inter-chromosome LD in farm animals than that in humans (Fig. 2c, d). By contrast, there was almost no inflation in the test statistics of LMM with different numbers of genotype PCs across all simulations (Fig. 2e, f, Supplementary Fig. 8c, d), demonstrating the robustness of LMM-based OmiGA in accounting for long-range LD due to complex relatedness and small effective population size.

We further assessed the statistical power of OmiGA for cis-eQTL discovery, compared to tensorQTL. Based on the same number of genotype PCs, we quantified the FDR at a given significance level obtained from OmiGA and tensorQTL. We detected eGenes from OmiGA using ACAT pvalue<0.05 at the gene level and those from tensorQTL using Storey’s qvalue<0.05, which is used in the human GTEx3,30. In the SIM-CLOSE population, OmiGA and tensorQTL exhibited mean FDRs of 4.3% and 19.7%, respectively, and mean TPRs of 73.6% and 71.3% under 5% FDR (Fig. 3a). By contrast, in the SIM-UNREL population, OmiGA exhibited a slight improvement than tensorQTL, with the mean FDR of 3.7% and 6.0%, respectively (Fig. 3b). For genes exhibiting varying degrees of cis-heritability, OmiGA controlled the FDR more effectively and efficiently than tensorQTL regarding eGene discovery, particularly in the SIM-CLOSE population (Fig. 3c, Supplementary Figs. 8e, f, 9). The FDR of OmiGA could be controlled well using various multiple-testing approaches in both populations (Supplementary Fig. 10). tensorQTL, using the approach of beta-approximate with Storey’s qvalue, could not control the FDR well in the SIM-CLOSE population (Supplementary Fig. 11). Fitting more genotype PCs in tensorQTL reduces the FDR (Supplementary Fig. 12), which was in line with the observation regarding the genomic inflation factor. However, the accuracy of effect estimates from tensorQTL were lower than those from OmiGA (Fig. 3d, e). Additionally, the mean Chi-square statistics of non-causal variants in tensorQTL were significantly higher than those in OmiGA (Fig. 3f).

To examine the characteristics of eGenes identified exclusively by OmiGA in the SIM-CLOSE population, we conducted a comparative analysis of the eGenes identified by both tools. To ensure a similar FDR between OmiGA and tensorQTL, we compared the results from tensorQTL (with 10 genotype PCs) under the Benjamini–Hochberg method and those from OmiGA (with 2 genotype PCs) based on ACAT pvalue<0.05 at the gene level. The majority of eGenes were shared between OmiGA and tensorQTL (i.e., ‘shared’), evidenced by that mean TPR of shared eGenes identified by both tools was 97.8% across 50 simulations (Fig. 3g). However, there was still a certain proportion of true eGenes that were not discovered by both tools (Fig. 3h). The average percentage of true positives among eGenes detected by OmiGA only (i.e., ‘OmiGA_only’) was notably higher than those detected by tensorQTL only (i.e., ‘tensorQTL_only’) (Fig. 3h). Furthermore, the OmiGA_only eGenes exhibited a higher cis-heritability and lower trans-heritability compared to tensorQTL_only eGenes (Fig. 3i). These results indicate that OmiGA effectively controlled for long-range confounding effects originating from trans-region and identified a greater number of true positives compared to tensorQTL.

In terms of dominant cis-eQTL discovery, consistent with the observations for additive effects, we found that LMMs exhibit a higher TPR than LM (Fig. 3j). The dominant model incorporating both additive and dominant polygenic effects (i.e., LMM(d+A+D)) shows greater robustness than other models across various simulations scenarios (Fig. 3j, Supplementary Figs. 10, 13). Additionally, we assessed the robustness of OmiGA in detecting context-dependent cis-eQTL effects (also named interaction eQTL, ieQTL) (Supplementary Figs. 14–16). We simulated cis-eQTL with interaction effects by incorporating an interaction term, as described in the Methods. We observed that OmiGA identified a higher number of true ieGenes, when maintaining the same FDR as those using tensorQTL in the SIM-CLOSE population (Supplementary Fig. 15), with mean TPRs of 69.0% and 64.7% for OmiGA and tensorQTL, respectively. Additionally, OmiGA demonstrated an advantage over tensorQTL in the discovery of ieQTLs with small interaction effects (explained \le 0.1 of genetic variance) (Supplementary Fig. 16a). In the SIM-UNREL, there were no obvious differences in association test statistics between OmiGA and tensorQTL (Supplementary Fig. 16b).

Application of OmiGA to real data
We employed OmiGA to perform eQTL analysis in three tissues (duodenum, liver, and muscle) of 300 pigs from the GIAD042 population and lymphoblastoid cell line (LCL) of 731 humans from the MAGE population. As expected, there was inflation in the test statistics of LM-based tools, but not for LMM-based OmiGA (Fig. 4a, Supplementary Figs. 17, 18), in line with our observations from simulations (Fig. 2e, f). OmiGA, MLMA (with the option ‘--mlma-no-preadj-covar’), and GEMMA exhibit extremely high stability of genomic control compared to popular GWAS tools (Fig. 4a). OmiGA is the fastest one among these three tools. This indicates that OmiGA can effectively correct for the inflation caused by relatedness among individuals.

We compared the discovery power of cis-eQTL between OmiGA and tensorQTL (Fig. 4b, c). Based on the experience of the simulation study above, to maintain a comparable FDR between OmiGA and tensorQTL for eGene identification, we employed the Benjamini-Hochberg method and qvalues with lambda=0.05 to detect eGenes in tensorQTL for GIAD042 and MAGE, respectively. Across the three tissues in GIAD042, OmiGA outperformed tensorQTL by identifying 1083, 550, and 841 more eGenes detected in duodenum, muscle, and liver, respectively (Fig. 4b, Supplementary Fig. 19). Furthermore, OmiGA detected 270 additional eGenes in MAGE compared to tensorQTL (Fig. 4c). Consistent with the simulation findings, eGenes detected by OmiGA alone exhibited higher cis-heritability and lower trans-heritability compared to those detected by tensorQTL only (Fig. 4d, e, Supplementary Fig. 19), suggesting that OmiGA identified more true eGenes in these datasets. Moreover, OmiGA detected additional ieGenes in both datasets compared to tensorQTL. There are 512, 800, and 842 ieGenes that were only detected by OmiGA in the duodenum, muscle, and liver for GIAD042, respectively, and 133 ieGenes in LCL for MAGE (Supplementary Fig. 20). Building on these findings, we re-conducted cis-eQTL mapping for 34 tissues from the pilot phase of PigGTEx project5 using OmiGA. OmiGA (ACAT pvalue<0.05) achieved a 64.5% increase in eGene discovery compared to tensorQTL (BH-corrected pvalue<0.05), with an average of 1349 additional eGenes identified across 34 tissues (Fig. 4f). As anticipated, the pattern of cis- and trans-heritability of these eGenes discovered by OmiGA only was similar to those observed in simulations and the two real datasets above (Supplementary Fig. 19).

We employed SuSiE-Inf31 to perform fine mapping analysis for all eGenes. We derived the posterior inclusion probability (PIP) of the lead cis-eQTL of each eGene from OmiGA and tensorQTL. An average of 53% and 44% of tested eGenes from OmiGA and tensorQTL had credible sets of causal variants, respectively. Notably, the PIPs of lead cis-eQTL and independent cis-eQTL from OmiGA were significantly higher than those from tensorQTL (Fig. 4g, Supplementary Fig. 21), indicating that the cis-eQTL identified by OmiGA are more likely to be causal variants.

We then investigated whether the cis-eQTL discovered by OmiGA enhanced the statistical power for interpreting the GWAS associations of complex traits, compared to those identified by tensorQTL. We utilized GWAS summaries of 25 complex traits from Pig Biobank32 to conduct colocalization analysis with cis-eQTL from duodenum, liver, and muscle using Coloc v5.2.333. The posterior probability of colocalization (PP.H4) using cis-eQTL from OmiGA was significantly higher than those from tensorQTL (Fig. 5a, b, Supplementary Fig. 22). For instance, cis-eQTL of the TRIM52 gene in the duodenum identified by OmiGA colocalized (PP.H4=0.84) with the GWAS signal for lean meat percentage in pigs, whereas the PP.H4 for tensorQTL results was only 0.05 (i.e., no colocalization) (Fig. 5b). The colocalization signal around rs324104102 in the GWAS was present in the OmiGA results but absent from the tensorQTL results. The Pearson’s correlation coefficient between the P-values of GWAS and cis-eQTL from OmiGA is 0.89, which was higher than that between GWAS and cis-eQTL from tensorQTL (Pearson’s r=0.79) (Fig. 5c, d). This discrepancy arises because the test statistics for variants downstream of rs324104102 were inflation in the tensorQTL results, thereby reducing the correlation between the test statistics of GWAS and cis-eQTL. We further utilized OmiGA to conduct cis-eQTL mapping in lymphoblastoid cell lines of all 731 individuals from the MAGE population. The performance of OmiGA was comparable to that of tensorQTL. This was because the MAGE cohort comprises unrelated individuals, rendering the use of LMM less necessary.

We applied OmiGA to detect the dominant cis-eQTL in both MAGE and GIAD042 datasets and observed that most (~71%) eGenes identified by the additive model harbor dominant cis-eQTLs (Fig. 6a & Supplementary Fig. 23a–c). Compared to the additive model, the dominant model detected additional 340 eGenes with dominant cis-eQTLs in the MAGE dataset (Fig. 6a). The dominant lead cis-eQTLs were located farther from the TSS of their target genes than additive lead cis-eQTLs (Fig. 6b, Supplementary Fig. 23d–f). There into, we identified a significant dominant cis-eQTL associated with PKD1 in LCLs, but no significant additive cis-eQTL was identified (Fig. 6c, d). Multiple studies have reported that PKD1 is linked to autosomal dominant Polycystic Kidney Disease 134,35. While most eGenes were detectable by both additive and dominant models, we observed dominant-specific cis-eQTLs within these shared eGenes. For instance, SNP rs150261774, the dominant lead cis-eQTL of TNFRSF8, is independent of the additive lead cis-eQTL (Fig. 6e). These results demonstrate that the dominant model implemented in OmiGA can provide new insights into the non-additive genetic regulatory of gene expression.

Full summary statistics for LMM-based cis-eQTL of gene expression levels in the three tissues from the GIAD042 population, 34 tissues from the PigGTEx project, and LCLs from the MAGE population can be accessed publicly through the OmiGA portal at https://omiga.bio.

Methods
Inclusion and ethics
Not applicable because no biological samples were collected, and no animal handling was performed for this study.

Linear mixed model for standard molecular QTL mapping
To detect molQTL with additive effects, we consider the following standard linear mixed model:
LMM(a+A): y=X\alpha+s_a\beta_a+g_a+e (1)
where y is an n\times1 vector of standardized molecular phenotypes, X is an n\times c matrix of customized covariates including a column of 1, \alpha is the fixed effect corresponding to covariates, s_a is an n\times1 vector of mean centered genotype values in additive coding (0/1/2 for AA/Aa/aa) for the genetic variant being tested, \beta_a is the variant’s additive genetic effect, g_a is an n\times1 vector of total additive genetic effects with g_a \sim N(0, A\sigma^2_a) where A= \frac{MM^0}{m} is an additive genetic relationship matrix (GRM), of which M is a matrix of mean centered genotypes for genome-wide m genetic variants and p_i is the MAF at ith variant, e is n\times1 vector of residuals with e \sim N(0, I_n\sigma^2_e) where I_n is an n\times n identity matrix. To obtain association statistical magnitude, we employed exact and accelerated methods described below.

AIDUL for association test
To accelerate the association test process, we propose an accelerated iterative dispersion update to fit linear mixed model (AIDUL) by modifying the IDUL (iterative dispersion update to fit linear mixed model) algorithm23. The computational complexities for these two algorithms are O(tpnc2+mpc2) and O(tmpnc2), respectively, with the iteration number denoted as t. We first randomly select a small subset of variants to estimate the variance components (\tau^{-1} and \eta=\sigma^2_a/\sigma^2_e) using the IDUL algorithm, where \tau^{-1} is the variance of the residual errors (\sigma^2_e), \eta is the ratio between the additive genetic variance (\sigma^2_a) and \sigma^2_e, MVN denotes multivariate normal distribution. The eigendecomposition of A can be written as QTQ0 with QQ0=I_n, where jth column of Q is an eigenvector, whose corresponding eigenvalue is the jth diagonal element of the diagonal matrix T. Let W=(X, s_a) and b=(\alpha, \beta_a), we can write the formula (1) as follows:
y=Wb+g_a+e (2)
We rotate the formula (2) by multiplying Q0 to get
y_Q \sim MVN(W_Q b, \tau^{-1}H) (3)
where y_Q=Q0y, W_Q=Q0W, and H=\eta T+I_n. We perform IDUL until \eta converges by repeating:
Step 0: assign a constant to \eta_i (e.g., 0.1) for the first iteration (i=1) of IDUL.
Step 1: fit formula (3) using weighted least squares with weight H^{-1} to obtain the residual r.
Step 2: fit r^2 \sim MVN(\mu+diag(T)\gamma, \tau^{-1}H^2) using weighted least squares with weight H^{-2} to obtain estimates \hat{\gamma} and \hat{\mu}.
Step 3: update \eta_{i+1} = \frac{\hat{\gamma}}{\hat{\mu}} and repeat steps 1–3 until \frac{|\eta_{i+1}-\eta_i|}{1+\eta_{i+1}} < 10^{-6} or the number of iterations achieve to the maximum iterations (e.g., 100).
After the \eta converges, we perform AIDUL process by obtaining the average estimates of \eta for the subset of variants. We then let \eta=\eta_{avg} in formula (3) and obtain the effect estimates \hat{\beta}_a and its variance V(\hat{\beta}_a), where \hat{b} = (W_Q^0 W_Q)^{-1} W_Q^0 y_Q of which W_Q = W_Q H^{-1} and y_Q = y_Q H^{-1}, V(\hat{b}) = \sigma^2_e(n-c-1)^{-1} diag((W_Q^0 W_Q)^{-1}). Under the null hypothesis, the Wald test statistics F_{Wald} \sim F(1, n-c-1) and the P value of effect (\beta_a) of tested genetic variant can be calculated.

Multiple-testing corrections for cis-molQTL mapping
To adjust the multiple testing across molecular phenotypes in cis-molQTL mapping, we implement two approaches including ACAT24 and Clipper25.
ACAT is a permutation-free and computationally efficient P-value combination method. We compute the ACAT test statistic using a linear combination of transformed P-values by T_{ACAT} = \frac{1}{k} \sum_{i=1}^k tan((0.5-p_i)\pi), where k is the total number of variants for target gene, p_i is the nominal P-value of ith variants tested in association test, and tan((0.5-p)\pi) is Cauchy distributed if p is from the null distribution. Therefore, we then used T_{ACAT} to calculate the P-value of each molecular trait using the Cauchy cumulative distribution function to obtain the trait-level P-value.
ClipperQTL requires the phenotype-level P-value significance in cis-QTL. ClipperQTL requires the experimental measurements under the experimental and background conditions to determine the phenotype-level FDR. We define the phenotype-level partial correlation between phenotypes and lead variants observed in nominal association test as experimental measurements, and those from permutation test as background measurements for Clipper, where the lead variant is the most significant/correlated variant for a given phenotype or permuted phenotype.
To obtain the experimental measurements, we compute the partial correlation (r_p) between y and s_a for lead variant using \hat{\beta}_a and V(\hat{\beta}_a) estimated from nominal association test (i.e., LMM(a+A) model (1)). The partial correlation r_p = \frac{\hat{\beta}_a}{\sqrt{\hat{\beta}_a^2 + (n-c-1)V(\hat{\beta}_a)}}. To obtain the background measurements, we implemented two approaches for permutation test including direct permutation and residual permutation schemes.
Direct scheme. We perform direct permutation test using the following model:
y^e=X\alpha^e+s_a\beta^e_a+g^e_a+e^e (4)
where y^e is an n\times1 vector of molecular phenotypes being shuffling, \alpha^e, \beta^e_a, g^e_a \sim N(0, A\sigma^2_{ae}), and e^e \sim N(0, I_n\sigma^2_e) are fixed effect of covariates, genetic effect of tested variant, total additive genetic effects, and residuals, respectively, corresponding to permuted phenotype values.
To eliminate the computationally intensive in the permutation process, we assign zero to \sigma^2_{ae} for permuted phenotypes assuming the genetic contribution is negligible. Therefore, the term g^e_a can be neglected and the approximate model for direct permutation is a simple linear regression model:
y^e=X\alpha^e+s_a\beta^e_a+e^e (5)
To eliminate the effects of covariates, we calculate the residual y^e_r = y^e - X(X0X)^{-1}X0y^e and s_{a,r} = s_a - X(X0X)^{-1}X0s_a. Then the partial correlation r_p = cor(y^e_r, s_{a,r}) = \frac{y^e_r, s_{a,r}}{\sqrt{y^e_r, y^e_r} \sqrt{s_{a,r}, s_{a,r}}}. The partial correlation can be simplified as y^e_r = y^e - X(X0X)^{-1}X0y^e assuming the effect of a single variant on permuted y^e is negligible. To compute the partial correlation, we standardize the y^e_r and s_{a,r} to have zero mean and unit sum of squares. Therefore, we simplify the calculation of the partial correlation as r_p = cor(y^e_r, s_{a,r}) = y^e_r, s_{a,r}.

Residual scheme. We perform residual permutation using ClipperQTL method25. Different from direct permutation scheme, residual permutation scheme first residualise phenotype's before shuffling step, and calculate r_p using y^e_r = shuffle(y-X(X0X)^{-1}X0y) as above. The residual scheme requires data sets with large sample sizes (>450 is recommended) and the direct scheme is suitable for a wider range of sample sizes.
After finishing permutation process, we can obtain a matrix P_r = [r_{p,1}, r_{p,2}, ..., r_{p,B}] consisting of the maximum APCC for corresponding conditions, where the first column is experimental measurements, the second and subsequent columns are background measurements, subscript P is the number of tested phenotypes, and superscripts B is the number of permutations. As recommend by Zhou et al.25, B is set at 1000 by default for the direct scheme and 20–100 for the residual scheme.
The variant-level significance threshold for a target molecular trait was defined as the 5th percentile of the null distribution of permuted P-values for the lead variant. This null distribution was generated through permutation analysis in direct permutation or residual permutation, as described above.

LMM for interaction molecular QTL mapping
Extending the standard LMM(a+A) model (1), we test interaction effects between target genetic variant and interaction factors (e.g., cell-types or ancestry component) on molecular phenotypes using the following model:
y=X\alpha+s_a\beta_a+(s_a \circ i)\beta_{ai}+g_a+e (6)
where \circ is Hadamard product, i is an n\times1 vector of the standardized interaction factor, (s_a \circ i) is the interaction term between genotype and interaction factor, \beta_{ai} is the interaction effect of the genetic variant being tested, rest of the term the same as in LMM(a+A) model (1). Under the null hypothesis, the Wald test statistics F_{Wald} = \frac{\hat{\beta}_{ai}^2}{V(\hat{\beta}_{ai})} \sim F(1, n-c-2) and the exact P value of effect (\beta_{ai}) of tested genetic variant can be calculated.
Permutation test for interaction cis-molQTL mapping
Direct the permutation test for interaction cis-molQTL mapping using the model (5), we follow the method for standard cis-molQTL mapping and define the interaction variant as the most significant \beta_a and/or \beta_{ai}. To eliminate the effects of covariates, we calculate the partial correlation (r_p) between y and (s_a \circ i) using the same methods as above. To avoid the phenotypes residualizing process repeatedly for each variant, we extend the APCC used in the model (5) for interaction analysis, that simplify it as y^e_r = y^e - X(X0X)^{-1}X0y^e assuming the effect of a single variant on permuted y^e is negligible. To compute the partial correlation, we standardize the y^e_r and (s_a \circ i)_r to have zero mean and unit sum of squares. Therefore, we simplify the calculation of the partial correlation as r_p = cor(y^e_r, (s_a \circ i)_r) = y^e_r, (s_a \circ i)_r.

LMM for dominant molecular QTL mapping
To identify molQTL with dominant effects, we implemented three dominant linear mixed model including LMM(d+A) model corrected with GRM(A), LMM(d+D) model corrected with dominant GRM(D), and LMM(d+A+D) model corrected with both additive and dominant GRMs (8, 9, 10):
LMM(d+A): y=X\alpha+s_d\beta_d+g_a+e (8)
LMM(d+D): y=X\alpha+s_d\beta_d+g_d+e (9)
LMM(d+A+D): y=X\alpha+s_d\beta_d+g_a+g_d+e (10)
where s_d is an n\times1 vector of standardized genotype values in dominant coding (0/1/0 for AA/Aa/aa) for the genetic variant being tested, \beta_d is the variant’s dominant genetic effect, g_d is an n\times1 vector of total dominant genetic effects with g_d \sim N(0, D\sigma^2_d) where D= \frac{SS0}{m} is a dominant GRM, of which S is a matrix of standardized genotypes for genome-wide m genetic variants and p_i is the MAF at ith variant, and the rest of terms the same as in (1).
To perform association test based on LMM(d+A) or LMM(d+D) models, we can efficiently calculate the \hat{\beta}_d and V(\hat{\beta}_d) using IDUL the same as formula (3). To perform association test based on LMM(d+A+D), we first estimate the additive (\sigma^2_a) and dominant (\sigma^2_d) genetic variances based on
y=X\alpha+g_a+g_d+e (11)
We obtain the estimates \hat{\sigma}^2_a, \hat{\sigma}^2_d, and \hat{\sigma}^2_e using AI-REML algorithm described in below. Then we can calculate the \hat{\beta}_d and V(\hat{\beta}_d) using IDUL the same as formula (3) on model (12) as
y=X\alpha+s_d\beta_d+g_c+e (12)
where g_c \sim N(0, K\sigma^2_c) is an n\times1 vector of total genetic effects combining additive and dominant effects where K=A\sigma^2_a+D\sigma^2_d is a combined GRM, and the corresponding P-value from the Wald test statistics F_{Wald} = \frac{\hat{\beta}_d^2}{V(\hat{\beta}_d)} \sim F(1, n-c-1). To adjust the multiple testing across molecular phenotypes in dominant cis-QTL mapping, we implement permutation test for \beta_d using the same methods as standard LMM(a+A) model (1) and define which phenotypes have significant dominant cis-QTL effects.

Global heritability
Additive global heritability. We implement the calculation of the additive global heritability (h^2_Ag) for molecular phenotypes using the following additive genetic model:
LMM(A_g): y=X\alpha+g_a+e (13)
where y is an n\times1 vector of standardized molecular phenotypes, X is an n\times c matrix of customized covariates including a column of 1, \alpha is the fixed effect corresponding to covariates, g_a is an n\times1 vector of total additive genetic effects with g_a \sim N(0, A\sigma^2_a) where A= \frac{MM^0}{m} is an additive GRM, of which M is a matrix of mean centered genotypes for genome-wide m genetic variants and p_i is the MAF at ith variant, e is n\times1 vector of residuals with e \sim N(0, I_n\sigma^2_e) where I_n is an n\times n identity matrix. If the sample size was sufficiently large (e.g., exceeding 1000), we use the IDUL algorithm to estimate the variance components \sigma^2_a and \sigma^2_e using the method described as follows.
The additive global heritability estimation model assumes y \sim MVN(X\alpha, \eta\tau^{-1}A), e \sim MVN(0, \tau^{-1}I_n), where \tau^{-1} is the variance of the residual errors (\sigma^2_e), \eta is the ratio between the \sigma^2_a and \sigma^2_e, MVN denotes multivariate normal distribution. The eigendecomposition of A can be written as QTQ0 with QQ0=I_n, where jth column of Q is an eigenvector, whose corresponding eigenvalue is the jth diagonal element of the diagonal matrix T. We rotate the formula (13) by multiplying Q0 to get
y_Q \sim MVN(X_Q \alpha, \tau^{-1}H) (14)
where y_Q=Q0y, X_Q=Q0X, and H=\eta T+I_n. We perform IDUL until \eta converges by repeating:
Step 0: assign a constant to \eta_i (e.g., 0.1) for the first iteration (i=1) of IDUL.
Step 1: fit formula (3) using weighted least squares with weight H^{-1} to obtain the residual r.
Step 2: fit r^2 \sim MVN(\mu+diag(T)\gamma, \tau^{-1}H^2) using weighted least squares with weight H^{-2} to obtain estimates \hat{\gamma} and \hat{\mu}.
Step 3: update \eta_{i+1} = \frac{\hat{\gamma}}{\hat{\mu}} and repeat steps 1–3 until \frac{|\eta_{i+1}-\eta_i|}{1+\eta_{i+1}} < 10^{-6} or the number of iterations achieve to the maximum iterations (e.g., 100).
After IDUL converges, we obtain the residual variance estimates \hat{\sigma}^2_e = \hat{\tau}^{-1}, the additive genetic variance \hat{\sigma}^2_a = \eta^*\hat{\tau}^{-1}, and the h^2_Ag = \hat{\sigma}^2_a / (\hat{\sigma}^2_a + \hat{\sigma}^2_e).

Dominant global heritability. We implement the calculation of the dominant global heritability (h^2_Dg) for molecular phenotypes using the following dominant genetic models:
LMM(D_g): y=X\alpha+g_d+e (15)
LMM(A_g+D_g): y=X\alpha+g_a+g_d+e (16)
where g_d is an n\times1 vector of total dominant genetic effects with g_d \sim N(0, D\sigma^2_d) where D= \frac{SS0}{m} is a dominant GRM. To estimate dominant genetic variance of LMM(D_g) model (15), we employed IDUL algorithm using method the same as for additive global heritability estimation (i.e., LMM(A_g) model (13)) and the h^2_Dg = \hat{\sigma}^2_d / (\hat{\sigma}^2_d + \hat{\sigma}^2_e). To estimate the \sigma^2_a and \sigma^2_d for LMM(A_g+D_g) model (16), we employed AI-REML algorithm optimized for molecular phenotypes. For LMM(A_g+D_g) model (16), the variance-covariance matrix of y is V = \sigma^2_a A + \sigma^2_d D + \sigma^2_e I. To reduce the repeated calculation across phenotypes for each iteration in AI-REML, we propose a pre-allocated Eigen Decomposition (PAED) method avoiding direct matrix inversion described as follows.
In PAED method, let P = \sigma^2_a A + \sigma^2_d D = \kappa(\omega_1 A + \omega_2 D) with \omega_1 + \omega_2 = 1 and \kappa is a scale factor, we can pre-calculate P with different ratios between \omega_1 and \omega_2 (e.g., \omega_1/\omega_2 \in \{1, 1.01, 1.02, ...\}). We pre-perform the Eigen Decomposition for all recalculated P above. In the iteration process of AI-REML, we can quickly extract the eigenvalues (E_{val}) and eigenvectors (E_{vec}) of P with matching \omega_1 and \omega_2 for calculating V^{-1}.
We can easily to calculate
V^{-1} = (\kappa P + \sigma^2_e I)^{-1} = E_{vec} diag(\kappa E_{val} + \sigma^2_e)^{-1} E_{vec}^0 with \kappa = (\sigma^2_a + \sigma^2_d) / (\omega_1 + \omega_2). Finally, we obtain the h^2_Dg = \hat{\sigma}^2_d / (\hat{\sigma}^2_a + \hat{\sigma}^2_d + \hat{\sigma}^2_e).

Cis and trans heritability
Additive cis heritability. For estimating the contribution of cis-region variants to molecular phenotypes, we implement
LMM(A_c): y=X\alpha+g_{a,cis}+e (17)
LMM(A_c+A_t): y=X\alpha+g_{a,trans}+g_{a,cis}+e (18)
LMM(A_g+A_c): y=X\alpha+g_{a,global}+g_{a,cis}+e (19)
where g_{a,cis} \sim N(0, A_c\sigma^2_{a,cis}), g_{a,trans} \sim N(0, A_t\sigma^2_{a,trans}), g_{a,global} \sim N(0, A_g\sigma^2_{a,global}) where A_c, A_t, and A_g are additive cis-GRM, trans-GRM, and global-GRM constructed using genotypes from cis-region variants, trans-region variants, and genome-wide variants, respectively. By default, we define the cis-region as the \pm 1Mb region surrounding the transcription start sites (TSS) of the target molecular phenotype, whereas the trans-region refers to chromosomes distinct from the one where the target molecular phenotype is situated.
For LMM(A_c) model (17), we estimate the residual variance \hat{\sigma}^2_e = \hat{\tau}^{-1}, the cis additive genetic variance \hat{\sigma}^2_{a,cis} = \eta^*\hat{\tau}^{-1} using IDUL the same as for LMM(A_g) model (13) and obtain the h^2_{Ac} = \hat{\sigma}^2_{a,cis} / (\hat{\sigma}^2_{a,cis} + \hat{\sigma}^2_e).
For LMM(A_c+A_t) model (18), we proposed a cis low-rank approximation (cisLRA) method to implement AI-REML optimized for molecular phenotypes. In cisLRA, we assume that the cis-GRM can be represented with low-rank matrices due to high linkage disequilibrium in cis-region. We decomposed cis-GRM into low-rank matrices multiplication using Partial Singular Value Decomposition (PSVD) A_c = U_k \Sigma_k \Lambda_k U_k^0, where k \ll n. A_c = U_k \Sigma_k \Lambda_k U_k^0, A_t = \sigma^2_t A_t + \sigma^2_e I, D = \sigma^2_d \Sigma, where A_c^{-1} can be quickly calculated using Eigen Decomposition. We calculate the variance-covariance matrix of y in the iteration process of the AI-REML as
V = \sigma^2_t A_t + \sigma^2_c A_c + \sigma^2_e I = \sigma^2_t A_t + \sigma^2_c U_k \Sigma_k \Lambda_k U_k^0 + \sigma^2_e I = \sigma^2_t A_t + U_k (\sigma^2_c \Sigma_k \Lambda_k + U_k^0 \sigma^2_e I) U_k^0.
where D^{-1} + \Lambda_k U_k^0 A_t^{-1} U_k is a k\times k low-rank matrix that can be inversed quickly. We obtain the h^2_{Ac} = \hat{\sigma}^2_{a,cis} / (\hat{\sigma}^2_{a,global} + \hat{\sigma}^2_{a,cis} + \hat{\sigma}^2_e).

Dominant cis heritability. For estimating the contribution of cis-region dominant genetic effects (h^2_{Dc}) to molecular phenotypes, we implement three genetic models described as follows:
LMM(D_c): y=X\alpha+g_{d,cis}+e (20)
LMM(A_t+D_c): y=X\alpha+g_{a,trans}+g_{d,cis}+e (21)
LMM(A_g+D_c): y=X\alpha+g_{a,global}+g_{d,cis}+e (22)
where g_{d,cis} \sim N(0, D_c\sigma^2_{d,cis}) where D_c is dominant cis-GRM constructed using genotypes from cis-region variants. We estimate the variance components for LMM(D_c) model (20), LMM(A_t+D_c) model (21), and LMM(A_g+D_c) (22) using the same methods as for LMM(A_c) model (17), LMM(A_c+A_t) model (18), and LMM(A_g+A_c) (19), respectively. We obtain the h^2_{Dc} = \hat{\sigma}^2_{d,cis} / (\hat{\sigma}^2_{d,cis} + \hat{\sigma}^2_e), the h^2_{Dc} = \hat{\sigma}^2_{d,cis} / (\hat{\sigma}^2_{a,trans} + \hat{\sigma}^2_{d,cis} + \hat{\sigma}^2_e), and the h^2_{Dc} = \hat{\sigma}^2_{d,cis} / (\hat{\sigma}^2_{a,global} + \hat{\sigma}^2_{d,cis} + \hat{\sigma}^2_e) for LMM(D_c) model (20), LMM(A_t+D_c) model (21), and LMM(A_g+D_c) (22), respectively.

Trans heritability. In addition to models (18) and (21), which include a trans genetic component, we implement the alternative models to calculate the additive and dominant trans heritability (h^2_{At} and h^2_{Dt}) for molecular phenotypes:
LMM(A_t): y=X\alpha+g_{a,trans}+e (23)
LMM(D_t): y=X\alpha+g_{d,trans}+e (24)
where g_{a,trans} \sim N(0, A_t\sigma^2_{a,trans}) and g_{d,trans} \sim N(0, D_t\sigma^2_{d,trans}) where A_t and D_t are the additive trans-GRM and dominant trans-GRM constructed using genotypes from trans-region variants of target molecular phenotype. We estimate the variance components of LMM(A_t) model (23) and LMM(D_t) model (24) using IDUL the same as for LMM(A_c) model (17). We obtain the h^2_{At} = \hat{\sigma}^2_{a,trans} / (\hat{\sigma}^2_{a,trans} + \hat{\sigma}^2_e) and h^2_{Dt} = \hat{\sigma}^2_{d,trans} / (\hat{\sigma}^2_{d,trans} + \hat{\sigma}^2_e) for LMM(A_t) model (23) and LMM(D_t) model (24), respectively.

Summary of real data
We collected three public datasets including GIAD04228, MAGE26, and GEUVADIS27 described as below.
GIAD042 comprises 300 pigs from three distinct sub populations: Duroc (n=100), Yorkshire (n=100), and Landrace (n=100). The Duroc (DUR) individuals are offspring from 33 sows and 10 boars, Yorkshire (YOR) individuals from 84 sows and 43 boars, and Landrace (LAN) individuals from 74 sows and 18 boars. The genotype data of all 300 individuals were obtained from the GigaScience GigaDB database (https://doi.org/10.5281/zenodo.10072081). Following quality control across all 300 pigs, we retained 15,495,927 biallelic variants with MAF \ge 5% and MCA \ge 6 in 18 autosomes. Additionally, we obtained 900 RNA-Seq samples from three tissues (duodenum, liver, and muscle). To exclude lowly expressed genes (i.e., TPM < 0.1 and/or raw read counts < 6 in more than 80% of samples within each of tissues), we retained a total of 14,007, 12,512, and 9164 expressed autosomal genes in duodenum, liver, and muscle, respectively. We normalized the raw counts by trimmed mean of M values (TMM) within each tissue and used the inverse normal transformed TMM for subsequent analyses.
MAGE (Multi-ancestry Analysis of Gene Expression) comprises 731 unrelated individuals from five distinct human populations including AFR (n=196), AMR (n=113), EAS (n=141), EUR (n=142), and SAS (n=139). We obtained the genotype data of all 731 individuals from the 1000 Genomes Project (1KGP) samples (http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_phased/). After quality control across all 300 individuals, we kept 14,267,797 biallelic variants with MAF \ge 1% and MAC \ge 6 in 22 autosomes. We obtained 731 RNA-Seq samples from EBV-transformed lymphoblastoid cells lines (LCLs). To eliminate the lowly expressed genes using the same rules as did in GIAD042, we obtained a total of 19,539 expressed autosomal genes and used their inverse normal transformed TMM for subsequent analyses.
GEUVADIS (Genetic European Variation in Disease) comprises 462 unrelated individuals from five distinct populations: CEU, FIN, GBR, TSI and YRI. The genotype data were sourced from the 1000 Genomes Project. After quality control across all 462 individuals, we retained 10,313,747 biallelic variants with MAF \ge 1% and MAC \ge 6 in 22 autosomes. We obtained 462 RNA-Seq samples from LCLs. To eliminate the lowly expressed genes using the same rules as did in MAGE, we obtained a total of 17,703 expressed autosomal genes and used their inverse normal transformed TMM for subsequent analyses.

Genome-wide linkage disequilibrium
To evaluate the pattern of genome-wide linkage disequilibrium, we split the genotype data into different sub-populations of GIAD042 and MAGE datasets. For each of all sub-populations, we kept SNPs with MAF > 5% for subsequent analysis. We computed the intra- and inter-chromosomal mean LD using GEAR56.

Genotype data used in simulations
To assess the performance of OmiGA in comparison to other methods, we conducted a comprehensive simulation study, utilizing both populations consisting of related individuals and populations comprising only unrelated individuals. We used the WGS data from two public populations including the GIAD04228 and MAGE26 datasets above mentioned, where the GIAD042 is a pig population with close relatedness and the MAGE is a human population without related individuals. We used PLINK v2.057 and BCFtools58 to carry out quality control on the datasets. For both datasets, we retained biallelic autosomal variants with call rate \ge 95%, minor allele frequency (MAF) \ge 5%, and minor allele count (MAC) \ge 6. We phased and imputed the missing genotypes using Beagle 5.246 without reference panel. We conducted LD pruning using the flag --indep-pairwise 50 50 0.25 in PLINK v2.057, resulting in 770,243 and 492,797 variants in the GIAD042 and the MAGE, respectively, for the downstream simulations. Given the significant difference of the degree of kinship among individuals, we take the simulations using the GIAD042 and the MAGE named as SIM-CLOSE and SIM-UNREL below, respectively.

Schemes for simulating gene expression
To evaluate the performance of OmiGA in cis-eQTL mapping, we simulated genes with causal variants within cis-region in SIM-CLOSE and SIM-UNREL. We defined 1Mb upstream and down stream of the transcription start site (TSS) of each gene as the cis-region. Based on the genome coordinate of genes from Ensembl, we used only those genes located on chromosomes 1 and 2 and have at least 30 cis-variants for the simulations.
Additive cis-eQTL simulation
Various study designs suggest that ~60–90% of the genetic variance in expression can be attributed to trans-acting variation59–62. However, trans-eQTLs are notoriously difficult to find in humans is mainly due to effect sizes are considerably small, making it more difficult to achieve statistical significance compared to cis-eQTLs63. Fisher’s polygenic model64 assumes that a trait is controlled by numerous loci, each of which exerting a weak effect. The framework can be utilized to simulate the causal effects within trans-regions. We can denote the contributions from cis (local) and trans (distant) regions to the gene expression (y) of an individual as:
y = \mu + \sum_{i=0}^{m_{cis}} c_i \beta_i + \sum_{j=0}^{m_{trans}} t_j \gamma_j + \epsilon (25)
where \mu is the overall mean, c and t represent the genotypes for causal variants in the cis and trans-regions, respectively, \beta_i and \gamma_j represent their corresponding effect size, m_{cis} and m_{trans} represent the number of causal cis-variants and trans-variants, respectively, and \epsilon is the residual error. In the trans-region, the cumulative effects of these causal variants can be approximated as a normally distributed random variable (polygenic_{trans}) in accordance with the Central Limit Theorem65.
Therefore, we simplify the formula (25) to simulate the expression levels for each gene as follows:
y = \mu + C\beta + polygenic_{trans} + \epsilon (26)
where y \in R^n is a vector representing the simulated gene expression levels across n individuals, C is the matrix the centered genotypes of m_{cis} causal variants in the cis-region of target gene, c_i = s_i - 2p_i, where s \in \{0,1,2\} and p represent the MAF of ith causal variants, the vector \beta represents the effect sizes of causal variants on gene expression, the polygenic_{trans} obey the multivariate normal distribution of N(0, \sigma^2_{trans} A_?) with A_? denotes the kinship matrix constructed using the LD-independent variants in trans-region and the \sigma^2_{trans} = \sigma^2_{genetic} (1 - h^2_{cis} - h^2_{trans}), where h^2_{cis} is the cis-heritability, and h^2_{trans} trans-heritability. The residual error \epsilon \sim N(0, \sigma^2_e I), with \sigma^2_e = \sigma^2_{genetic} (1 - h^2_{cis} - h^2_{trans}) and polygenic_{trans} = \sigma^2_{genetic} \times A_? with A_? denotes the kinship matrix constructed using the LD-independent variants in trans-region and scale it by multiplying the factor \sigma^2_{trans} = Var(polygenic_{trans}).
Based on the strategies above, we simulated 3556 eGenes and 1879 non-eGenes in the SIM-CLOSE, and 2171 eGenes and 1144 non-eGenes in the SIM-UNREL, where the eGenes and the non-eGenes are the genes simulated with at least one causal cis-variant and zero causal cis-variant, respectively. All the simulations were repeated 50 times.

Dominant cis-eQTL simulation
Similar to the additive cis-eQTL simulation, we simulated genes with causal dominant variants within the cis-region in SIM-CLOSE and SIM-UNREL. The expression levels for each gene are as follows:
y = \mu + C\gamma + polygenic_{trans} + \epsilon (27)
where \mu, polygenic_{trans}, and \epsilon is the same as those in formula (26), C is the matrix the standardized genotypes of m_{cis} causal variants in the cis-region of target gene, c_i = s_i - 2p_iq_i, where s \in \{0,1,0\} and p_i represent the MAF of ith causal variants, and q_i = 1 - p_i, the vector \gamma_i represents the dominance effect sizes of casual variants on gene expression and \gamma_i from a normal distribution N(0, \frac{\sigma^2_{cis}}{2p_iq_i(1-2p_iq_i)}), where \sigma^2_{cis} is the cis genetic variance. We then scaled the sampled effect \gamma_i of ith casual variant by multiplying the factor \sigma^2_{cis} = Var(C\gamma). We generated the cis-heritability from a scaled Gamma distribution of 0.1*rand(\Gamma(0.4, 1.6)) and ranged from 0.005 to 0.3. We generated the polygenic_{trans} under three scenarios: 1) additive polygenic effects: K equals to the additive GRM (A); 2) dominant polygenic effects: K equals to the dominant GRM (D); and 3) mixed polygenic effects: K = r_1 A + r_2 D, of which the ratio r = r_1/r_2 was sampled from the empirical ratios between the additive and dominant genetic variances that estimated from the MAGE dataset.

Interaction cis-eQTL (cis-ieQTL) simulation
To evaluate the performance of OmiGA in context-dependent cis-eQTL effect (i.e., cis-ieQTL) mapping, we simulated genes with ancestry-specific effects for causal variants within cis-region in SIM-CLOSE and SIM-UNREL. We first estimated the proportion of the Yorkshire ancestry of 300 individuals from GIAD042 populations and the AFR ancestry of 731 individuals from MAGE populations using ADMIXTURE v1.3.0 based on genotype data.
To simulate genes with ancestry-specific effects for causal variants, we used the following formula to generate the expression levels (y) for each gene:
y = \mu + W\beta + polygenic_{trans} + \epsilon (28)
where \mu, polygenic_{trans}, and \epsilon is the same as those in formula (26), W is the matrix of interaction term between genotype and the proportion of ancestry with W = C \circ a, of which C is the matrix the standardized genotypes of m_{cis} causal variants in the cis-region of target gene, c_i = s_i - 2p_i, where s \in \{0,1,2\} and p_i represent the MAF of ith causal variants, a is the vector of the proportion of ancestry across individuals, \beta is the sampled effect sizes of causal variants follows the normal distribution of N(0, 1). We further scale the sampled effect sizes \beta by multiplying the factor \sigma^2_{int} = Var(W\beta), where \sigma^2_{int} is the genetic variance explained by the interaction between genotype and ancestry. We sample the number of causal variants on the basis of a prior probability table (Supplementary Data 2) and restrain the proportion of phenotypic variance explained by interaction term (h^2_{int}) to a range of 0.005–0.3. All other settings are consistent with formula (26), wherein h^2_{int} substitute for h^2_{cis}.
Based on the strategies above, we simulated 2459 ieGenes and 2976 non-ieGenes in the SIM-CLOSE, and 1513 ieGenes and 1802 non-ieGenes in the SIM-UNREL, where the ieGenes and the non-ieGenes are the genes simulated with at least one causal cis-variant and zero causal cis-variant, respectively. All the simulations were repeated 50 times.

Performance comparison for molQTL mapping
To quantify the performance of different methods, we computed several metrics including the genomic inflation factor (\lambda), false discovery rate (FDR), F1-score, and true positive rate (TPR). \lambda is defined as the median chi-squared statistic of genome-wide associations divided by its expected value for null variants. The FDR is defined as the proportion of simulated non-eGenes within these genes that were reported as significant (for example, with a P-value < 0.05). The F1-score is computed by the formula F = 2 * precision * recall / (precision + recall), where precision = TP / (TP+FP) and recall = TP / (TP+FN), in which TP, FP, and FN represent the number of true positive, false positive, and false positive genes, respectively. The TPR is defined as the proportion of simulated eGenes that can be discovered under 5% FDR.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors have developed a tool named OmiGA for ultra-efficeint molQTL mapping analysis, whilte accounting for population structure and relatedness. Both the paper and the software are well-constructed and clearly presented. I have some suggestions for the authors to consider:\n1. The work mostly compared OmiGA with tensorQTL, however, other tools such as fastGWA (REML-based method) that has been developed which is computationally efficient and can handle population structure and sample relatedness. Another very commonly used tool is plink2 and more recently Regenie. It would be very helpful for users to assess and understand your tool by including the comparisons with these other tools that are already well-implemented by other researchers. For example, to include the comparisons in simulation studies or data applications.\n2. The simulation study compared the inclusion of 0-10 PCs in tensorQTL, while the purpose is to demonstrate how robust the performance of OmiGA is even not including PCs as covariates in the model, it is a standard practice to adjust for at least 10 PCs in QTL analysis. So it is not clear how this simulation demonstrates the better performance compared to tensorQTL (if suppose we always do the standard practice). It would be helpful to further investigate/explain this.\n3. How is the scalability of OmiGA? For example, in human studies, there are large sample sizes in biobanks that have proteomics data, such as the UK Biobank, which may have tens of thousands of individuals (30K, many of them with some levels of relatedness as well).\n4. The method is capable of detecting non-additive effects. In the data applications, were any interesting non-additive effects observed?\n(Remarks on code availability)\nThe source codes are available on Github. I was able to download the summary statistics as well. It would be great to have a better documentation of the summary statistics (what each column represent). The website is nicely built.", "reviewer_2": "Summary:\nThis paper introduces OmiGA, a toolkit claimed for ultra-efficient molecular trait analysis in populations with complex relatedness. It is claimed that “existing tools often rely on simple linear models (LM)”, which fail to account for interindividual relatedness and population structure, leading to inflated false positives. OmiGA addresses this by leveraging linear mixed models (LMM) and algorithms (e.g., AIDUL, PAED) to improve computational efficiency and statistical power. Through simulations and real-data analyses (pig and human datasets), OmiGA demonstrates superior performance in controlling false discovery rates (FDR), enhancing eGene discovery, and refining causal variant localization.\nOverall, the paper is well-written and the github is well maintained. However, the novelty is over-claimed. The authors have neither used existing state-of-the-art algorithms of genomic LMM in the design, nor did they use modern GPU acceleration for implementation. As such, this reviewer is not convinced that the presenting work will be broadly used for a long time.\nMajor comments:\n(1) Despite the main claim of the work being the speed, there is a lack of incorporation of state-of-the-art algorithms that are frequently adapted in developing LMM tools for huge genomic datasets. In the manuscript, it is stated that the factorization of the GRM needs a complexity of O(n^3) (Line 112)-- this is a very old way to solve GRM, which leads to the main bottleneck of the computation. For instances, both FaST-LMM’s O(m^2*n) (Christoph Lippert et al, Nat Methods, 2011) and Bolt-LMM’s O(n*m) (Po-Ru Loh et al, Nat Genet, 2015) are way faster than O(n^3). (m is the number of SNP and n is the sample size.) They both did not factorize GRM using Eigen decomposition, instead they factorized the genotype matrix using SVD.\n(2) Not using GPU-based acceleration, which represents the current trend of new genomic tools development. Indeed, there is a lack of comparisons to existing GPU-accelerated tools, for instance tensorQTL (Amaro Taylor-Weiner et al, Genome Biol, 2019), which was actually mentioned by the authors as their Ref [9].\nMinor comments:\n(1) Overclaim on: “all the existing molQTL tools (Line 81)”.\n(2) Figure 1b: It claims the novelty of using a genetic background term in the regression formula – this is not true. Please refer to Kang HM et al Genetics 2008 EMMA as the first one that uses LMM for eQTL analysis.\n(3) Limited Validation of Non-Additive Effects: While models for dominant effects (e.g., LMM(d + A + D)) are proposed, their performance lacks empirical support. Additional simulations or case studies are needed.\n(4) Comparison of Multiple Testing Methods: A deeper comparison between ACAT/ClipperQTL and other methods (e.g., Bonferroni, adaptive FDR) would strengthen the argument for their superiority."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThe authors have developed a tool named OmiGA for ultra-efficeint molQTL mapping\nanalysis, whilte accounting for population structure and relatedness. Both the paper and\nthe software are well-constructed and clearly presented. I have some suggestions for\nthe authors to consider:\n1. The work mostly compared OmiGA with tensorQTL, however, other tools such as\nfastGWA (REML-based method) that has been developed which is computationally\nefficient and can handle population structure and sample relatedness. Another very\ncommonly used tool is plink2 and more recently Regenie. It would be very helpful for users\nto assess and understand your tool by including the comparisons with these other tools\nthat are already well-implemented by other researchers. For example, to include the\ncomparisons in simulation studies or data applications.\nAU: Thank you very much for the constructive suggestions. As suggested, we have\ncarefully compared OmiGA with the commonly used association analysis tools, including\nfastGWA1, PLINK22, Regenie3, GCTA-MLMA4, and GEMMA5, in a real gene expression\ndataset of pigs (i.e., the duodenum tissue of 300 pigs), which shows a much higher family\nrelatedness than that of human populations. As these tools are developed and optimized\nfor genome-wide association studies (GWAS) of complex traits, we benchmarked them\nwith OmiGA by conducting GWAS for gene expression. We found that OmiGA, GCTA-\nMLMA (no-preadj), and GEMMA outperformed other methods in controlling the global\npopulation structure (measured by lambda) in such a family-based or highly correlateddata,\nwhich is very common in livestock populations (Response Figure 1). Although OmiGA,\nGCTA-MLMA (no-preadj) and GEMMA had consistent results (Pearson’s correlations of\nassociation statistics between them close to 1) (Response Figure 2), OmiGA required\nsignificantly less runtime compared to GCTA-MLMA (no-preadj) and GEMMA under default\nsettings. More specifically, among these three tools, OmiGA was approximately 9 times\nfaster than GCTA-MLMA (no-preadj) and about 55 times faster than GEMMA. To speed up\ncomputations, fastGWA employs covariates-adjusted phenotypes and sparse genetic\nrelationship matrix (GRM) for genetic relatedness (after neglecting the lowest terms; e.g.,\nlower than 0.05), which works well in unrelated human populations but is less well-adapted\nto the complex and highly correlated livestock data, as demonstrated by our own results\n(reduced statistical power, Response Figure 1) and previous studies in dairy cattle6,7.\nIn addition, we would like to discuss that tools such as fastGWA, PLINK2, and Regenie\nare highly optimized for genome-wide association studies (GWAS) of complex traits, they\nare not directly suitable and user-friendly for molecular QTL (molQTL) mapping due to the\nbig differences between molQTL mapping and GWAS as below:\na) Differences in the analytical frameworks of GWAS and molQTL mapping8.\nGWAS is typically applied to non-molecular traits (e.g., body weight, height, and BMI),\nand a fixed P-value threshold (e.g., 5◊10-8) is generally adopted to determine\ngenome-wide significant associations. In contrast, molQTL mapping focuses on\ngenetic associations with molecular phenotypes (e.g., gene expression and\nalternative splicing) and typically identifies cis-QTL (SNP residing 1Mb around\ntranscript start site, TSS) of thousands of genes rather than conducting genome-\nwide screening due to limited sample sizes. Therefore, molQTL mapping requires\nmultiple testing correction at two levels (i.e., trait-level across the genome and\nvariant-level within each trait), which is different from that used in GWAS.\nb) Conventional linear mixed model (LMM)-based GWAS tools lack optimization\nfor high-throughput molecular phenotypes. These LMM-GWAS tools are not that\neasy to adapt in molQTL mapping for millions of molecular phenotypes. Moreover,\nkey aspects of molQTL mapping, such as eGene detection, independent test,\nGenotype × Environment interaction, and non-additive cis-QTL, as well as heritability\nestimation (including both cis- and trans-heritability), are often not integrated in a\nsingle toolkit.\nIn summary, we have incorporated these results and discussed them carefully in the\nrevised manuscript (Lines 338-340, 429-453, Revised Fig. 4a, Extended Data Fig. 6).\nResponse Fig. 1 | Genomic control inflation factor from genome-wide eQTL mapping\nusing different association tools. Analyses were performed based on 100 randomly\nselected genes in the duodenum of 300 pigs. MLMA (default): GCTA-MLMA with default\nparameters. MLMA (no-preadj): GCTA-MLMA with the ‘--mlma-no-preadj-covar’ option.\nLM: linear model. LMM (preadj): linear mixed model using covariates-adjusted phenotypes.\nLMM (raw): linear mixed model using raw phenotypes. ML: machine learning.\nResponse Fig. 2 | Correlations between the association P values of OmiGA and other\ntools. We computed Pearson’s correlations (r) based on the -log10(P) of\nENSSSCG00000038931 gene between different tools in the duodenum of 300 pigs. MLMA\n(default): GCTA-MLMA with default parameters. MLMA (no-preadj): GCTA-MLMA with the\n‘--mlma-no-preadj-covar’ option. The red line is a diagonal line with a slope of 1.\n2. The simulation study compared the inclusion of 0-10 PCs in tensorQTL, while the\npurpose is to demonstrate how robust the performance of OmiGA is even not including\nPCs as covariates in the model, it is a standard practice to adjust for at least 10 PCs in\nQTL analysis. So it is not clear how this simulation demonstrates the better performance\ncompared to tensorQTL (if suppose we always do the standard practice). It would be\nhelpful to further investigate/explain this.\nAU: Thanks for this insightful comment. As suggested, we have conducted additional\nanalyses using simulated datasets to compare OmiGA and tensorQTL with different\ngenotype PCs (0, 2, 5, 10) as shown in Response Fig. 3. (1) Accuracy of effect estimates\n(Response Fig. 3a): The absolute error of effect estimates from tensorQTL with different\nPCs was significantly higher than that from OmiGA (nPC=0), indicating that OmiGA\nprovided more precise effect size estimates. (2) Correlation with true effects (Response\nFig. 3b): Pearson’s correlation between estimated effects and true effects was significantly\nlower in tensorQTL (with different PCs) compared to OmiGA (nPC=0), further supporting\nOmiGA’s stronger ability to capture the true genetic effects. (3) Control of non-causal\nvariants (Response Fig. 3c): The mean Chi-square statistics of non-causal variants in\ntensorQTL (with different PCs) were significantly higher than those in OmiGA (nPC=0),\nsuggesting that tensorQTL still exhibited residual inflation in non-associated regions even\nafter adjusting for different PCs, whereas OmiGA better avoided such spurious signals. (4)\nThe optimal number of PCs: Although including more genotype PCs in tensorQTL\nenhances the power of eQTL mapping, it requires extensive testing to determine the\noptimalnumbers of PCs in the real-world data, as done in human GTEx study9. For\ndatasets with large sample sizes and complex relatedness, the process of selecting the\noptimal number of PCs might be computationally intensive in tensorQTL.\n(5) Causal\nprobability of fine-mapping for independent cis-eQTL from OmiGA was significantly higher\nthan tensorQTL in real data (Response Fig. 4).\nAltogether, these results indicate that OmiGA is user-friendly, robust in selecting\ngenotype PCs, and outperforms tensorQTL in the accuracy of effect estimation, particularly\nin the livestock populations with high family relatedness. We have incorporated these\nresults and carefully discussed them in the revised manuscript (Lines 297-300, Revised\nFig. 3d-f, Revised Extended Data Fig. 8).\nResponse Fig. 3 | Comparison of effect estimates in tensorQTL and OmiGA based\non different numbers of genotype principal components (PCs). a, Absolute error of\neffect estimates from tensorQTL and OmiGA. b, Pearson’s correlation between estimated\nand true effects. c, The mean Chi-square statistics of causal and non-causal variants.\nResponse Fig. 4 | Comparison of causalposterior inclusion probability (PIP) of\nindependent cis-eQTL obtained from OmiGA and tensorQTL in different tissues of\npig dataset (GIAD042). PIP was computed using SuSiE-infand only these genes with 95%\ncredible set were retained. The point is the average PIP. The P-value was calculated from\na two-sided Student’s t-test.\n3. How is the scalability of OmiGA? For example, in human studies, there are large sample\nsizes in biobanks that have proteomics data, such as the UK Biobank, which may have\ntens of thousands of individuals (30K, many of them with some levels of relatedness as\nwell).\nAU: OmiGA supports analysis of the biobank-scale dataset efficiently via the ‘--chunk-size’\nparameter. As suggested by Reviewer #2 as well, we have extended OmiGA to support the\nGPU acceleration, thus further speeding up the processing of biobank-scale datasets. In\naddition, we evaluated OmiGA’s runtime on simulated large-scale datasets of 10,000,\n20,000 and 30,000 individuals, including ~20,000 expressed genes and ~10 million genetic\nvariants.\nOur results demonstrated that OmiGA exhibits strong scalability (Response Fig. 5): In\noutperformed the LM-based tensorQTL in terms of runtime when processing this large\ndataset. With GPU, OmisGA can further speed up. These findings indicate that OmiGA,\nparticularly with GPU, is well-suited for the scale of human studies with tens of thousands\nof individuals (e.g., proteomics data in UK Biobank), even when accounting for relatedness\namong individuals.\nWe incorporated these results in the revised manuscript (Lines 198-205, Revised\nExtended Data Fig. 1).\nResponse Fig. 5 | Runtime of cis-eQTL mapping using OmiGA and tensorQTL in\ndatasets with different sample sizes. LM: linear model. LMM: linear mixed model. All\ntests were performed in the same computing environment: 120 GB of memory, 10 CPU\ncores (Intel Xeon CPU Max 9462), 1 GPU (GeForce RTX 4090 GPU) in one computer\nnode.\n4. The method is capable of detecting non-additive effects. In the data applications, were\nany interesting non-additive effects observed?\nAU: In this revision, we have added an additional case study to present new findings\nregarding non-additive effects in Lines 390-404 (Response Fig. 4-5). In both human and\npig datasets, we observed that most (～71%) eGenes identified by the additive model also\nharbor dominant cis-eQTL (Response Fig. 6a & 7a-c). Compared to the additive model,\nthe dominant model detected additional 340 eGenes with dominant cis-eQTL only in the\nMAGE dataset (Response Fig. 6a). The dominant lead cis-eQTL were located farther from\nthe TSS of their target genes than additive lead cis-eQTL (Response Fig. 6b, 7d-f). For\ninstance, we identified a significant dominant cis-eQTL associated with PKD1 in LCLs, but\nno significant additive cis-eQTL was identified (Response Fig. 6c-d). Previous studies\nhave reported that PKD1 is linked to autosomal dominant Polycystic Kidney Disease 1 10,11.\nWhile most eGenes had both additive and dominant effects, we observed dominant-\nspecific cis-eQTL within these genes. For instance, SNP 1:11927646:A:T, the dominant\nlead cis-eQTL of TNFRSF8, is independent of the additive lead cis-eQTL (Response Fig.\n6e). These results demonstrated that the dominant model implemented in OmiGA can\nprovide novel insights into the non-additive genetic regulatory of gene expression.\nWe incorporated all these results in the revised manuscript (Lines 390-404, Revised Fig.\n6, Revised Extended Data Fig. 10).\nResponse Fig. 6 | Dominant cis-eQTL discovery in MAGE dataset. a, -log10-scaled P-\nvalues of additive and dominant lead cis-variants detected in LCL from MAGE dataset. b,\nDistance to the transcription start site (TSS) of target genes for additive and dominant lead\ncis-eQTL. The P value was obtained by the two-sided Wilcoxon rank-sum test. c,\nNormalized expression of PKD1 for three genotypes of 16:1829064:A:T. d, -log10 P-values\nof cis-variants from dominant (Dom, top) and additive (Add, bottom) models for PKD1. e, -\nlog10 P-values of cis-variants from dominant (Dom, top) and additive (Add, bottom) models\nfor TNFRSF8.\nResponse Fig. 7 | Dominant cis-eQTL discovery in GIAD042 pig dataset. a-c, -log10-\nscaled P-values of additive and dominant lead cis-variants detected in duodenum (a),\nliver (b), and muscle (c). d-f, Distance to the transcription start site (TSS) of target\ngenes for additive and dominant lead cis-eQTL in duodenum (d), liver (e), and muscle (f). The P\nvalue was obtained by the two-sided Wilcoxon rank-sum test.\nRemarks on code availability:\nThe source codes are available on Github. I was able to download the summary statistics\nas well. It would be great to have a better documentation of the summary statistics (what\neach column represent). The website is nicely built.\nAU: Thank you very much for the valuable suggestion. We have added a detailed\nexplanation of the summary statistics (https://omiga.bio/#/molQTL?id=output-field-\ndefinition).\nReviewer #2 (Remarks to the Author):\nSummary:\nThis paper introduces OmiGA, a toolkit claimed for ultra-efficient molecular trait analysis in\npopulations with complex relatedness. It is claimed that “existing tools often rely on simple\nlinear models (LM)”, which fail to account for interindividual relatedness and population\nstructure, leading to inflated false positives. OmiGA addresses this by leveraging linear\nmixed models (LMM) and algorithms (e.g., AIDUL, PAED) to improve computational\nefficiency and statistical power. Through simulations and real-data analyses (pig and\nhuman datasets), OmiGA demonstrates superior performance in controlling false discovery\nrates (FDR), enhancing eGene discovery, and refining causal variant localization.\nOverall, the paper is well-written and the github is well maintained. However, the novelty is\nover-claimed. The authors have neither used existing state-of-the-art algorithms of\ngenomic LMM in the design, nor did they use modern GPU acceleration for implementation.\nAs such, this reviewer is not convinced that the presenting work will be broadly used for a\nlong time.\nAU: Thank you very much for the valuable and constructive comments. As suggested, we\nhave carefully revised and improved the novelty of OmiGA in this revised version. More\nspecifically, we have implemented the state-of-the-art algorithms of genomic LMM, such\nas IDUL (iterative dispersion updates to fit LMMs)12, and also adapted the modern GPU\nacceleration to OmiGA (https://omiga.bio/#/Installation?id=gpu-acceleration). As also\nsuggested by Reviewer #1, we systematically compared OmiGA with widely used\nassociation study tools, including fastGWA, PLINK2, Regenie, GCTA-MLMA, GEMMA, in\nthe GIAD042 dataset of 300 pigs. As shown in Response Fig. 1, we demonstrated that the\nOmiGA can outperform most of them regarding computational efficiency and estimate\naccuracy. In addition, through comparing OmiGA and the GPU-based tensorQTL, we\nsimulated large datasets of 10,000, 20,000 and 30,000 individuals, and demonstrated the\nefficiency of the algorithm and scalability of OmiGA for processing large-scale datasets\n(Response Fig. 5).\nAs reviewed previously8, there are some differences between molQTL mapping and GWAS\nof complex traits such as number of tested phenotypes (millions of molecular phenotypes\nvs. a few complex traits), and definition of significance (two-steps’ multiple-testing\ncorrection methods in molQTL vs. fixed cutoff e.g., 5×10-8). Therefore, we have carefully\ndiscussed the limitations of most existing LMM-based tools, such as EMMA (Kang HM et\nal., 2008), EMMAX (Kang HM et al., 2010), fastGWA, Regenie, and GEMMA, in molQTL\nmapping, as they were initially developed for GWAS of complex traits:\na) Conventional LMM-based GWAS tools lack optimization for high-throughput\nmolecular phenotypes. As a result, using these tools for molQTL mapping results in\na substantial number of redundant computations across tens of millions of molecular\nphenotypes, including the factorization of the GRM. Moreover, many key aspects of\nmolQTL mapping (as listed in refs8), such as conditionally independent test, genotype-\nby-environment interaction, and non-additive cis-QTL, as well as heritability (including\nboth cis- and trans-heritability), have not been integrated in a single analysis toolkit as\nwe did in OmiGA. As suggested by Reviewer #1, we also systematically compared the\nperformance of OmiGA with these commonly used GWAS tools regarding the\ncomputational efficiency and estimate accuracy in GWAS of gene expression, and\ndemonstrated that the OmiGA outperformed most of them in both simulated and real-\nworld data (Lines 338-340, Response Fig. 1).\nb) Substantial differences exist between the analytical frameworks of GWAS and\nmolQTL mapping8. GWAS is typically applied to organismal traits (e.g., body weight,\nheight, and BMI), and a fixed P-value threshold (e.g., 5◊10-8) is generally adopted to\ndetermine genome-wide significant associations. In contrast, molQTL mapping focus\non genetic associations with millions of molecular traits (e.g., gene expression and\nalternative splicing) and typically identifies cis-QTL rather than conducting genome-\nwide screening due to limited sample size. Therefore, molQTL mapping requires\nmultiple testing correction at two levels (i.e., trait-level across the genome and\nvariant-level within each trait), which is distinct from that used in GWAS.\nAs outlined in the Farm animal Genotype-Tissue Expression (FarmGTEx) perspective\npaper, in the coming years, many large-scale molQTL studies will be conducted across a\nrange of livestock species13, where high genetic relatedness among such animals is very\ncommon due to the small effective population size. In addition, when coming to biobank-\nscale datasets in humans, molQTL mapping also has to address the cryptic relatedness\namong individuals to increase statistical power14–16. Therefore, we believe that OmiGA will\nbe broadly used in the near future.\nIn summary, we have incorporated this information in the revised manuscript (Lines 338-\n340, 429-452, Revised Fig. 4a, Revised Extended Data Fig. 1).\nResponse Fig. 1 | Genomic control inflation factor from genome-wide eQTL mapping\nusing different association study tools. Analyses were performed based on 100\nrandomly selected genes in the duodenum of GIAD042. MLMA (default): GCTA-MLMA with\ndefault parameters. MLMA (no-preadj): GCTA-MLMA with the ‘--mlma-no-preadj-covar’\noption. LM: linear model. LMM (preadj): linear mixed model using covariates-adjusted\nphenotypes. LMM (raw): linear mixed model using raw phenotypes. ML: machine learning.\nOmiGA is approximately 9 times faster than MLMA (no-preadj) and about 55 times faster\nthan GEMMA.\nResponse Fig. 5 | Runtime of cis-eQTL mapping using OmiGA and tensorQTL in\ndatasets with different sample sizes. LM: linear model. LMM: linear mixed model. All\ntests were performed in the same computing environment: 120 GB of memory, 10 CPU\ncores (Intel Xeon CPU Max 9462), 1 GPU (GeForce RTX 4090 GPU) in one computer\nnode.\nMajor comments:\n(1) Despite the main claim of the work being the speed, there is a lack of incorporation of\nstate-of-the-art algorithms that are frequently adapted in developing LMM tools for huge\ngenomic datasets. In the manuscript, it is stated that the factorization of the GRM needs a\ncomplexity of O(n^3) (Line 112)-- this is a very old way to solve GRM, which leads to the\nmain bottleneck of the computation. For instances, both FaST-LMM’s O(m^2*n) (Christoph\nLippert et al, Nat Methods, 2011) and Bolt-LMM’s O(n*m) (Po-Ru Loh et al, Nat Genet, 2015) are way faster than O(n^3). (m is the number of SNP and n is the sample size.) They\nboth did not factorize GRM using Eigen decomposition, instead they factorized the\ngenotype matrix using SVD.\nAU: Thank you for the valuable comments. As suggested, we have replaced the Eigen\nDecomposition with SVD algorithm in OmiGA (Line 112). Additionally, we demonstrated\nthat the impact of GRM factorization on the overall computational time is negligible in\nOmiGA as below.\nThe factorization of the GRM is a one-time operation in OmiGA. While its time complexity\nis O(n^3), the one-time GRM factorization does not dominate the overall running time of\nthe entire analysis. Our empirical tests further confirmed this: when factorizing the GRM of\n10,000 samples (a sample size exceeding that of most current studies) using SVD, the\nprocess took only ~30 seconds with 10 CPU threads (Intel(R) Xeon(R) CPU Max 9462)\n(see Response Table 1). Even for a large-scale dataset (n=30,000), GRM factorization time\nremains manageable in comparison to the time required for subsequent iterative analyses.\nHowever, we will keep OmiGA updated timely. When sample size of molQTL mapping\nreaches >100,000 individuals, we will utilize more efficient algorithms in the future, such\nas OCMA17.\nWe acknowledge the reviewer's concern again, and have carefully discussed this in the\nrevised manuscript (Lines 429-452).\nResponse Table 1 | Runtime (second) of GRM factorization in different sample sizes\nusing SVD in OmiGA. Analyses were performed in the same computing environment: 120\nGB of memory, 10 CPU cores (Intel Xeon CPU Max 9462) in one computer node.\nN 100 200 500 1,000 2,000 5,000 10,000 20,000 30,000\nSec. 0.001 0.004 0.022 0.081 0.338 5.083 27.833 189.356 532.289\n(2) Not using GPU-based acceleration, which represents the current trend of new genomic\ntools development. Indeed, there is a lack of comparisons to existing GPU-accelerated\ntools, for instance tensorQTL (Amaro Taylor-Weiner et al, Genome Biol, 2019), which\nwas actually mentioned by the authors as their Ref [9].\nAU: Thank you for the valuable suggestion. As suggested, we have adapted the GPU-\nbased acceleration to OmiGA and compared it with GPU-accelerated tensorQTL carefully\nusing simulated large-scale datasets of 10,000, 20,000 and 30,000 individuals. When using\noutperformed the LM-based tensorQTL in terms of runtime; When using linear regression\n(LM) model, both CPU and GPU versions of OmiGA were faster than the GPU-accelerated\ntensorQTL (Response Fig. 5).\nWe have incorporated this information in the revised manuscript (Lines 198-205 and\nRevised Extended Data Fig. 1).\nResponse Fig. 5 | Runtime of cis-eQTL mapping using OmiGA and tensorQTL in\ndatasets with different sample sizes. LM: linear model. LMM: linear mixed model. All\ntests were performed in the same computing environment: 120 GB of memory, 10 CPU\ncores (Intel Xeon CPU Max 9462), 1 GPU (GeForce RTX 4090 GPU) in one computer\nnode.\nMinor comments:\n(1) Overclaim on: “all the existing molQTL tools (Line 81)”.\nAU: Thank you for pointing this out. We have revised the sentence as: “we demonstrated\nthat OmiGA outperformed LM-based tensorQTL regarding the statistical discovery power\nof molQTL and computational efficiency, …”. Please see Lines 82-83.\n(2) Figure 1b: It claims the novelty of using a genetic background term in the regression\nformula – this is not true. Please refer to Kang HM et al Genetics 2008 EMMA as the first\none that uses LMM for eQTL analysis.\nAU: Thank you for pointing this out. We have removed this in the revised Figure 1, and\nrevised the manuscript accordingly (see Revised Fig. 1).\n(3) Limited Validation of Non-Additive Effects: While models for dominant effects (e.g.,\nLMM(d + A + D)) are proposed, their performance lacks empirical support. Additional\nsimulations or case studies are needed.\nAU: Thank you for the constructive suggestions. In this revision, we have added additional\nsimulations and case studies to demonstrate the performance of OmiGA in non-additive\nmolQTL discovery.\nSimulations (Lines 315-319, 1249-1268): In terms of dominant cis-eQTL discovery,\nconsistent with the observations for additive effects, we found that LMMs exhibited a higher\ntrue positive rate (TPR) than LM (Response Fig. 8). The dominant model incorporating\nboth additive and dominant polygenic effects (i.e., LMM(d+A+D)) showed a greater\nrobustness than other models (simple linear model, LMM(d+A), and LMM(d+D)) across\nvarious simulation scenarios as shown in Response Fig. 8.\nCase studies (Lines 389-403): We have added case studies to present findings\nregarding non-additive effects (Response Fig. 9-10). In both human and pig datasets, we\nobserved that average of 71% eGenes that were identified by the additive model also\nharbored significant dominant cis-eQTL (Response Fig. 9a & 10a-c). Compared to the\nadditive model, the dominant model detected additional 340 eGenes with only dominant\ncis-eQTL in the MAGE dataset (Response Fig. 9a). The dominant lead cis-eQTLs were\nlocated farther from the TSS of their target genes than additive lead cis-eQTLs (Response\nFig. 9b, 10d-f). For instance, we identified a significant dominant cis-eQTL associated with\nPKD1 in LCLs, but no significant additive cis-eQTL was detected (Response Fig. 9c-d).\nMultiple studies have reported that PKD1 is linked to autosomal dominant Polycystic\nKidney Disease 1 10,11. While most eGenes were detectable by both additive and dominant\nmodels, we observed dominant-specific cis-eQTL within these eGenes. For instance, SNP\n1:11927646:A:T, the dominant lead cis-eQTL of TNFRSF8, is independent of the additive\nlead cis-eQTL (Response Fig. 9e). These results demonstrated that the dominant model\nimplemented in OmiGA can provide new insights into the non-additive genetic regulatory\nof gene expression.\nWe incorporated these results in the revised manuscript (Lines 315-319, 389-403, 1249-\n1268, Revised Fig. 3j&6, Revised Extended Data Fig. 10, Revised Supplementary Fig. 8).\nResponse Fig. 8 | True positive rate (TPR) under 5% FDR for dominant cis-eQTL\ndiscovery using different statistical models in OmiGA. a,b, Analyses were performed\nbased on SIM-CLOSE (a) and SIM-UNREL (b) datasets. Additive, dominant, and mixed\npolygenic effects represent the simulated gene expression with different background\npolygenic effects, which contain additive only, dominant only, and both additive and\ndominant effects, respectively.\nResponse Fig. 9 | Dominant cis-eQTL discovery in MAGE dataset. a, -log10-scaled P-\nvalues of additive and dominant lead cis-variants detected in LCL from MAGE dataset. b,\nDistance to the transcription start site (TSS) of target genes for additive and dominant lead\ncis-eQTL. The P value was obtained by the two-sided Wilcoxon rank-sum test. c,\nNormalized expression of PKD1 for three genotypes of 16:1829064:A:T. d, -log10 P-values\nof cis-variants from dominant (Dom, top) and additive (Add, bottom) models for PKD1. d, -\nlog10 P-values of cis-variants from dominant (Dom, top) and additive (Add, bottom) models\nfor TNFRSF8.\nResponse Fig. 10 | Dominant cis-eQTL discovery in the GIAD042 pig dataset. a-c, -\nlog10-scaled P-values of additive and dominant lead cis-variants detected in duodenum\n(a), liver (b), and muscle (c). d-f, Distance to the transcription start site (TSS) of target\ngenes for additive and dominant lead cis-eQTL in duodenum (d), liver (e), and muscle (f).\nThe P value was obtained by the two-sided Wilcoxon rank-sum test.\n(4) Comparison of Multiple Testing Methods: A deeper comparison between\nACAT/ClipperQTL and other methods (e.g., Bonferroni, adaptive FDR) would strengthen\nthe argument for their superiority.\nAU: Thank you for the constructive suggestion. As suggested, we have conducted a deep\ncomparison between ACAT/ClipperQTL and four additional methods—adaptive FDR\n(Storey), Bonferroni (Bonf), Benjamini & Hochberg (BH), and Benjamini & Yekutieli (BY) —\nfor multiple-testing correction, in the simulated datasets. We have provided these methods\nin OmiGA for users to choose flexibly. We summarized the F1-scores of different\napproaches in Response Fig. 11. In additive and dominant cis-eQTL mapping, we\nobserved all tested multiple testing methods resulting in an FDR<0.05. Under 5% FDR,\nClipperQTL and ACAT exhibited the highest F1-scores in both simulated datasets (SIM-\nCLOSE and SIM-UNREL). For interaction cis-eQTL mapping, the multiple-testing methods\nshould be carefully selected, while ACAT was recommended as the default method in\nOmiGA.\nWe incorporated and discussed these results in the revised manuscript (Lines 292, 497-\n502 and Revised Extended Data Fig. 4).\nResponse Fig. 11 | Comparison of multiple testing approaches for eQTL mapping\nusing OmiGA. a-f, F1-scores of different approaches for additive cis-eQTL, dominant cis-\neQTL, and interaction cis-eQTL mapping in SIM-CLOSE (a-c) and SIM-UNREL (d-f)\nscenarios. *:p: Gene-level P-values were calculated using different methods and the\nsignificant threshold was defined as gene-level P-values < 0.05. *:q: Benjamini & Hochberg\n-corrected gene-level P-values < 0.05 across all tested genes. Specially for Clipper,\ngene-level P-values (Clipper:p) were calculated using its permutation algorithm and\ncorrected by ClipperQTL (Clipper:q). ACAT: Aggregated Cauchy Association Test, the gene-level P-\nvalues were computed via the ACAT method using all cis-variants within each gene.\nStorey-*: Storey’s q-values with different lambda parameters. BH: Benjamini & Hochberg\nmethod. Bonf: Bonferroni correction. BY: Benjamini & Yekutieli method. For approaches\nprefixed with ‘Storey’, ‘BH’, ‘Bonf’, or ‘BY’: Gene-level P-values were defined as the\ncorrected P-value of the lead variant (within each gene) using the respective method.\nFDR>5% and ≤5% indicates whether the observed false discovery rate is significantly\ngreater than 5% or not based on one-side Student’s t-test."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have fully addressed all my comments. The toolkit is well developed and I have no further questions.\n(Remarks on code availability)\nThe documentations are clear and the website for the tool is well constructed and easy to navigate. The source code is also available on Github. I was able to follow the documentation to download and install the tool in Linux environment and it worked smoothly.", "reviewer_2": "The authors have largely addressed my previous comments, adding necessary method comparisons and additional simulations and case studies. Overall, these improvements strengthen the manuscript.\n(Remarks on code availability)\nThe code is well maintained."}, "author_response": ""}]

Decision letter
A version of this paper was originally rejected for publication by Nature Communications, however that decision was reconsidered after appeal by the authors.