{'fraction_non_alphanumeric': 0.03612817244176614, 'fraction_numerical': 0.021395874377100475, 'mean_word_length': 5.010186313773289, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 18, 'https://': 2, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Common genetic variation drives molecular heterogeneity in human iPSCs

Induced pluripotent stem cell (iPSC) technology has enormous potential to provide improved cellular models of human disease. However, variable genetic and phenotypic characterisation of many existing iPSC lines limits their potential use for research and therapy. Here, we describe the systematic generation, genotyping and phenotyping of 711 iPSC lines derived from 301 healthy individuals by the Human Induced Pluripotent Stem Cells Initiative (HipSci: http://www.hipsci.org). Our study outlines the major sources of genetic and phenotypic variation in iPSCs and establishes their suitability as models of complex human traits and cancer. Through genome-wide profiling we find that 5-46% of the variation in different iPSC phenotypes, including differentiation capacity and cellular morphology, arises from differences between individuals. Additionally, we assess the phenotypic consequences of rare, genomic copy number mutations that are repeatedly observed in iPSC reprogramming and present a comprehensive map of common regulatory variants affecting the transcriptome of human pluripotent cells.


Introduction
Induced pluripotent stem cells (iPSCs) are important model systems for human disease1. A major open question is whether iPSCs can be used to study the functions of genetic variants associated with complex traits and normal human phenotypic variation. Previous work has suggested that individual iPSC lines are highly heterogeneous2-5, although some of these differences may arise due to genetic background of the donor6,7. Nonetheless, high variability could make iPSCs unsuitable cellular models for genetic variants with small effects. Existing iPSC lines also frequently have limited genetic and phenotypic data of variable quality, or are derived from individuals with severe genetic disorders, limiting their utility for studying other phenotypes.
The Human Induced Pluripotent Stem Cells Initiative (HipSci: www.hipsci.org) was established to generate a large, high-quality, open-access reference panel of human iPSC lines. A major focus of the initiative is the systematic derivation of iPSCs from hundreds of healthy volunteers using a standardised and well-defined experimental pipeline. The lines are extensively characterised and available to the wider research community along with the accompanying genetic and phenotypic data. Here, we report initial results from the characterization of the first 711 iPSC lines derived from 301 healthy individuals. We provide approach10 that can detect genetic abnormalities of >200 Kb occurring in 20% or more cells. We identified trisomies in 4% of lines (none of the selected lines), and 41% of lines (18% of the selected lines) harboured one or more CNAs of, on average, 7.15 Mb in length with duplications outnumbering deletions by 2.8 to 1 (Fig. 1e, Supplementary Table 2). Although the majority of CNAs were unique to single iPSC lines, 22% were also observed in at least one replicate line from the same donor (at least one base pair overlap), and 15% were identified in all replicates (Fig. 1f). We found no significant association between the number of CNAs and either passage number, donor age, gender or PluriTest score of a line (P > 0.09, Fig. 1g, Extended Data Fig. 3).
CNAs observed in pluripotent stem cells (PSCs) are known to recur at certain genomic locations11-13. We observed 35 regions where CNAs occurred significantly more often than expected under a uniform genomic distribution, including whole chromosome duplication of the X chromosome (P = 1.5x10 -9 ), 20 sub-chromosomal duplications, 11 deletions and three regions with both duplications and deletions (Fig. 2a, Supplementary Table 2). The three most frequent CNAs (X trisomy, chromosome 17 and 20) have been previously observed in PSCs12,14,15, but others are newly identified, to our knowledge.
Although recurrent CNAs could be due to mutational hotspots we did not find a significant overlap between our recurrent CNA set and annotated chromosome fragile sites16 (17% overlap, P = 0.075). Recurrent CNAs could also arise if duplication or deletion of specific genes led to a selective advantage. To identify potential targets of selection, we defined peak regions of amplification (regions of maximum recurrence e.g. Fig. 2c and Extended Data Fig. 4) within each CNA and identified expressed genes (FPKM>0 in >10% of lines). Fourteen candidate regions contained fewer than six expressed genes including genes with establised roles in cancer progression (DOCK1, FATS, WWOX, STAG2 and XIAP) 17-21. In regions with larger numbers of genes we searched for: (i) significant differential expression between lines with copy number 2 and 3 (ii) reported oncogenes from COSMIC20 and (iii) high scoring genes (top 2%) in a genome-wide siRNA screen for hESC cell proliferation22 (Fig. 2c, Extended Data Fig. 4; Supplementary Table 2). This approach identified BCL2L112 on chr20q11.21, EIF4A3, NOL11 and seven other genes on chr17q and UTP6 and SUZ12 on chr17q11.2. One candidate, EIF4A3, scored more highly than BCL2L1 in reducing ESC proliferation (Fig. 2c, top 0.1% of genes), was highly expressed in iPSCs, and over-expressed in lines with increased copy number, at both the mRNA (Q = 2x10 -5 ) and protein level (Extended Data Fig. 5). Finally, we compared lines from the same donor with and without CNAs to test for genome-wide effects on gene expression levels and, in a subset of cases, for effects on cell growth, proliferation and apoptosis (Fig. 2b, Extended Data Fig. 5). The recurrent duplication on chromosome 17 was associated with the largest number of changes in gene expression, including 1,098 genes (FDR < 1%) in trans located on other chromosomes, which were enriched for 'Neural Crest Differentiation' and 'DNA strand elongation' pathways (PathCards23, Supplementary Table 2). We also detected significant increases and decreases in cell growth rate associated with CNAs on chromosome 17 and 20 (Extended Data Fig. 5).
We further partitioned iPSC gene expression variation using the Tier 1 expression array data, the assay with the largest number of donors and lines. Of the 25,434 probes analysed (16,829 genes) (Supplementary Table 3), donor effects explained the largest proportion of variation in 46.4% of probes (53.3% of genes), substantially more than any other factor, including copy number status (23.4%), culture conditions (26.2%), passage (2%) and gender (1.9%, Fig. 3d). Donor effects were common, and consistent across large numbers of genes, while others such as CNA status had larger effects on a smaller number of genes (Fig. 3d). We observed minor effects of gender and line passage number on RNA-seq, methylation and protein immunofluorescence (Fig. 3d, Extended Data Fig. 6). Likewise, we did not observe substantial changes in PluriTest scores, or pluripotency marker expression across passages (P > 0.3, Extended Data Fig. 6), reflecting that pluripotency was maintained during culture. In principle, the estimated donor variation could arise due to shared reprogramming environment because lines were derived from the same population of fibroblast cells. However, expression quantitative trait locus (eQTL) effect sizes mapped using Tier 1 expression data (Supplementary Table 4, Extended Data Fig. 6) revealed that higher donor variation was associated with larger effect sizes of lead eQTL variants (Fig. 3e), suggesting that donor variance primarily reflects genetic differences.

Identification and characterization of iPSC-specific regulatory variants
Using RNA-seq data from 166 unrelated donors (median sequencing depth 38M reads), we next mapped eQTLs in a 1 Mb cis-window from the gene start. We identified 6,631 genes with an eQTL (FDR 5%, hereafter 'eGenes'), 598 of which had a significant secondary eQTL (Supplementary Table 4). Power to discover eGenes in iPSCs was comparable to that in somatic tissues24 given our sample size, and iPSC eQTLs showed similar genomic properties to eQTLs in cell lines and tissues (Extended Data Fig. 7; Supplementary Table 5).
As many eQTLs are shared among tissues9,25,26, we sought to place iPSC eQTLs in the broader context of somatic tissues. We assessed iPSC eQTL replication across 44 GTEx tissues (lead eQTLs and proxy variants, r 2 > 0.8, defining replication as P < 0.01/45; Methods), revealing 2,131 eQTLs that were specific to iPSCs (Fig. 4a). We also considered secondary eQTLs, identifying a similar proportion of iPSC-specific genetic effects (both 32%). Most tissue-specific signals (72%) occurred in genes with at least one GTEx eQTL that was not in high linkage disequilibrium (LD) with the lead iPSC eQTL variant, suggesting that iPSC-specific eQTLs are frequently driven by alternative regulatory variants. Only 11% of the iPSC-specific eQTLs could be attributed to tissue-specific gene expression (Fig. 4b), despite greater numbers of expressed genes in iPSCs compared with somatic tissues (Extended Data Fig. 7). Similarly, most somatic tissue-specific eQTLs were also driven by alternative regulatory variants, with only testis showing a substantial fraction (16%) of eQTLs attributable to tissue-specific gene expression (Fig. 4b). Using alternative methods for eQTL detection and assessing the extent of sharing between eQTLs in iPSCs and GTEx tissues, we confirmed that our conclusions were robust to methodological differences between GTEx and our study (Extended Data Fig. 8). However, due to variation in sample size, we cannot rule out that a fraction of the iPSC-specific eQTLs may have a weak effect on gene expression in some somatic tissues.
The transcriptional regulatory networks that maintain pluripotency are unique to stem cells. We next investigated how common genetic variants modulate these networks to produce iPSC-specific genetic effects on expression. We used chromatin state annotations from 127 reference epigenomes from the Roadmap Epigenomics Project27 to quantify the fold enrichment of iPSC-specific and nonspecific eQTL sets across 25 chromatin states (using matched null variants; Methods). iPSC-specific eQTLs were enriched in active enhancers and poised promoters in PSCs and PSC-derived cell types, while shared eQTLs were enriched for active promoters and transcribed regions in somatic tissues (Fig. 4c). iPSCspecific eQTLs were also enriched for binding sites of NANOG, POU5F1 (OCT4), and multiple other pluripotency factors9,28 ( Fig. 4d; Extended Data Fig. 8). Our results suggest that common genetic differences between individuals may affect expression regulation during early stages of development.

iPSC eQTLs tag common disease variants
We next identified iPSC eQTLs that may be associated with disease. iPSC eQTLs tagged 322 variants associated in genome-wide association studies with 145 different disease traits, corresponding to a 1.4-fold global enrichment over control variants (Fisher's exact P = 1.4x10 -6 ), and trait-specific enrichments for seven traits (Supplementary Table 6), a comparable level of enrichment to eQTLs from most somatic tissues (Extended Data Fig. 9). We also observed that iPSC eQTLs tagged a larger number of known cancer genes (COSMIC cancer census 27/04/201620) than somatic tissue eQTLs, with only cancer eQTLs tagging more (Extended Data Fig. 9).
Next, we used statistical colocalisation29 to identify loci where the same causal variant appeared to be driving both an iPSC eQTL and an association with one of 14 complex traits, identifying 233 loci where the posterior probability of a joint association exceeded 0.5 (Supplementary Table 6). Of these, 45 were iPSC-specific, including PTPN2, an iPSCspecific eQTL that strongly colocalised with risk variants for four autoimmune disorders (Fig. 5a). Previous eQTL studies in both immune cells30-32 and GTEx tissues have not identified a PTPN2 eQTL (Extended Data Fig. 9), suggesting that disease risk variants at PTPN2 may function in stem cells, or early development.
Statistical colocalisation analysis is limited to instances where full summary statistics are available for both traits. For other disease traits available in the GWAS catalogue, we searched for sharing of lead iPSC-specific eQTL and GWAS SNPs. We found six variants where the lead eQTL variant was identical to a catalogued GWAS variant, with no other common variants in LD (r 2 < 0.8). One example was rs10069690, the lead eQTL variant for the TERT (Telomerase Reverse Transcriptase) gene (Fig. 5b). Although this variant is associated with germline predisposition to seven cancers33-35, this eQTL is not reported in cancer eQTL studies36-38 nor in any GTEx tissue. Previous studies have reported aberrant splicing of TERT caused by rs1006969039. We quantified TERT intron retention rates and found that the minor allele of rs10069690 increased the fraction of TERT transcripts in which intron four is retained (P = 1.7x10 -9 , Bonferroni adjusted) (Fig. 5d, Extended Data Fig. 10). Somatic TERT promoter mutations only manifest in differentiated cells, resulting in increased telomerase activity40. We speculate that the germline TERT eQTL we identified in iPSCs results in genotype-dependent variability in telomerase activity in somatic cell types, leading to differential cancer susceptibility.

Discussion
Here we present the most comprehensive analysis yet of genetic and phenotypic data from human iPSC lines. Our study substantially extends on previous work6,7 by demonstrating widespread functional consequences of genetic variation for many molecular and cellular phenotypes in human pluripotent stem cell lines, including in the efficiency with which iPS cells differentiate41-43. This is potentially a consequence of variation in core components of the regulatory networks controlling cellular differentiation and responses to external environmental stimuli, as observed previously in hematopoietic cells and mouse and fly embryos44-46.
We have also created a high-resolution map of recurrent genetic abnormalities in hiPSCs and identified plausible candidate targets of selection. The majority of these recurrent loci are rare and were not reliably identified in previous studies with smaller sample sizes. Compared to previous work11,12, we observed substantially lower levels of genetic aberrations. One possible explanation is that access to donor-matched reference samples helped us more accurately identify germline CNAs that would otherwise have inflated our estimates, while previous studies in ESCs were unable to perform similar comparisons.
Our study provides the highest resolution map to date of common regulatory variation in human PSCs. We show that variation in local gene regulation in iPSCs is similar to that in somatic tissues, with eQTLs driving cell-type specific expression profiles through distal tissue-specific regulatory elements. We have identified eQTLs that function primarily in pluripotent cells, a subset of which tag loci associated with disease. These loci may drive disease-susceptibility through molecular changes early in development or, more generally, in cells with 'stem-like' characteristics, which are not well captured by studies of differentiated primary tissues from adult individuals. A compelling example of this is the iPSC-specific eQTL regulating TERT expression. In human tissues, telomerase activity is mainly restricted to stem cells, with most somatic tissues silencing TERT expression. However, cancer cells bypass this tumour suppressive mechanism by reactivating telomerase activity47. This result highlights how iPSCs could be used to study the genetic effects of diseases that manifest in transient states during cellular growth and differentiation, including in cancer48.
The analysis of the recurrence of CNAs and the eQTL map we present are based on a large sample, providing a high-confidence map of molecular associations in iPSCs. We have presented preliminary experimental characterisation of some of the CNAs and eQTLs we detected, however our results are inconclusive. An important next step will be to perform more extensive functional characterisation to understand how iPSC cellular phenotypes are influenced by CNAs and iPSC eQTLs. We anticipate that the lines and data we have generated here will be a valuable starting point for future studies to understand how germline and somatic genetic variation influences iPSC growth and differentiation.
In summary, our study provides a detailed picture of the genetic and phenotypic variability in human pluripotent stem cells, including the major drivers of this variation. Data and cell lines from this study are being made available through www.hipsci.org, the European Collection of Authenticated Cell Cultures (ECACC) and the European Bank for Induced Pluripotent Stem Cells (EBiSC). As the HipSci resource continues to expand in sample size and assays, it will enable the study of subtler genetic effects, under a wider range of conditions, in an increasing range of disease-relevant differentiated cell types.

Generation of iPSC lines
All samples for the HipSci resource were collected from consented research volunteers recruited from the NIHR Cambridge BioResource Fibroblast isolation-Primary fibroblasts were derived from 2 mm skin punch biopsies from each donor. Biopsies were collected in fibroblast growth medium (Advanced DMEM, 10% FBS, 1% L-Glutamine, 0.007% 2-mercaptoethanol and 1% Pen/Strep) in falcon tubes at room temperature. Biopsies were manually dissected using a microscope under a drop of fibroblast medium using sterile scalpels. The biopsy fragments were transferred onto a 60 mm Petri dish containing several drops of fibroblast growth medium. Sterile cover slips were placed onto the dissected pieces of tissue to hold them in place against the bottom of the plate. The explants were cultured for five days and the spent media was removed and replaced with a few drops of media (1 ml) to prevent dehydration. The explants were fed every five days with 1 ml fibroblast media until fibroblast outgrowths appeared. The explants were screened for presence of mycoplasma using a standard PCR kit (EZ-PCR Kit, Gene flow (41106313-001)). On average outgrowths appeared within 14 days, with a small fraction of samples failing to produce outgrowths (12% of cases). Failures were due to contamination (0.5%) or lack of observed outgrowths after 30 days (11%). Approximately 30 days post dissection, when the fibroblasts had reached confluence, the culture was trypsinized and passaged into a 25 cm 2 tissue culture flask. When 80-90% confluent, the fibroblasts were further passaged into a 75 cm 2 flask. Cells were then expanded to confluency in 225 cm 2 flasks (at a split ratio of 1:3) and either cryopreserved at 1-2 million cells per vial in FBS and 10% DMSO or seeded immediately for reprogramming. iPSC derivation-Fibroblasts were transduced in one well of a six-well plate using Sendai vectors expressing hOCT3/4, hSOX2, hKLF4, and hc-MYC 50 (CytoTune™, Life Technologies, Cat. no. A1377801). The transduced cells were cultured on an irradiated mouse embryonic fibroblast (MEF-CF1) feeder layer on a 10 cm 2 tissue culture dish in iPSC medium consisting of Advanced DMEM (Life technologies, UK) supplemented with 10% Knockout Serum Replacement (KOSR, Life technologies, UK), 2 mM L-glutamine (Life technologies, UK) 0.007% 2-mercaptoethanol (Sigma-Aldrich, UK), 4 ng/mL of recombinant Zebrafish Fibroblast Growth Factor-2 (CSCR, University of Cambridge), and 1% Pen/Strep (Life technologies, UK). Cells with an iPSC morphology appeared approximately 25 to 30 days post-transduction. The undifferentiated colonies (six per donor) were picked between days 30-40, transferred onto 12-well MEF-CF1 feeder plates and cultured in iPSC medium with daily media change until ready to passage. Cells were passaged every five to seven days, depending on the confluence and morphology of the cells, at a maximum 1:3 split ratio until established -usually at passage five or six. Once the iPSC lines were established in culture, three of the six lines were selected based on morphological qualities (undifferentiated, roundness and compactness of colonies) and expanded for banking and characterisation.
Transfer to feeder-free culture-Between passages four to eight, selected feederdependent iPSC lines were transferred to feeder-free culture. The feeder-dependent iPSC lines were split and passaged onto both feeder-dependent and feeder-free conditions. The feeder dependent lines continued to be cultured on MEF-CF1 feeder plates in iPSC medium, whilst the feeder free lines were cultured in Essential 8 (E8) medium on tissue culture dishes coated with 10 µg/ml Vitronectin XF (StemCell Technologies, UK, 07180). E8 complete medium consists of basal medium DMEM/F-12(HAM) 1:1(Life technologies, UK, A1517001) supplemented with E8 supplement (50X) (Life technologies, UK, A1517001) and 1% Pen/Strep (Life technologies, UK, 15140122). Media was changed daily. To passage feeder-free iPSC lines, cells were washed with PBS and incubated with PBS-EDTA solution (0.5 mM) for 5-8 minutes. PBS-EDTA solution was removed, and cells were resuspended in E8 medium and seeded at split ratios ranging from 1:3 to 1:6 onto Vitronectin coated tissue culture dishes. Cells were passaged every four to seven days (depending on the confluence and morphology of the cells). Once the feeder-free iPSC lines were established in culture, the cells were expanded for banking and characterisation.
iPSC line selection and molecular assays-Each iPSC line was passaged on average 16 times before being expanded for the collection of initial molecular data for quality control ('Tier 1 assays'). These included genotyping ('gtarray'), gene expression data ('gexarray'), and an assessment of the pluripotency and differentiation potential of each line ('Cellomics'). Pluripotency of the lines was additionally verified in silico, using the PluriTest assay 9. Following Tier 1 assays, one or two lines were selected from a subset of donors (hereafter 'selected lines') and further expanded to enable collection of a richer set of molecular data ('Tier 2 assays'). The criteria for line selection were: (i) level of pluripotency, as determined by the PluriTest assay (ii) number of copy number abnormalities and (iii) ability to differentiate into each of the three germ layers. These included proteomics, DNA methylation ('mtarray'), RNA-sequencing and high-content cellular imaging. Once Tier 1 genotyping data were collected (see below) cell lines originating from the same donor were checked for possible sample swaps using BCFtools (bcftools gtcheck -G1).

Nucleic acid extraction
DNA and RNA from iPSC lines were extracted using Qiagen Chemistry on a QIAcube automated extraction platform. Sample volume was checked using a BioMicroLab automated volume check system. The PicoGreen assay was used to measure the concentration of the samples, using both Beckman FX liquid handling platforms and Molecular Devices plate readers. Invitrogen E-Gels were run to check sample integrity; the loading of these gels was automated using Beckman FX/NX liquid handling platforms. A standard Fluidigm genotyping assay containing 24 SNPs (22 autosomal and 2 gender markers) was performed to produce a fingerprint of the samples, which was used to confirm sample identity after sequencing or genotyping. The gender markers also allowed for sample swaps and plate orientation issues to be identified prior to downstream analysis. Samples that passed quality control were quantified to 50 ng/µl by the onsite sample management team prior to submission for sequencing.

Genotyping ('gtarray')
Experimental processing of arrays-Samples were hybridised to the Illumina HumanCoreExome-12 Beadchip according to the manufacturer's guidelines. Four microlitres (200 ng) of DNA is required for the pre-amplification reaction using a Tecan Freedom Evo. The process is automated except for manual agitation/centrifugation step midway through and at the end of the process. Post-amplification processes (fragmentation, precipitation, resuspension, hybridization to beadchip and xStaining) were completed over three days as per Illumina protocol. Following the staining process, beadchips were coated for protection and dried completely under vacuum before scanning on the Illumina iScan paired with Illumina Autoloader 2.x. Prior to downstream analysis, all samples were subjected to initial quality control to establish that the assay was successful. Sample callrates below 92.5% were flagged before loading samples into Illumina's GenomeStudio software. Using Illumina's QC dashboard, sample performance was assessed by measuring dependant and non-dependent controls that are manufactured onto each beadchip during production.
Genotype calling and imputation-After primary quality control, the Genotyping (GT) module of the GenomeStudio software (Illumina, CA, USA) was used to call the genotypes. For each probe, the GT module estimates the Log R ratio and B-allele frequency for each sample using a clustering model applied to the distribution of signal intensities. These statistics are used internally by GenomeStudio to assign the sample genotypes for each marker. Variant coverage was further increased using statistical imputation and phasing. We constructed a reference panel of haplotypes from a combination of SNPs and small insertions and deletions (indels) in the UK10K cohorts and 1000 Genomes Phase 1 data 51,52. Samples were independently imputed using IMPUTE2 v2.3.1 53 and subsequently phased using SHAPEIT v2.r790 54. This analysis was done in chunks of average 5 Mb, with 300 Kb buffer regions on each side. IMPUTE2 was used with its default MCMC options (-Ne 20000 -k 80) for autosomes and -Ne 15000 -k 100 for X chromosome. SHAPEIT was run without MCMC iteration (-no-mcmc) so that each sample was phased independently using the reference panel as the haplotype scaffold regardless of the phasing of the other samples. Single-sample VCFs were merged together and INFO scores were re-calculated from genotype posterior probabilities (GPs). Variants with INFO score less than 0.4 were excluded from further analysis. Cell lines originating from the same donor were checked for possible sample swaps using BCFtools (bcftools gtcheck -G1). Swapped samples typically had large number of discordant genotypes (>20%), whereas in samples from the same donor the number of discordant genotypes was low (<0.3%), even in the presence of large copy number variation.

Gene expression arrays ('gexarray')
Experimental processing of arrays-500 ng of total RNA for each sample was amplified and purified using the Illumina TotalPrep-96 RNA Amplification kit (Life Technologies, UK), according to the manufacturer's instructions. Biotin-Labelled cRNA was then normalized to a concentration of 150 ng/ul and 750 ng was hybridised to Illumina Human-12 v4 BeadChips (Illumina, CA, USA) for 16 hours (overnight) at 58°C. Following hybridisation, BeadChips were washed and stained with streptavidin-Cy3 (GE Healthcare, UK). BeadChips were then scanned using the BeadArray reader and image data was then processed using GenomeStudio software (Illumina, CA, USA).
Re-mapping of array probes-Prior to analysis, array probe sequences (N probes = 47,230; length 50 bp) were re-mapped against the human genome build 37 using BWA version 0.7.5 55.
We first mapped the sequences allowing no mismatches (-n 0, seeding disabled) and kept uniquely mapping probes with a minimum mapping quality (MAPQ) of 10 (-q 10). These sequences were then mapped again, this time allowing one mismatch (-n 1). Again, only uniquely mapping probes with MAPQ > 10 were retained, resulting in a total of N = 37,740 probes. We further removed all probes that overlapped with any variant with a minor allele frequency greater than 0.05 in the main imputed dataset (N lines = 858). Remaining probes were annotated with Gencode version 19 gene annotations 56 and only probes mapping uniquely to a single gene were kept (final probeset N probes = 25,604, representing 17,116 unique genes of which 14,569 are protein coding. Pre-processing and normalization of data-Gene expression profiles were measured with Illumina HumanHT-12 v4 Expression BeadChips. After limiting the dataset to iPSC lines derived from fibroblast of healthy donors, we obtained data from 711 iPSC lines and 301 somatic fibroblast lines. Probe intensity estimates were normalised separately for the two cell types using the variance-stabilizing transformation implemented in the R/ Bioconductor vsn package 57. After normalization, the datasets were limited to the final remapped set of probes (N probes = 25,604). We refer to this version of the "gexarray" data by vsn log2 (iPSC/somatic). Kilpinen

Cellular differentiation assay ('Cellomics')
Differentiation potential assay-Selected iPSC lines were assessed for their pluripotency and differentiation properties by culturing the cells under conditions favouring the formation of the three embryonic germ layers, and subsequent immunostaining with markers specific for pluripotency and differentiation. Differentiation was performed as described previously 58. Briefly, iPSCs grown in feeder-dependent or feeder-free conditions were harvested using either collagenase and dispase or EDTA, respectively. Colonies were collected, washed in media and mechanically broken up before being re-plated onto 24-well mouse embryonic fibroblast (MEF) feeder plates or pre-coated gelatine/FBS plates. For pluripotency assays, feeder-dependent colonies were seeded on MEF feeder plates and feeder-free colonies onto Vitronectin plates. For the differentiation assay, colonies were grown on gelatine/FBS plates. Prior to differentiation to mesoderm (dME), endoderm (dEN), and neuroectoderm (dEC), cells were cultured overnight in pre-differentiation media CDM-PVA supplemented with recombinant Activin-A (10 ng/ml; CSCR, University of Cambridge) and zebrafish FGF2 (12 ng/ml; CSCR, University of Cambridge).
Immunostaining for pluripotency and differentiation markers-For the detection of pluripotency and differentiation markers, cells grown in 24-well plates were fixed with 4% paraformaldehyde for 20 minutes. Cells were permeabilized and blocked with 10% donkey serum and 0.1% Triton X-100 in PBS. Subsequently, cells were stained with primary antibodies overnight at 4°C and finally incubated with fluorochrome-labeled secondary antibodies (Invitrogen, UK). The primary antibodies used for detecting pluripotency markers were: anti-OCT4 (SC-5279, Santa Cruz Biotech, USA), anti-SOX2 (AF2018, R&D, UK), Kilpinen  Additionally, DAPI staining was used to label cell nucleus in order to facilitate cell segmentation.
Images were captured and quantified using a Cellomics Array Scan imaging system. Briefly, images were taken in 24-well plates. Individual plates were used to either measure pluripotency markers or markers to assess differentiation for one of the germ layers. Each plate contained cells from one or two cell lines, as well as technical replicates for each measurement. Three types of plate layouts were considered throughout the project: Twochannel, three-channel, and three-channel with single staining. For all layouts, the signal from the DAPI staining was read in the first channel. The first columns of each plate were used for marker staining; subsequent columns (one or two) were stained with the secondary antibody to measure background signal (Extended Data Fig. 1).
Processing of images on the Cellomics instrument-Individual wells in the plate were imaged consecutively, either until the whole plate was imaged or until 10,000 individual cells were detected. Cell detection was performed based on nucleus segmentation from DAPI staining. All considered markers except for CXCR4 are nuclear markers, so their signal intensities were measured in the segmented nucleus area. The cell surface marker CXCR4 was quantified in a circle around the segmented nucleus. For each cell and marker, we used the average intensity within the respective quantification area as final readout. Each batch of lines for staining included the reference line ('CTRL0214pf-iely'). This reference line was used to determine parameter values for cell size (usually around 30-400) and an approximate intensity threshold for detecting responding cells.
To quantify Cellomics phenotypes, we fit a Gamma mixture model to the Cellomics raw intensities ( Supplementary Information). Briefly, this model was fit to primary wells as well as background wells (Extended Data Fig. 1), thereby estimating both the proportion of responding cells as well as the overall intensity (expression) of the corresponding cells.
For downstream analyses, technical replicates on each plate were aggregated using average values. Analogously to the processing steps for gene expression arrays, we regressed out batch (derived from the date of staining), media type, gender, passage number, plating technician, fixation technician, and the technician in charge of the staining. Analyses of proportions of responding cells are based on estimates of the proportion of responding cells (Fig. 1c,d). Analyses that consider intensities were based on quantitative expression estimates, averaged across individual markers for a given layer (Extended Data Fig. 1c).

Proteomics
Sample preparation-Frozen iPSC pellets were thawed and washed with PBS twice prior to lysis. The protein content of the cells was extracted by re-dissolving the pellets in 8 M urea, 100 mM TEAB, pH 8.5 and mixing at room temperature for 15 minutes. Next, the DNA content of the cells was sheared using ultrasonication. The protein amount was determined using a fluorescence based assay (EZQ, Life Technologies) prior to double digestion using mass spectrometry grade lysyl endopeptidase (Wako, Japan) and trypsin (Pierce) in a substrate-to-enzyme ratio of 1:50; w:w, at a final urea concentrations of 2 M and 0.8 M, respectively. The digested proteins were desalted using sepak vacuum cartridges (waters) and dried in vacuo. The desalted peptides were redissolved in (10 mM borate at pH 9.3 : acetonitrile; 80:20; v:v) for hSAX fractionation using a 40 minute gradient. A total of 16 fractions were collected, desalted and dried. The hSAX fractions were redissolved in 5% formic acid for label-free LC-MS analysis. In addition to individual samples, a composite reference sample ('HPSI_composite_1503') was constructed by pooling together protein lysates from 43 iPSC lines. 2 mg of protein was used for each. All samples in this reference were of fibroblast origin and reprogrammed with sendai virus.
For Tandem Mass Tag (TMT)-based quantification, the dried peptides were re-dissolved in 100mM TEAB (50 µL) and their concentration was measured using a fluorescent assay (CBQCA) (Life Technologies). 100 µg of peptides from each cell line to be compared, in 100 µL of TEAB, were labelled with a different TMT tag (20 µg ml -1 in 40 µL acetonitrile) (Thermo Scientific), for two hours at room temperature. After incubation, the labelling reaction was quenched using 8 µl of 5% hydroxylamine (Pierce) for 30 minutes and the different cell lines/tags were mixed and dried in vacuo.
The TMT samples were fractionated using off-line high pH reverse phase chromatography: samples were loaded onto a 4.6 x 250 mm Xbridge™ BEH130 C18 column with 3.5 µm particles (Waters). Using a Dionex bioRS system, the samples were separated using a 25minute multistep gradient of solvents A (10 mM formate at pH 9) and B (10 mM ammonium formate pH 9 in 80% acetonitrile), at a flow rate of 1 ml/min. Peptides were separated into 48 fractions, which were consolidated into 24 fractions. The fractions were subsequently dried and the peptides re-dissolved in 5% formic acid and analysed by LC-MS.

LC-MS/MS
Label-free analysis: RP-LC was performed using a Dionex RSLC nano HPLC (Thermo Scientific). Peptides were injected onto a 75 μm × 2 cm PepMap-C18 pre-column and resolved on a 75 μm × 50 cm RP-C18 EASY-Spray temperature controlled integrated column-emitter (Thermo) using a four-hour multistep gradient from 5% B to 35% B with a constant flow of 200 nL min-1 as described previously 59,60. The mobile phases were: 2% ACN incorporating 0.1% FA (Solvent A) and 80% ACN incorporating 0.1% FA (Solvent B). The spray was initiated by applying 2.5 kV to the EASY-Spray emitter and the data were acquired on a Q-Exactive Orbitrap (Thermo Scientific) under the control of Xcalibur software in a data dependent mode selecting the 15 most intense ions for HCD-MS/MS. TMT-based analysis: 5% of the material was analysed using an orbitrap fusion tribrid mass spectrometer (Thermo Scientific), equipped with a Dionex ultra high-pressure liquid chromatography system (nano RSLC). RP-LC was performed using a Dionex RSLC nano HPLC (Thermo Scientific). Peptides were injected onto a 75 μm × 2 cm PepMap-C18 precolumn and resolved on a 75 μm × 50 cm RP-C18 EASY-Spray temperature controlled integrated column-emitter (Thermo), using a four-hour multistep gradient from 5% B to 35% During the MS3 analysis for more accurate TMT quantifications, 5 fragment ions were coisolated using synchronous precursor selection using a window of 2 Th and further fragmented using HCD collision energy of 55%. The fragments were then analysed in the orbitrap with a resolution of 60,000. The AGC target was set to 1.0 e5 and the maximum injection time was set to 105 ms.

RNA-sequencing
Library preparation and sequencing-mRNA in total RNA was isolated and converted into non-stranded or stranded libraries. Non-stranded libraries were produced manually using reagents provided in the Illumina TruSeq RNA Sample Preparation Kit v2 in accordance with manufacturer's recommendations. The protocol was modified to produce size-selected libraries by modifying the fragmentation conditions and using a Caliper LabChip XT instrument. Stranded libraries were prepared using a NeoPrep Library Prep System and the reagents provided in the Illumina TrueSeq Stranded mRNA Library Preparation kit. The stranded library prep workflow is similar to the non-stranded workflow, except that it involves additional ribosomal reduction chemistry to maximise the percentage of uniquely mapped reads. Following purification, the RNA was fragmented and synthesised into cDNA using a reverse transcriptase process. The products were then enriched with PCR (maximum 10 cycles) to create the final cDNA library. Enriched libraries were subjected to 75 base paired-end sequencing using Illumina HiSeq 2000 v3 kits following manufacturer's instructions. Pre-processing of sequence data-Raw RNA-seq reads were aligned using STAR v2.4.0 63 against the 1000 Genomes Phase2 reference genome assembly that integrates the GRCh37 primary assembly with the human decoy sequence 37d5. Exon-intron junctions derived from Gencode v19 transcript annotations 56 were used to improve the alignments.
The same approach was taken to re-align data for two tissues from the GTEx Project (Extended Data Fig. 8). Raw fastq-files were obtained from dbGap (accession phs0004242.v6.p1.c1) for Adrenal Gland (N=126 samples) and Esophagus Gastroesophageal Junction (N=127 samples; limited to 126 unique samples used in the GTEx V6p map).
Quantification and normalization of data-Mapped reads were quantified on the level of genes using HTSeq version 0.6.1p1 64 and annotations from Gencode v19 56. We used the 'union' method of 'htseq-count' for unstranded libraries (-s no) and considered only uniquely mapping reads (-a 255; with 255 indicating uniquely mapped reads from the STAR aligner). Of note, STAR only outputs properly paired reads. Raw gene counts were scaled across individuals with scaling factors obtained with DESeq 65. The same approach was used to generate 'probe-level' counts using the final re-mapped set of gexarray probes. These were used to filter expressed probes in the CNA analysis. Finally, an alternative set of gene-level quantifications was generated for quality control purposes using RNA-SeQC v1.1.8 66. To match the original GTEx v6p quantifications as closely as possible, we ran RNA-SeQC with the -strictMode flag and used custom exon annotations generated and used by GTEx (gencode.v19.genes.v6p_model.patched_contigs.gtf.gz) to obtain gene RPKMs. The same RNA-SeQC quantification pipeline was applied to the two re-mapped GTEx tissues.

High-content cellular imaging
Sample preparation and cellular imaging-Each line was cultured and plated as previously described 67. Briefly, 96-well plates were coated with three concentrations of Fibronectin in alternated columns in a randomised fashion. Cell lines were seeded also in rows in a randomised fashion. 3,000 cells were plated and fixed after 24 hours. EdU was incorporated 30 minutes before fixation. Plates were then fixed and stained with DAPI and cell mask and EdU staining. Images were acquired using the Operetta (Perkin Elmer) high content device. Using the Harmony software, measurements were derived for each cell. Measurements included intensity features (DAPI, EdU), morphology features (cell area, cell roundness, cell width to length ratio, nucleus area, nucleus roundness, nucleus width to length ratio) and context features related with cell adhesion properties (number of cells per clump). Processing quantification and normalization of data was performed as previously described 67. The Cellomics fluorescence imaging data was used to quantify the fraction of cells expressing each protein marker independently. In the absence of co-staining information we sought to use the marginal fractions of cells expressing each marker to calculate lower and upper bounds for the fractions of cell estimating all markers simultaneously. Let be the fraction of cells expressing protein marker and the fraction of cells expressing all n markers simultaneously, then this value is bounded by: Variance component analysis-Feature (gene, protein, or probe) intensity estimates for each of the assays were pre-processed and normalised as described in the individual assay sections and subsequently transformed into a standard normal distribution across lines. For each feature in each assay, variance was partitioned using a linear mixed model (implemented in the lme4 R package) fitted with all metadata variables as random effects.
Only lines with complete metadata information were included in each one of the analyses, these numbers are shown in parenthesis in Fig. 3a-c. The variance components were normalized to sum to one and subsequently averaged across the different sets of features considered (Fig. 3). The fraction of non-technical variance explained by each biological or experimental factor refers to the variance explained by the factor, divided by the total variance minus the variance explained by assay batches (see below for a definition of experimental and assay batch factors). Confidence intervals for the Cellomics variance components were obtained with the 'profile' method implemented in the 'confint' function of the lme4 package. The following random effects were included for each assay: • Methylation -Donor. Experimental factors (summed up to produce Fig. 3a): gender, passage interval at time of assay. Assay batches (summed up to produce Fig. 3a): sentrix id, sentrix array (within sentrix id), bisulfite conversion plate, year and month of assay, year and month of Tier 1 assays.

CNA analysis
Pairwise fibroblast-iPSC CNA detection-Copy number differences between fibroblast and iPSC lines from the same donor were checked using a HMM algorithm implemented in BCFtools/cnv for this purpose 10. In order to distinguish between normal and novel copy number variation as well as to reduce the number of false calls, the program was run in the pairwise mode (bcftools cnv -c <donor> -s <derived>) with default parameters. The CNA calls were filtered to exclude calls with quality score smaller than 2, deletions with fewer than 10 markers, and duplications with fewer than 10 heterozygous markers. Three sets of CNA calls were generated: a more lenient set containing all calls >= 0.2 Mb in length, a set with all calls >= 0.5 Mb and a stricter subset of the previous with calls >= 1 Mb.
Statistical significance of recurrent CNAs was estimated from the complementary cumulative distribution function of the binomial distribution, and the significance of subchromosomal events was estimated using a permutation test (Supplementary Information).
Overlap with annotated regions-To assess the significance of the overlap between CNAs and annotated regions (namely chromatin fragile sites 16 and recurrent somatic copy number altered regions in cancer 68 we randomly generated a set of 2,000 matched control regions for each CNA. Control regions were generated so that they had the same size as the CNA, did not overlap with telomeres or centromeres, and did not overlap the original CNA.
Overlaps were determined between the CNA and the annotated regions and between the matched control regions and the annotated regions to calculate an empirical P-value.
Association between CNAs and gene expression-To determine the functional consequence of CNAs with regards to gene expression we first selected CNAs for which five (~1%) or more of the cell lines had a copy number different from two (regardless of copy number in the corresponding somatic cells). For each of these CNAs we took the copy number at the region of peak coverage for each cell line (see CNA coverage plots in Fig. 2c and Extended Data Fig. 4). We then defined the set of expression array probes to test by choosing only expressed probes. Here a probe was defined as expressed if the number of RNA-seq fragment counts (normalised between samples for sequencing depth) mapping to the genomic regions targeted by the probe is greater than 0 in 10% or more of the lines. Finally, we used a linear mixed model to independently test for association between copy number of each CNA and the intensity of each probe. We included culture condition, gender and an interaction between copy number and culture condition as fixed effects; and used donor and assay batch as random effects. Q-values were obtained for each CNA using Benjamini Hochberg to adjust for multiple testing. The same approach was employed to test for association between X chromosome copy number and gene expression, but we limited the tests to female samples and to probes on the X chromosome.
Gene set enrichment analysis-Pathway enrichment analysis of the genes regulated by chr17 was performed with GSEA 69 on the full list of genes ordered by effects size of association between gene expression with copy number and with a custom set of pathways. The custom set of pathways considered comprised 1,156 super pathways from the PathCards Europe PMC Funders Author Manuscripts database 23 filtered to exclude pathways relating to infectious diseases and pharmacokinetics. Multiple testing correction was performed as described in 69.

Expression quantitative trait loci (eQTL) eQTL mapping
Pre-processing of genotype data: Variants in the original VCF files were renamed to format "type_chr_pos" (e.g. 'snp_1_236887241' or 'indel:2D_1_18847945'), and filtered for polymorphic and bi-allelic sites with VCFtools v.0. Pre-processing of expression data: 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 prenormalized 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).

Linear mixed model:
Expression quantitative trait loci (eQTL) were identified using a linear mixed model implemented in LIMIX 73,74. eQTLs were mapped in cis, considering a window of 1Mb around the gene start (as defined by Gencode v19 annotations). We modelled the genotype as a fixed effect, with population structure included as random effect. Population structure was modelled with a kinship matrix, calculated as the dot product of the genotypes in trans for each cis window (realized relationship).
eQTL mapping was performed with the following datasets (Supplementary Table 4

Europe PMC Funders Author Manuscripts
When multiple lines were available for a donor, the mean expression value of lines was used.
In the array-based iPSC dataset, which had the largest number of replicate lines available per donor (246 donors with multiple lines), we additionally mapped the eQTLs using two randomly drawn sets of individual lines per donor to assess the replicability of the iPSC eQTLs (Extended Data Fig. 6f-h). In the main RNA-seq based map, we identified both primary and secondary eQTL effects. To identify secondary effects, we repeated the mapping with the genotypes of the lead eQTL variant included as a covariate in the model. Finally, an alternative version of the main eQTL map was generated using a pipeline matched with GTEx V6p eQTLs for quality control purposes (See Supplementary  Information).

Multiple testing correction:
For cis eQTLs (primary and secondary), to adjust for multiple testing, we permuted genotype sample labels in each cis window 10,000 times, keeping everything else in the model constant. To derive an empirical P-value distribution, the test statistic of the most significant variant in each permutation round was stored. A region-wise adjusted cis P-value was derived from the proportion of permuted test statistics that were larger than the most significant observed test statistic in the region. These threshold P-values were further adjusted for genome-wide analysis using the Benjamini-Hochberg (BH) correction. A gene was considered an eGene if its final BH-adjusted P-value was less than 0.05.  Supplementary Table 5a. For each discovery tissue, we tested for the replication of all lead eQTL effects (lead eQTL variant -target eGene; hereon referred to as 'ePair') originally reported (FDR 5% in both HipSci and GTEx; Supplementary Table 5a). Ideally, the exact same ePair would have been tested across all tissues. However, due to differences in genotyping methods and allele-frequency due to sample size it was not always possible to query the exact same ePair in all tissues. To account for this, if the original lead variant was not available to test in the query tissue, a proxy variant was tested instead. We note that the selected proxy may differ across the replication tissues. The approach to define the proxy variant is summarized in the Supplementary Information.

Identification of tissue-specific cis eQTLs-In
Replication was defined as the query variant (original lead or proxy) having a nominal eQTL P < 2.2x10 -04 for the same eGene (corresponding to P = 0.01 / 45, where 45 refers to the total number of tissues tested). A lenient threshold was chosen in order to rule out any evidence of replication. High-LD proxies for a lead were defined as having r 2 > 0.8 in the UK10K European reference panel and located in the same cis window. All available LDproxies were tested for replication and the variant with the most significant eQTL P-value in the discovery tissue was stored as the proxy. If no LD-proxies could be tested, the cis variant with the most significant eQTL P-value in the discovery tissue overall was selected. Overall, the same lead variant was available to test in 95% of tests across tissues (median of discovery tissues, with each discovery tissue represented as the median of tests across all replication tissues; 90% for iPSC eQTLs). A high-LD proxy for the lead was tested 3.6% of the time (7.3% in iPSC), while the best available cis variant was tested only 1.1% of the time (2.7% in iPSC). Of the rare cases when the best available variant was tested, the selected variant was independent of the original eQTL effects (r 2 < 0.1) 0.2% of the time (1.2% in iPSC), indicating that in the vast majority of cases, the same eQTL effect was tested for replication and the choice of variant is unlikely to have a marked effect on the results. Statistics of the replication of the iPSC eQTLs, including numbers of selected proxy variants, are provided in Supplementary Table 5b. If a gene was not tested for eQTLs in the query tissue (for e.g. due to gene not expressed), replication information was defined as missing. A replication profile was derived for each eGene in the discovery tissue, indicating whether the lead eQTL effect replicated (yes | no) or could not be tested ('NA'). We then extracted eGenes, for which the lead eQTL effect did not show evidence of replication in any other tissue (P > 2.2x10 -04 ) or could not be tested (hereon referred to as 'tissue-specific eQTLs'). Of note, in this analysis, for an eQTL effect to replicate, it has to affect the expression of the same gene in both tissues. It the same variant is an eQTL for two different genes in the two tissues, it is not considered replicating. We also investigated the impact of the specific replication threshold, considering a threshold of P<0.01 and P<0.05 (Extended Data Fig. 8c). This analysis showed that the ability to replicate an eQTL signal in a second tissue is primarily determined by the sample size of the replication tissue and not the specific choice of threshold. This dependency needs to be taken into account when interpreting tissue-specific eQTL effects.
As an alternative strategy for comparing eQTLs among different tissues, we also calculated the π 1 statistic (π 1 = 1 -π 0 ; 75) for all pairs of tissues (Extended Data Fig. 8b) using the qvalue package in R. The π 1 statistic provides a global measure of similarity between a pair of tissues by estimating the proportion of eQTL signal discovered in one tissue that shows evidence of replication in a second tissue. For this analysis, for each discovery tissue, we queried all significant variants per eGene in the full cis output of all other tissues and estimated q 0 with the 'bootstrap' method of the qvalue package.

Functional annotation of eQTLs
Matched sets of variants: For all functional enrichment analyses, 100 matched sets of variants (hereon 'control variants') were used as the null. These sets were generated with SNPsnap 76 using unique lead eQTL variants from different datasets as the input. Variants were matched for minor allele frequency, number of SNPs in LD ('LD buddies'; r 2 > 0.5), distance to the nearest gene, and gene density, allowing for maximum deviation of +/-50% for each criterion. HLA SNPs (defined as falling between positions 25,000,000 and 35,000,000 on chromosome 6) were excluded from the analysis and all matched sets were non-overlapping with the input variants. 1000 Genomes Phase 3 European population was used as the genotype reference panel. Both target and control sets of variants were further expanded with their high-LD proxies (r 2 > 0.8) derived from the European UK10k reference panel. These expanded sets form the basis for all subsequent enrichment analyses.
Overlap with chromatin state annotations: Functional enrichment of eQTLs was assessed using chromatin state data from the Roadmap Epigenomics Project 27. The data comprised of 25 chromatin states derived from reference epigenomes from 127 cell types. We overlapped target eQTL variants separately with each chromatin state and cell type. We also overlapped 100 sets of control variants with the same annotations and derived an empirical P-value for each enrichment. This was defined as the number of control variant sets (N) that showed a higher overlap with the target annotation than the eQTL lead variants (P = N/100). The empirical P-values were further adjusted for the number of tests (25 states x 127 cell types) within each eQTL set using the Q-value 75 package in R. Annotations with a Q-value < 0.05 were considered significantly enriched. We tested two sets of eQTLs for enrichment: iPSC-specific eQTLs (N = 2,131) and non-specific eQTLs (N = 4,500).
For visualization purposes we aggregated the 127 cell types into five clusters, using k-means clustering. The number of clusters was chosen based on the number of different sample types annotated by Roadmap (primary cell, primary culture, primary tissue, cell line, ESC derived). A heatmap of the difference in fold enrichment between iPSC-specific and nonspecific eQTLs (DIFF = FE specific -FE nonspecific ) was generated with the pheatmap package in R (Fig. 4d) in order to assess how well our definition of iPSC-specific enriches for functional elements active in stem cells.
Overlap with transcription factor binding sites: Functional enrichment of eQTLs was additionally assessed using ChIP-seq based transcription factor (TF) binding sites (TFBS) from the ENCODE Project 49. Specifically, we used a set of proximal and distal TFBSs, where proximal is defined as a 2,000 bp window centered on Gencode v19 annotated transcription start sites (TSS). These sets were constructed by overlaying TF ChIP-seq peaks from all available cell types with DNase peaks. We limited our analysis to those peaks that overlapped with H1-ESC derived TF peaks. As with chromatin state annotations, we overlapped target and control eQTL variants with binding sites for each individual factor (all sites where the given factor is bound or co-bound). An empirical P-value for the enrichment was derived based on the control sets and adjusted for multiple testing using the Benjamini-Hochberg method. Factors with an adjusted P < 0.05 were considered enriched. Plotted in Fig. 4 are factors with at least 10 observed overlaps.
Overlap with the GWAS catalogue: The NHGRI-EBI GWAS catalogue was downloaded on 2016-04-18 (release 2016-04-10). Entries with missing positional or P-value information were removed and the positions of the remaining entries were converted to hg19 with the UCSC liftOver function in the "rtracklayer" R/Bioconductor package 77 and the 'hg38ToHg19' chain file. This resulted in 24,861 catalogue entries, corresponding to 18,446 unique variants of which 6,681 were significantly associated to a trait (P < 5x10 -08 ). However, many studies included in the complete catalogue are not reliable, so we parsed the sample size information in the catalogue and excluded studies that did not report effect sizes (odds ratio or regression coefficient), had a sample size below 1,000, or assayed fewer than 100,000 variants. Then, following the approach taken in 78, we further filtered the set of remaining associations to retain only traits that had at least six significantly associated variants, and kept all associations with P < 1x10 -06 for these traits. This approach yielded a filtered set of 9,562 associations for 6,059 variants with 358 diseases and traits. Each variant position was then parsed for trait overlap (all traits associated to the given variant). If a variant -trait association was reported multiple times (e.g. by different studies), the most significant association was kept.
We first tested whether eQTLs in iPSCs and somatic tissues (lead variants and their high-LD proxies, control variants as before) showed global enrichment in this final set of diseaseassociated variants compared with matched sets of variants (Extended Data Fig. 9b). Fold enrichments were derived from a comparison to 100 sets of matched controls (mean of control overlaps per tissue). For iPSCs, we additionally tested enrichment for individual traits, deriving an empirical P-value. Traits with an adjusted P < 0.05 (Benjamini-Hochberg) and minimum five observed overlaps were considered enriched (Supplementary Table 6c).
Lastly, we parsed each disease variant position for overlap with iPSC eQTLs, again considering lead variants and their high-LD proxies for each tissue as follows. We report all disease variants that were tagged by iPSC eQTLs (lead or proxy) and the eGene(s) regulated by the eQTL lead/proxy (Supplementary Table 6a). We additionally report a subset of these variants that are exact matches with iPSC lead eQTL variants and highlight the number of high-LD proxies (r 2 > 0.8) each variant has (Supplementary Table 6b).  Fig. 9a).

Splicing of TERT:
Alternative splicing of the TERT gene was analysed using Leafcutter 91, which focuses on introns and quantifies both known and novel alternative splicing events by quantifying reads mapping to exon-exon junctions. Annotations were derived from Gencode v19. Introns supported by fewer than 30 reads (-m 30; default) across all samples were removed. We obtained quantifications for eight intron clusters within TERT. After removing individual introns with a mean intron usage of zero, we had a total of 22 introns to test. We used the intron excision proportions to assess genotype-dependent effect of rs10069690 on TERT splicing (linear model between genotype and excision proportion, Bonferroni correction of P-values for the total number of introns tested). One intron showed evidence of a splicing-QTL effect (P < 0.05, Bonferroni adjusted; Extended Data Fig. 10).

Data availability
The assay data used in this publication are listed in the Biostudies archive (https:// www.ebi.ac.uk/biostudies/studies) with accession identifier S-BSMS5. All data can be accessed via the HipSci data portal (http://www.hipsci.org), which references to EMBL-EBI archives that are used to store the HipSci data.

Code availability
Scripts that were used to process the raw data and for implementing the statistical analyses presented are available from https://github.com/hipsci/Nature2017.

Extended Data
Extended Data Figure 1. Overview of the Cellomics assay.
(a) Example plate layout for the cellular differentiation assay. Images are shown for the pluripotency markers (Oct4, Sox2, and Nanog) as they are measured in the Cellomics imaging device. Each line is measured in two rows of the same plate as technical replicates. The secondary antibody used for each marker is shown in parenthesis. Each plate also has measurements for staining with the secondary antibody only, which serves as a means to Kilpinen et  Europe PMC Funders Author Manuscripts assess background fluorescence. The red channel shows the signal from the DAPI staining, the green channel the marker signal. As expected, there is only little signal from the green channel in the wells stained only for the secondary antibody. Image acquisition stops as soon as 10,000 cells have been detected. (b) Detailed variance components of the Cellomics markers (Methods). Substantial proportions of the marker variance could be attributed to batch factors, including staining, technician effects and antibody lots. These effects mean that the fraction of cells expressing particular markers need to be interpreted with caution (Fig. 1c,d). (c) Pairwise correlation between quantitative expression scores derived from immunostaining for pluripotency and differentiation and the PluriTest score. Kilpinen   Top panel shows genomic location versus number of lines with CN three (grey) and with a CNA (black). Bottom panel shows the NAV gene score from ref22 and log 2 gene expression fold change between the iPSC lines with CN two and three (color scale), in the region highlighted in red in the top panel. Highlighted genes are up-regulated when copy number increases, known onco/tumour-suppressor genes and/or genes with NAV score in the top 2%. Extended Data Figure 5. Functional assessment of CNAs using growth assays.
Cell growth rate (a), proliferation (b) and apoptosis (c) in cell lines with copy number two ("wild type", blue dots) or copy number three ("mutant", red dots) in a recurrently duplicated region in iPSCs on chromosome 1, 17 or 20. Plot titles show the donor name and the genomic coordinates of the CNA. (a) Shown are cell counts taken on successive days in culture, for pairs of lines (one mutant, one wild type) grown on the same 24-well plates. Star symbols denote significance levels for statistical interactions between day and copy number in a linear mixed model, using fixed effects to fit day and copy number, and random effects Europe PMC Funders Author Manuscripts to account for culture plate effects. "EIF4A3" denotes whether a copy number variant overlaps one of the suspected candidate genes on chromosome 17. * -P < 0.05; ** -P < 0.01; *** -P < 0.001. (b) Protein expression level measured using Tandem Mass Tag (TMT)-based quantitation on the Q-exactive plus (labelled "QE Plus") orbitrap and a fusion (labelled "Fusion") orbitrap MS platforms. (c) Estimated fraction of fluorescing nuclei following EdU assay in mutant and wild type lines, following exposure to mitomycin ("Treated"), or in a control sample ("Untreated"). (d) Estimated fraction of fluorescing nuclei following Terminal deoxynucleotidyl transferase dUTP nick end labelling assay (TUNEL) in mutant and wild type lines, following exposure to mitomycin ("Treated"), or in a control sample ("Untreated"). Solid trend lines are least squares regression fits. P-values in b and c denote the significance of statistical interactions between copy number and mitomycin treatment condition ("Treated" or "Untreated"). Kilpinen  Association was also not deemed significant when including gender and media as fixed effects and batch variables and donor as random effects (score P = 0.3, novelty P = 0.14). (c) Passage number versus log10 RNA-seq expression of pluripotency factors Nanog and Pou5f1 (Oct4) shows no significant association between passage number and pluripotency. Trend lines are fit using linear regression of log10 expression on passage number (Nanog P = 0.5, Pou5f1 P = 0.15). Association was also not deemed significant when considering the  lead eQTLs that replicate in each tissue, with replication defined using two different replication thresholds (TH1: nominal eQTL P < 0.01/N_tissues; TH5: P < 0.05/N_tissues; plotted as dots and triangles, respectively). (d) Enrichment of alternative iPSC eQTLs ('GTEx-like") at promoter proximal and distal (defined as less than or greater than 2 Kb from the transcription start site) transcription factor binding sites (TFBS) in H1-hES cells from the ENCODE Project 49. Fold enrichments per factor are shown for iPSC-specific and non-specific eQTLs (minimum 10 observed overlaps) (Methods). Pluripotency-associated factors are indicated with an asterisk. The profile of enrichments is comparable to that obtained with the standard HipSci pipeline (Fig. 4d). Kilpinen  and identified one alternative splicing event associated with rs10069690, the lead iPSC eQTL variant for TERT (Fig. 5b). Shown is TERT intron 4 retention ratio (PSI, percent spliced in) in iPSC lines of all individual donors stratified by their genotype at rs10069690. This variant affects the splicing of the intron where it is located, with the minor allele (T) increasing the fraction of TERT transcripts in which intron 4 is retained (P = 1.7x10 -9 , Bonferroni adjusted linear regression).

Supplementary Information
Refer to Web version on PubMed Central for supplementary material.    (a) Proportion of tissue-specific eQTLs in iPSCs and 44 GTEX tissues24. (b) Most likely source of tissue-specific eQTLs in iPSCs (lead and secondary), testis and somatic tissues in GTEx (averaged; including cell lines, excluding testis). Breakdown: gene not expressed (red); gene expressed but no eQTL (blue); eQTL effect is driven by distinct lead variants (r 2 < 0.8; green). (c) Heatmap of the fold enrichment (FE) difference between iPSC-specific and non-specific eQTLs at chromatin states from the Roadmap Epigenomics Project27, shown for five aggregated clusters representing 127 cell types (SOM, somatic; PSCd, PSC-derived). Colouring: enriched for iPSC-specific eQTLs (blue), enriched for non-specific eQTLs (red).
(d) Enrichment of iPSC eQTLs at promoter proximal and distal transcription factor binding sites in H1-hES cells from the ENCODE Project49. Fold enrichments per factor are shown for iPSC-specific and non-specific eQTLs. Pluripotency-associated factors are indicated with an asterisk.  (a) Colocalised association signal for iPSC expression of PTPN2 (top) and five common diseases (bottom; inflammatory bowel disease, IBD; rheumatoid arthritis, RA; Crohn's disease, CD; celiac disease, CEL; and type 1 diabetes, T1D). PP4 is the posterior probability that the disease and gene expression associations are driven by the same causal variant29. (b) An iPSC-specific eQTL for TERT (rs10069690) that is associated with risk for breast,