For Illumina arrays, the raw unprocessed expression matrix was exported from GenomeStudio. First, two PCs were calculated on the quantile-normalized and log2 transformed expression data and plotted to identify and exclude outlier samples. The data were normalized in several steps: quantile normalization, log2 transformation, probe centering and scaling by the equation ExpressionProbe,Sample = (ExpressionProbe,Sample - MeanProbe) / Std.Dev.Probe. Genes showing no variance were removed. Next, the first four multidimensional scaling (MDS) components, calculated based on non-imputed and pruned genotypes using plink v1.0757, were regressed out of the expression matrix to account for population stratification. We further removed up to 20 of the first expression-based PCs that were not associated with any SNPs, as these capture non-genetic variation in expression. After regressing out these covariates, the residual gene expression matrix was used for eQTL mapping. Each cohort also ran MixupMapper58 software to identify incorrectly labeled genotype–expression combinations and resolved any sample mix-ups (https://github.com/molgenis/systemsgenetics/wiki/Resolving-mixups).
