or aggregating measurements across multiple tissues (Flutre et al., 2013). In practice, this requires millions of association tests in order to scan all possible phenotype-variant pairs in cis (i.e. variants located within a specific window around a phenotype), resulting in millions of nominal P-values. Matrix eQTL (Shabalin, 2012) has recently emerged as a ‘gold standard’ for this task (GTEx Consortium, 2015; Lappalainen et al., 2013) by taking advantage of efficient matrix operation implementations to perform the many association tests in acceptable running times. Due to the large number of tests performed per phenotype, multiple testing has to be accounted for to assess the significance of any discovered candidate QTL. A first naive solution to this problem is to correct the nominal P-values for the number of tested variants using the Bonferroni method. However, due to the specific and highly variable nature of each genomic region being tested in terms of allele frequency and linkage disequilibrium (LD), the Bonferroni method usually proves to be overly stringent and results in many false negatives. To overcome this issue, a commonly adopted approach (Montgomery et al., 2010) is to analyze thousands of permuted datasets for each phenotype in order to empirically characterize the null
