Scaled gene counts were filtered for missing values (maximum 90% missing values, i.e. zero counts, allowed per gene). Zero values were offset by 1, after which the data was log10 transformed and quantile normalized across individuals using R limma normalizeQuantiles function 71. We then ran PEER 72 with full pre-normalized dataset with the following parameters: K=30; covariates = gender, iPSC growth condition (feeder-dependent/E8), mean expression (‘addMean=True’ in PEER); maximum iterations = 10,000. Residuals for each gene were gaussianised, i.e. converted to the quantiles of a standard normal distribution and finally mean centered and standardized prior to mapping. In total, we had 26,936 and 17,116 genes available for mapping (RNA-seq and ‘gexarray’, respectively).