Association of expression levels with probabilities of imputed genotypes were tested in samples of related individuals using a two-step mixed model-based score test developed in the works of Aulchenko et al. and Chen and Abecasis 33 (link),34 (link) and implemented in the GenABEL/ProbABEL packages 35 (link),36 (link). Briefly, the approach is an approximation of a full linear mixed model where the first step includes a mixed model containing all terms but those involving SNPs fitted by maximum-likelihood i.e. fixed effects as well as kinship matrix based on genomic data. Fixed effects included age and experimental batch in the adipose and LCL analysis, while age, batch and sample processing were used in the skin analysis. This step is performed using the GenABEL software 35 (link)using the polygenic() function. The resulting object contains the inverse variance-covariance matrix of the estimates and expression trait residuals which are used in the second step together with posterior genotypic probabilities performing the score test in ProbABEL 36 (link), using the “—mmscore” option. In total, 776 adipose, 777 LCL and 667 skin samples had both expression profiles and imputed genotypes and were included in the analysis. Cis analysis was limited to SNPs located within 1MB either side of the transcription start or end site or within the gene body. False discovery rate (FDR) for the cis analysis was calculated from the complete list of p values using the qvalue package20 (link) implemented in R2.11 37 . In order to characterize likely independent regulatory effects, the identified cis-eQTLs were mapped to recombination hotspot intervals 38 (link). For each gene, the most significant SNP per hotspot interval were selected followed by additional LD filtering (for each remaining SNP pair with D’ > 0.5 and r2>0.1 the least significant SNP was ignored).
Trans analysis was limited to SNPs located on a different chromosome than the tested transcript. Post-QC analysis of the trans-eQTLs revealed 52 probes with extreme outlier effects which were filtered from further trans analysis. Transcripts associated with a trans-SNP at P<5×10−8 were used for calculations of transcript-wise FDR from the complete list of p values using the qvalue package20 (link) implemented in R2.11 37 .
The score test is known to slightly underestimate additive effect sizes 34 (link), so the top association per probe was validated with a linear mixed-effects model in R, using the lmer() function in the lme4 package 39 , fitted by maximum-likelihood(Supplementary Fig. 2). The linear mixed-effects model were adjusted for both fixed (age, experimental batch effect and sample processing effect (skin tissue only)) and random effects (family relationship and zygosity). A likelihood ratio test was applied to assess the significance of the SNP effect. The p-value of the SNP effect in each model was calculated from the Chi-square distribution with 1 degree of freedom using -2log(likelihood ratio) as the test statistic.