of $k$ transcripts, proteins or metabolites in $n$ individuals are organized in an $s\times n$ genotype matrix ${\mathbf{G}}$ and $k\times n$ expression data matrix ${\mathbf{X}}$. Genetic markers take values $0,1,\dots,\ell$, where $\ell$ is the maximum number of alleles ($\ell=2$ for biallelic markers), while molecular traits take continuous values. In the linear model, a linear relation is tested between the expression level of gene $i$ and the genotype value (i.e. the number of reference alleles) of SNP $j$. The corresponding test statistic is the Pearson correlation between the $i$th row of ${\mathbf{X}}$ and the $j$th row of ${\mathbf{G}}$, for all values of $i$ and $j$. Standardising the data matrices to zero mean and unit variance, such that for all $i$ and $j$, $$\begin{aligned}
  \sum_{l=1}^nX_{il} = \sum_{l=1}^nG_{jl} = 0 \quad\text{and}\quad
   \sum_{l=1}^n X_{il}^2 = \sum_{l=1}^n G_{jl}^2 = n ,     \end{aligned}$$ it follows that the correlation values can be computed as $$\begin{aligned}
  R_{ij} = \sum_{l=1}^nX_{il}G_{jl} = ({\mathbf{X}}{\mathbf{G}}^T)_{ij},\end{aligned}$$ where ${\mathbf{G}}^T$ denotes the transpose of ${\mathbf{G}}$. Hence, a single matrix multiplication suffices to compute the test statistics for the linear model for all pairs of traits and SNPs.

The ANOVA models test if expression levels in different genotype groups originate from the same distribution. Therefore, ANOVA 