1. Computational and Systems Biology
  2. Genetics and Genomics
Download icon

Systematic genetic analysis of the MHC region reveals mechanistic underpinnings of HLA type associations with disease

  1. Matteo D'Antonio
  2. Joaquin Reyna
  3. David Jakubosky
  4. Margaret KR Donovan
  5. Marc-Jan Bonder
  6. Hiroko Matsui
  7. Oliver Stegle
  8. Naoki Nariai
  9. Agnieszka D'Antonio-Chronowska
  10. Kelly A Frazer  Is a corresponding author
  1. University of California, San Diego, United States
  2. Rady Children’s Hospital, University of California, San Diego, United States
  3. European Molecular Biology Laboratory, European Bioinformatics Institute, United Kingdom
Research Article
  • Cited 0
  • Views 856
  • Annotations
Cite this article as: eLife 2019;8:e48476 doi: 10.7554/eLife.48476

Abstract

The MHC region is highly associated with autoimmune and infectious diseases. Here we conduct an in-depth interrogation of associations between genetic variation, gene expression and disease. We create a comprehensive map of regulatory variation in the MHC region using WGS from 419 individuals to call eight-digit HLA types and RNA-seq data from matched iPSCs. Building on this regulatory map, we explored GWAS signals for 4083 traits, detecting colocalization for 180 disease loci with eQTLs. We show that eQTL analyses taking HLA type haplotypes into account have substantially greater power compared with only using single variants. We examined the association between the 8.1 ancestral haplotype and delayed colonization in Cystic Fibrosis, postulating that downregulation of RNF5 expression is the likely causal mechanism. Our study provides insights into the genetic architecture of the MHC region and pinpoints disease associations that are due to differential expression of HLA genes and non-HLA genes.

Introduction

The 4 Mb major histocompatibility (MHC) region on chromosome 6p21.3, which encodes the human leucocyte antigen (HLA) gene complex, is highly polymorphic, gene dense, and has strong linkage disequilibrium (LD) (Dendrou et al., 2018; Fehrmann et al., 2011). Genome-wide association studies (GWAS) have identified many autoimmune and infectious diseases that are more prevalent in individuals carrying certain HLA types (alleles), whereas other HLA types are protective against disease (Blackwell et al., 2009; Gough and Simmonds, 2007; Matzaraki et al., 2017). The molecular mechanisms by which certain HLA types predispose or protect against disease are currently unclear. Several studies have suggested that the mechanistic underpinnings underlying associations with some autoimmune diseases may be due to certain HLA types being more likely than other HLA types to atypically present self-antigens (Klein et al., 2014; Oldstone, 1998; Yin et al., 2013). Expression differences between different HLA types of the same HLA gene could also contribute to disease associations. Finally, some of the genetic associations with specific HLA types could be due to causal variants on the same haplotype, that affect the expression or function of non-HLA genes in the MHC region (Holoshitz, 2013). Hence there are complex interactions between genetic variation, gene expression and function of both HLA and non-HLA genes, and disease, rendering the MHC region an important but difficult interval to genetically interrogate.

Here, we used deep whole genome sequence (WGS) from 419 individuals to call full resolution (eight-digit) HLA types and RNA-seq data from 361 matched induced pluripotent stem cells (iPSCs) to create a comprehensive map of regulatory variants in the MHC region. We identified eight-digit HLA types that are strongly associated with expression differences of their cognate HLA gene, and conducted eQTL analyses using both single-variants and eight-digit HLA types, which represent haplotypes of coding and regulatory variants. We built upon the comprehensive MHC regulatory map by functionally exploring GWAS signals for 4083 complex traits, finding strong evidence of colocalization for 180 disease loci with eQTLs. Our data are a valuable resource for both the genetics and disease-focus communities, providing insights into the genetic architecture of the MHC region that have not previously been described. To demonstrate the utility of our resource, we examined the common Caucasian 8.1 ancestral haplotype (8.1AH), which spans the entire MHC region and is associated with many autoimmune disorders and other diseases (Gambino et al., 2018), but is protective against bacterial colonization in Cystic Fibrosis (CF) patients (Tomati et al., 2015). Here, we show that 8.1AH is associated with decreased expression of RNF5, and previous studies have shown that inhibition of RNF5 results in rescue of F508del mutant cystic fibrosis transmembrane conductance regulator (CFTR) function (Sondo et al., 2018). Based on these observations, we postulate that downregulation of RNF5 expression is the causal mechanism underlying the association between 8.1AH and delayed colonization in CF patients.

Results

Eight-digit HLA typing

We used HLA-VBSeq (Nariai et al., 2015) to estimate HLA types at eight-digit resolution for 30 genes [six HLA class I genes (HLA-A, -B, C, -E, -F, -G); eight HLA Class I pseudogenes (HLA-H, -J, -K, -L, -P, -T, -V, -W); twelve HLA class II genes (HLA-DMA, -DMB, -DOA, -DOB, -DPA1, -DPA2, -DPB1, -DPB2, -DQA1, -DQB1, -DRA, -DRB1); and four non-HLA genes (MICA, MICB, TAP1, and TAP2)] in 650 WGS samples (from 419 unique individuals) and identified 526 unique alleles (Figure 1, Supplementary file 1, Supplementary file 2). The 650 WGS samples were obtained by combining data from two large induced pluripotent stem cell (iPSC) resources: 1) iPSCORE (273 individuals) (Panopoulos et al., 2017), and 2) HipSci (377 samples from 146 individuals) (Kilpinen et al., 2017a; Kilpinen et al., 2017b; Streeter et al., 2017). All the WGS samples had deep coverage and the expected high polymorphism rate across the MHC region (Figure 2A–B). The 146 individuals in HipSci were of European descent while the 273 iPSCORE individuals are of diverse ethnic backgrounds (190 European, 30 Asian, 20 Multiple ethnicities, 18 Hispanic, 7 African American, 5 Indian and 3 Middle Eastern) and thus harbored HLA alleles from multiple human populations (Supplementary file 3). While all 146 HipSci donors were unrelated, there were 183 donors in iPSCORE that are part of 56 unique families (2–14 individuals/family), including 25 monozygotic (MZ) twin pairs. For 84 HipSci individuals, we HLA typed 231 matched fibroblast-iPSC WGS pairs (one fibroblast sample and two or more iPSC clones per donor).

HLA nomenclature.

(A) Description of an HLA type (allele) using current standard nomenclature. The first two digits typically correspond to the serological antigen. The third and fourth digits distinguish HLA types that differ by one or more missense variants. The fifth and sixth digits distinguish HLA types that differ by one or more synonymous variants. The seven and eight digits distinguish HLA types that differ by noncoding variants. (B) The location of the MHC region on chromosome 6 (left). The 8.1 ancestral haplotype (8.1AH) is comprised of six HLA types (corresponding to six HLA genes) typically reported at four-digit resolution (middle). For each of the six HLA genes, we identified one to three eight-digit resolution HLA types that can be collapsed into their 8.1AH four-digit counterparts (right).

Figure 2 with 1 supplement see all
WGS coverage in MHC region and HLA types.

(A) SNP density distribution in consecutive 10 kb bins in the HLA region (chr6: 29690551, 33102442). The dashed lines represent the chromosome-wide SNP density in the iPSCORE (blue) and HipSci (red) WGS samples. The 30 HLA genes analysed in this study are represented by diamonds with green representing HLA Class I, purple HLA Class I pseudogenes, orange HLA Class II and yellow non-HLA genes. The four-digit HLA types for the six genes composing the 8.1AH are shown in red. (B) Read depth coverage of the HLA region for the iPSCORE WGS samples (top, red) and the HipSci WGS samples (bottom, blue) was calculated using consecutive 10 kb bins. Each line represents a single WGS sample, the black dashed lines represent the average coverage for the iPSCORE WGS samples (52X) and the HipSci WGS samples (38X). (C, D) Fraction of the 273 iPSCORE (C) and 377 HipSci (D) WGS samples with homozygous, heterozygous and undetermined alleles for each of the 30 HLA genes. For each sample, if one of the two alleles for an HLA gene was undetermined, we considered the HLA genotype as ‘undetermined’. (E) Twin concordance (iPSCORE: 25 monozygotic twin pairs), (F) Mendelian inheritance concordance (iPSCORE: 17 trios) and (G) Fibroblast-iPSC concordance (HipSci: 231 fibroblast-iPSC pairs) of HLA types at eight-digit resolution. Genes were sorted based on their genomic position and colored according to their class, as shown in panel A. The red dashed line represents the median across all genes.

HLA types have high recall rates and reproducibility

For each HLA gene, we initially calculated recall rate, measured as the fraction of individuals that could be HLA-typed. We observed a high recall rate for all 30 genes (mean = 98.5% for both iPSCORE and HipSci samples, Figure 2C,D) and for only three HLA genes in iPSCORE (HLA-H, HLA-T and HLA-K) and four genes in HipSci (HLA-H, HLA-T, HLA-K and HLA-DRB1) fewer than 95% of individuals were HLA-typed. The fraction of heterozygous and homozygous HLA types across the 30 genes was highly similar in the iPSCORE and HipSci individuals (r = 0.92). These results show that we were able to obtain high recall rates for all 30 MHC region genes.

To examine the reproducibility of the eight-digit HLA types, we determined HLA type concordance across the 25 MZ twin pairs (defined as twin concordance), Mendelian inheritance concordance from 17 non-overlapping trios and HLA type concordance of the 231 fibroblast-iPSC pairs (defined as HipSci concordance). Overall, we found that the median concordance across for all 526 HLA alleles was very high; however, both the twin pairs (96.7%) and fibroblast-iPSC pairs (95.7%) had slightly higher concordance than the Mendelian inheritance (90.9%, Figure 2E,F,G). These results show that HLA-VBSeq produces highly reproducible results when samples with exactly the same two HLA alleles are analyzed; but has slightly lower reproducibility calling the same HLA alleles when in different diplotypes.

To estimate the accuracy of the 526 called HLA alleles we tested for Hardy-Weinberg equilibrium (HWE) in the unrelated individuals and found that the vast majority (511, 97.1%) were in HWE (p>1×10−6, Figure 2—figure supplement 1A–C), consistent with the reproducibility values in twin pairs (96.7%) and in fibroblast-iPSC pairs (95.7%). We also conducted a comparative analysis with HLA alleles in 3.5 million individuals with diverse genetic backgrounds from the Allele Frequency Net Database (AFND) (González-Galarza et al., 2015) and observed high correlation between HLA allele frequencies in AFND and iPSCORE (r = 0.921) and HipSci (r = 0.908, Figure 2—figure supplement 1D). Overall, HLA-VBSeq (Nariai et al., 2015) was able to accurately detect eight-digit resolution HLA alleles in the 650 WGS samples as evidenced by high recall rates, reproducibility and allele frequencies consistent with Hardy-Weinberg law.

Eight-digit HLA types associated with expression of cognate HLA gene

To analyze gene expression, we used RNA-seq data from 446 iPSC samples (derived from 361 of the HLA typed individuals Supplementary file 2). There were 146 expressed genes in the MHC interval (TPM >2 in at least 10 samples), of which 24 were HLA genes, including six HLA class I genes, 3 HLA class I pseudogenes, 11 HLA class II genes, MICA, MICB, TAP1 and TAP2. For each individual, we estimated expression levels of the 24 HLA genes by aligning reads to cDNA sequences specific for the HLA types called by HLA-VBSeq (Nariai et al., 2015) to avoid alignment biases (Aguiar et al., 2019; Gensterblum-Miller et al., 2018; Lee et al., 2018; Panousis et al., 2014) (Figure 3—figure supplement 1, Figure 3—figure supplement 2, Supplementary file 4). The HLA class I genes tended to be expressed at high levels, consistent with their ubiquitous expression in all cell types, and the HLA class II genes tended to be expressed at lower levels, as expected due to their primary role in immune cells (Matzaraki et al., 2017). To examine the extent that HLA types at eight-digit resolution were associated with gene expression levels, we compared the expression levels of each HLA type against all other alleles from the cognate HLA gene. We quantile-normalized expression of each gene and compared the TPM distributions between all samples carrying a specific HLA allele and all the other samples (Figure 3A,B). We observed high variability in the expression levels of the different HLA alleles, with 189 of 428 HLA types (44.2%) present in at least two individuals showing a significantly different expression level of the cognate HLA gene compared with the samples not carrying the HLA type (t-test, FDR < 5%, Supplementary file 5). Five HLA genes (HLA-A, HLA-C, HLA-DRB1, HLA-DQB1 and HLA-DPB1) showed more than four-fold differences between the least expressed and the most expressed alleles. We examined the HLA types corresponding to the 8.1AH at eight-digit resolution (Figure 1B), and determined that the class I genes tended to be expressed higher than the other HLA types of the cognate gene, while two of the class II genes (HLA-DQA1 and HLA-DQB1) were significantly expressed lower and one (HLA-DRB1) was not differentially expressed (Figure 3C–H). Overall, eight-digit HLA types were strongly associated with expression differences of the cognate HLA gene, suggesting that the associations between HLA types and disease may in part be driven by differential allelic expression.

Figure 3 with 2 supplements see all
Allele-specific expression of HLA genes in iPSCs.

(A, B) Mean TPM expression across HLA types: (A) shows the four HLA genes with highest expression levels (mean TPM >20); and (B) shows lower expressed HLA genes (mean TPM <20). Each point represents one allele. Point size represents the number of individuals carrying each allele. All HLA types are shown, although differential gene expression analysis (Supplementary file 5) was performed only on those present in at least two individuals. Eight-digit HLA-types significantly associated with the expression of their cognate HLA gene are shown in dark blue. (C–H) Boxplots showing mean expression of the eight-digit 8.1AH HLA types (Figure 1) compared with the mean expression of all the other HLA types for the cognate gene. Red stars indicate that the mean expression of the ancestral haplotype HLA type is significantly different (Supplementary file 5).

eQTLs associated with the expression of 56 eGenes in the MHC region

We performed an eQTL analysis by testing each of the 146 expressed genes and the genotypes of 45,245 variants (MAF >1%) spanning the MHC region (±1 Mb flanking each side) using a linear mixed model. We detected associations between 56 eGenes and 72,473 eQTLs (Benjamini-Hochberg-corrected p-value<0.1, Supplementary file 6) corresponding to 26,363 distinct variants; hence 36.4% of the variants in the MHC region were associated with the expression of one or more gene. Due to the low recombination rate (DeGiorgio et al., 2014; Trowsdale and Knight, 2013) and complex LD structure (Jensen et al., 2017; Miretti et al., 2005) in the MHC region, we investigated the distributions of all primary eQTLs rather than solely focusing on lead variants (Figure 4). The 56 eGenes tended to be associated with high numbers of primary eQTLs (median = 842), ranging from 26 (HLA-E) to 4911 (HLA-DQB1). Many of the eGenes were also associated with primary eQTLs spanning relatively large distances (median = 2.1 Mb), including ZFP57, HLA-C, HLA-B, MICA and RNF5, whose associated primary eQTLs nearly spanned the entire 4 Mb MHC region, which is consistent with the strong LD of ancestral haplotypes. While genome-wide eQTL analyses usually test associations between gene expression and neighboring variants (within 500 kb or 1 Mb), our results show that gene expression in the MHC is associated with cis regulatory variants that span regions several times larger than those normally tested.

Figure 4 with 5 supplements see all
Associations between variants and gene expression levels in the MHC region.

eQTL associations for 56 genes in the MHC region. Colors represent whether each eQTL is primary or conditional: ‘conditional 1’ refers to eQTLs independent from the lead eQTL and ‘conditional 2’ refers to eQTLs independent from the lead eQTL and the lead ‘conditional 1’ eQTL. TSS for each gene is shown in red.

To identify independent eQTLs, we repeated the analysis conditioning gene expression on the lead eQTL for each of the 56 eGenes. For 42 eGenes, we observed secondary eQTLs (conditional on the lead eQTL), and for 36 eGenes we identified a tertiary eQTL upon conditioning on the top two independent eQTLs (Figure 4, Figure 5). Overall, for the 56 eGenes, we observed 76,809 eQTLs (72,473 primary, 2424 conditional 1 and 1912 conditional 2, Bonferroni-adjusted p-value<0.1) with genome-wide significance corresponding to 27,200 distinct variants. These results show that, within the gene-dense MHC region, a large fraction (56; 38.4%) of the expressed genes are eGenes, many of which (42; 75%) are associated with multiple independent regulatory variants that span long distances.

Associations between HLA types and gene expression levels in the MHC region.

(A) Heatmap showing the associations between HLA types at 8-digit resolution (Y axis) and the expression of 114 genes in the MHC region (X axis). For eGenes associated with multiple HLA-type eQTLs, only the most significant p-value is shown. Self-associations between the expression of each HLA gene and its HLA types were not considered (grey squares). Four different groups of genes with shared HLA-eQTLs are outlined in different colors. (B) Barplot showing for each of the 114 genes in the heatmap in (A), the number of primary and conditional associations (adjusted based on the lead eQTL, the top two independent eQTLs, or the top three independent eQTLs) between HLA types and gene expression levels. ‘Conditional 1’ refers to HLA-eQTLs conditional on the lead primary SNP eQTL, ‘Conditional 2’ refers to HLA-eQTLs conditional on the lead primary and lead conditional 1 eQTLs and ‘Conditional 3’ refers to HLA-eQTLs conditional on the lead primary, conditional 1 and conditional 2 eQTLs. (C) Scatterplot showing the number of primary eQTLs (X axis) and the number of primary HLA-type eQTLs (Y axis) for the 52 genes that had both eQTLs and HLA-type eQTLs.

Given the complexity of conducting genetic association studies in the MHC region, we set out to validate the eQTLs that we discovered by showing that they have characteristics of regulatory variants. We initially investigated the distance between lead eQTLs and the transcription start site (TSS) of their associated eGene and found that in general they were close (median distance = 18.3 kb, Figure 4—figure supplement 1), and for 13 eGenes (23.2%), the lead eQTL was localized within 5 kb of their TSS. We next visually examined the association between the expression levels of each of the 56 eGenes and genotypes of the lead eQTL and observed a strong correlation (Figure 4—figure supplement 2). We examined if the 76,809 primary and conditional eQTLs (27,200 distinct variants) were enriched in regulatory elements by testing their overlap with chromatin states in iPSC lines from the Roadmap Epigenomics Consortium (Kundaje et al., 2015). We observed that the eQTLs in the MHC region with the largest effect sizes were enriched for overlapping regulatory elements active in iPSCs including active transcription start sites (TSS), regions flanking active TSS, UTRs (defined as regions transcribed at 5’ and 3’) and transcribed regions (Figure 4—figure supplement 3). To further characterize the associations between eQTLs and gene expression, we tested if the 56 eGenes were more likely to have allele-specific expression (ASE) than the other 90 expressed MHC genes; and found that eGenes were significantly more likely to have ASE [mean allelic imbalance fraction (AIF) = 67.3%, compared with 61.9% in the other genes, p=4.38×10−64, Mann-Whitney U test, Figure 4—figure supplement 4, Supplementary file 7]. We also observed that eGenes in samples with heterozygous eQTLs (primary and conditional) have stronger ASE compared with eGenes in samples homozygous for eQTLs (39,891 eQTLs in total examined, corresponding to 21,739 variants; Figure 4—figure supplement 5). These results show that the primary and conditional eQTLs in the HLA region have characteristics of regulatory variants including proximity to the TSS of their associated eGene, enrichment in active chromatin marks, and association with allele-specific expression of their eGene.

Regulatory variants in the MHC region play important roles in complex traits

To characterize the associations between genetic variation in the MHC region, gene expression and disease, we investigated if eGenes and complex traits shared the same causal variants using a colocalization approach (Giambartolomei et al., 2014). We tested the eQTL signals (primary and conditional) for all 56 MHC eGenes against the GWAS signals for 4083 complex traits from the UK Biobank and as expected identified many strong associations (1480 loci with posterior probability of association (PPA) >0.2 for sharing the same causal variant, including 180 with PPA >0.8, Supplementary file 8).

To demonstrate the utility of using our resource, we explored colocalization between eGenes and complex traits that contribute to wide phenotypic variation in Cystic Fibrosis (CF) patients including pulmonary diseases (Stoltz et al., 2015), weight (Leung et al., 2017; Morison et al., 1997) as well as liver, biliary tract and gallbladder diseases (Assis and Debray, 2017; Diwakar et al., 2001; Herrmann et al., 2010) (Figure 6). We identified eight eGenes whose eQTL colocalized with GWAS signals for pulmonary diseases, of which three, HLA-B, HLA-DQA1 and HLA-DQB1, had been previously associated with Asthma in GWAS (Kontakioti et al., 2014; Li et al., 2012; Mahdi et al., 2018) and in the UK Biobank (Vicente et al., 2017). However, the other five eGenes, including HLA-DRA (the most strongly associated), APOM, DXO, RNF5 and TAP1 were novel, that is had not been identified in GWAS. Interestingly, the association of RNF5 was consistent with its function as a regulator of the CFTR protein and suggested role in cystic fibrosis (Sondo et al., 2018; Tomati et al., 2015). We next identified four novel eGenes (EHMT2, GPANK1, HLA-C and RNF5) associated with weight phenotypes. The two eGenes most strongly associated, GPANK1 and RNF5, had opposite associations. GPANK1 was only associated with fat-free mass, suggesting that the protein product is involved in muscle or bone mass regulation. Conversely, RNF5 was more strongly associated with fat mass than fat-free mass, suggesting that its altered expression can contribute to CF phenotypic variation, as these individuals are known to have bile and other digestive issues, including difficulties in assimilating fat (Bijvelds et al., 2005; Freudenberg et al., 2008; Morison et al., 1997; Wilke et al., 2011). Finally, eQTLs for HLA-DQB1, PSORS1C2 and RNF5 colocalized with liver and gallbladder GWAS signals. These findings suggest that altered gene expression levels in the MHC region are associated with phenotypic variation in CF patients. They also demonstrate the utility of our data to serve as a resource to narrow down the location of causal regulatory variants underlying known associations between genes in the MHC region and complex traits, as well as to identify biologically meaningful novel associations.

Colocalization results.

Heatmap showing the posterior probability of association (PPA) for colocalization (Giambartolomei et al., 2014) between an eQTL for the indicated eGene (columns) and a GWAS signal for the indicated trait (rows). GWAS traits in Supplementary file 8 were filtered to include those having at least one eGene with PPA > 0.5 and the following terms in their description: ‘asthma’, ‘pulmonary’, ‘mass’, ‘weight’, ‘fat’, ‘liver diseases’ and/or ‘biliary tract diseases’. Associations with PPA > 0.2 for all 4,083 UKBB traits are reported in Supplementary file 8.

HLA-type eQTLs associated with the expression of 114 eGenes in the MHC region

Given that eight-digit HLA types represent a haplotype of many regulatory and coding SNPs (Figure 1), we examined if conducting an eQTL analysis using HLA types would detect regulatory signals independent from those detected in the single-variant eQTL analysis. We calculated associations between each of the 283 common eight-digit HLA-types (MAF ≥1% in 361 individuals with expression data) and the 146 genes expressed in the HLA region, a linear mixed model (LMM) (see Materials and methods). We first examined self-associations between the expression of each HLA gene and its HLA types (for example, we tested the associations between HLA-A expression and HLA-A:01:01:01:01). We found that for all 24 expressed HLA genes, there were at least two alleles significantly associated (referred to as HLA-type eQTLs hereafter, Benjamini-Hochberg-corrected p-value<0.05) with its expression (median = 4, range 2 to 11; Supplementary file 9), which is consistent to what we observed in the allele-specific expression analysis (Figure 3). Self-associations were not considered for the remaining HLA-type eQTL analyses. We detected 1139 associations between 250 HLA-type eQTLs and 114 eGenes (Supplementary file 9, Figure 5A,B). Associations primarily occurred between HLA-type eQTLs and the expression of genes located in their proximity (median distance = 146.2 kb for significant associations compared with 1.35 Mb for non-significant associations, p=2.7×10−67, Mann-Whitney U test), resulting in HLA-type eQTLs largely clustering into four distinct genomic regions. Interestingly, RNF5 showed strong associations with multiple HLA types in clusters 1, 2 and 3, confirming the observation that it had long-range associations with eQTLs spanning the entire MHC region (Figure 4). In total, there were 62 eGenes identified only by the HLA-type eQTL analysis, four identified only by the single-variant eQTL analysis and 52 identified by both analyses (44.4% all 118 eGenes; Supplementary file 10). These results show that the HLA-type eQTL analysis had substantially greater power to identify eGenes than the single-variant eQTL analysis.

We sought to determine if HLA-type eQTLs provided independent signals from the single-variant eQTL analysis in the 52 overlapping eGenes. We observed that the number of HLA-type eQTLs was highly correlated with the number of eQTLs (r = 0.829, p=7.9×10−31, Figure 5C). To identify HLA-type eQTL independent signals, we conditioned gene expression on the lead eQTLs from the single-variant analysis. We observed 406 HLA-type eQTLs conditional on the lead eQTLs (49 eGenes), 316 HLA-type eQTLs conditional on the top two independent eQTLs (27 eGenes) and 150 HLA-type eQTLs (13 eGenes) conditional on the top three independent eQTLs (Figure 5B, Supplementary file 9). Of note, the conditional HLA-eQTLs were largely confined to the four genomic regions containing clusters of genes with shared HLA-eQTLs and the RNF5 gene. These findings show that the HLA-eQTL analysis identifies a different set of regulatory variants than the single-variant eQTL analysis even for the 52 eGenes discovered by both analyses.

HLA-type eQTLs capture regulatory variation associated with expression of neighboring genes

We investigated the molecular underpinnings of the clustering of the most significant HLA-type eQTLs within the four distinct genomic regions (Figure 5A): 1) chr6:29640168–30152231, including HLA class I genes (HLA-A, HLA-F and HLA-G) and HLA class I pseudogenes; 2) chr6:31171745–31512238, including HLA class I genes (HLA-B, HLA-C), MICA and MICB; 3) chr6: 32361463–32789609, including HLA class II genes (HLA-DRA, HLA-DRB1, HLA-DQA1, HLA-DQB1 and HLA-DOB), TAP1 and TAP2; and 4) chr6:32916389–33115544, including HLA class II genes (HLA-DPA1, HLA-DPB1, and HLA-DPB2). We observed that different eGenes in each cluster were typically associated with the genotypes of different HLA types, suggesting that the observed clustering was not due to co-expression (i.e., correlated expression levels) of neighboring HLA genes (Figure 7). There was one exception in cluster 3, with the expression of both HLA-DRB5 and HLA-DRB1 associated with many of the same HLA types (Figure 7C). For this reason, we examined the co-expression of the eGenes in each of the four clusters. Only five eGenes in cluster 3 (HLA-DRA, HLA-DRB5, HLA-DQA1, HLA-DQB1 and HLA-DRB1) displayed moderate correlated expression (Figure 7—figure supplement 1), despite having their expression levels largely associated with different HLA types of neighboring genes (Figure 7). These findings suggest that the HLA-type eQTLs capture regulatory variants associated with the expression of one or a few neighboring genes. Furthermore, unlike the eQTLs from the single-variant analysis that tend to span long distances most likely due to the strong LD in the MHC region (Figure 4), the HLA-type eQTLs capture regulatory variants that act over relatively short distances.

Figure 7 with 1 supplement see all
HLA-type eQTLs for HLA genes in the four clusters shown in Figure 6.

Heatmaps showing, for each of the four clusters identified in Figure 6, the associations between HLA-type eQTLs and gene expression. Each column represents an eGene, while rows represent all HLA-type eQTLs associated with other genes in the cluster. 8.1AH HLA types are highlighted in red. The heatmaps show that the HLA-type eQTLs capture regulatory variants associated with the expression levels of neighboring genes.

RNF5 expression levels underlie association between 8.1AH and colonization in CF

While we were able to colocalize the single-variant eQTLs and GWAS signals (Figure 6), this was not possible for HLA-type eQTLs, because GWAS for complex traits have not yet been conducted using eight-digit HLA types. Therefore, we decided to perform focused analyses of associations between HLA-type eQTLs, RNF5, and the 8.1AH to gain insight into the molecular mechanisms underlying CF phenotypic variation. The expression of RNF5 was negatively associated with eight-digit HLA types for all six genes composing the 8.1AH (Figure 8A–F, Figure 8—figure supplement 1, Supplementary file 6). While we discovered several eight-digit resolution HLA types for some of the HLA genes composing 8.1AH (Figure 1), only one eight-digit HLA type for each of the six genes was associated with RNF5 expression. It has been previously proposed that the association between 8.1AH and delayed bacterial colonization in CF patients could be due to the impact of the ancestral HLA-types on the microbiota composition in the lungs (Laki et al., 2006). However, our findings, combined with the fact that RNF5 is a current drug target in CF patients because its downregulation contributes to stabilizing and rescuing the function of mutant CFTR proteins (Sondo et al., 2018; Tomati et al., 2015), suggests that the association between 8.1AH CF carriers and delayed bacterial colonization could be due to decreased expression of RNF5 (Lyczak et al., 2002; Mall and Hartl, 2014; McNicholas, 2017; Pier et al., 1996). In this proposed model, CF patients carrying 8.1AH express RNF5 at lower levels than non-carriers, resulting in lower protein levels of RNF5 and thereby less degradation of the misfolded mutated CFTR protein, improved Cl- secretion, lower mucus secretion and delayed colonization by S. aureus and P. aeruginosa (Figure 8F).

Figure 8 with 1 supplement see all
Associations between HLA-type eQTLs, RNF5 expression and Cystic Fibrosis.

(A–F) Distributions of RNF5 expression levels in samples that are homozygous (red) and heterozygous (yellow) for significant HLA-type eQTLs corresponding to 8.1AH HLA alleles. RNF5 expression distributions for all samples that do not carry the indicated 8.1AH HLA alleles are shown in grey. RNF5 expression distributions for all significant HLA-type eQTLs that are not part of the 8.1AH are shown in Figure 8—figure supplement 1. (G) Cartoon depicting proposed model of the molecular mechanisms underlying the associations between 8.1AH, RNF5 expression and CFTR function in CF. On the left, in healthy individuals, CFTR is correctly folded in the airway epithelium and Cl- ions can be secreted. On the right, in CF patients not carrying 8.1AH, CFTR is misfolded and high levels of RNF5 result in its degradation which, in turn, results in decreased Cl- secretion, mucus hypersecretion and colonization by S. aureus and P. aeruginosa. In the middle, in CF patients carrying 8.1AH, RNF5 is expressed at low levels, resulting in the lower levels of RNF5, less degradation of the misfolded mutated CFTR protein, improved Cl- secretion, lower mucus secretion and delayed colonization by S. aureus and P. aeruginosa.

Discussion

The vast majority of risk loci identified through GWAS are located in non-coding regions and the causal variants are implicitly assumed to affect gene regulation. Colocalization methods integrating GWAS and eQTL studies can identify causal variants that are shared between a complex trait locus and altered gene expression; and thereby elucidate the mechanism underlying GWAS non-coding variants and identify their target genes. Although the human MHC region on 6p21.3 has been associated with many autoimmune and infectious diseases, it is typically excluded from eQTL studies due to its high SNP frequency and complex LD structure (Lam et al., 2017). To address this gap, we developed a comprehensive map of regulatory variation in the MHC region using deep WGS from 419 individuals and RNA-seq data from matched iPSCs. We have previously shown that iPSCs are well-powered for eQTL mapping (Bonder et al., 2019; Jakubosky et al., 2019a; Jakubosky et al., 2019b), have a distinct regulatory landscape relative to somatic tissues (DeBoever et al., 2017), and thereby provide insights into regulatory variants that exert their effects during early development. While we demonstrate the feasibility of generating a map of regulatory variation for the MHC region and its utility for examining molecular mechanisms underlying disease associations in the interval, future studies using eQTLs from other cell types could increase the number of complex trait loci in the interval that are functionally annotated.

We constructed a large panel of HLA types at eight-digit resolution for 30 genes in the MHC region using a state-of-the-art computational algorithm, HLA-VBSeq (Nariai et al., 2015), and the comprehensive IPD-IMGT/HLA database (Robinson et al., 2016). To examine the extent that eight-digit HLA types were associated with expression levels of their cognate HLA gene, we estimated expression levels of the HLA genes by aligning RNA-seq reads to cDNA sequences specific for the HLA types called by HLA-VBSeq. Overall, we observed that the eight-digit HLA types were strongly associated with expression differences of the cognate HLA gene. This finding suggests that the associations between certain HLA types and disease may not only be driven by polymorphic amino acid epitopes resulting in differential binding affinities to specific peptides, but also in part by different expression levels.

We identified 56 eGenes in the MHC region that were associated with high numbers of primary eQTLs (median = 842) spanning relatively large distances (median = 2.1 Mb), which is consistent with the strong LD in the interval. We also conducted conditional eQTL analyses and discovered that the majority of these eGenes were associated with multiple independent regulatory variants. Given the complexity of conducting genetic association studies in the MHC region, we examined the primary and conditional eQTLs and showed that they have characteristics of regulatory variants including proximity to the TSS of their associated eGene, enrichment in active chromatin marks, and association with allele-specific expression of their eGene. We then tested the eQTL signals (primary and conditional) for all 56 eGenes against GWAS signals for 4083 complex traits from the UK Biobank and found strong evidence of colocalization with 180 disease loci. Given that the 8.1 ancestral haplotype (8.1AH) is an important genetic modifier of lung disease in Cystic Fibrosis (CF) patients (Laki et al., 2006), we examined if any of the MHC eGenes colocalize with loci for complex traits contributing to phenotypic variation in CF, including pulmonary diseases (Stoltz et al., 2015), weight (Leung et al., 2017; Morison et al., 1997), as well as liver, biliary tract and gallbladder diseases (Assis and Debray, 2017; Diwakar et al., 2001; Herrmann et al., 2010). Our findings suggest that altered gene expression levels in the MHC region contributes to phenotypic variation in these complex traits and demonstrates the utility of our data to identify shared causal variants underlying GWAS and eQTL signals in the interval.

We hypothesized that HLA types, which represent a haplotype of many regulatory and coding SNPs, could exert a greater effect on gene expression than single variants. We conducted an HLA-type eQTL analysis and identified 114 eGenes, showing that indeed taking HLA type haplotypes into account substantially increased power to identify eGenes over single-variant eQTL analysis. HLA-type eQTLs largely clustered into four distinct genomic regions. We showed that this clustering was due to the fact that HLA-type eQTLs capture regulatory variants that act over relatively short distances and tend to be associated with the expression of only one or a few neighboring genes. Unfortunately, we were not able to conduct a colocalization analysis of HLA-type eQTLs and GWAS signals, because GWAS for complex traits have not yet been conducted using eight-digit HLA types. However, in the past few years, novel methods have been developed to impute HLA types from variants assayed by SNP arrays (Jia et al., 2013; Zheng et al., 2014), which makes conducting GWAS using HLA types feasible for currently existing datasets. Based on our findings, future GWAS directly testing associations between HLA types and complex traits are likely to uncover many novel associations.

We examined the association between the 8.1AH and delayed colonization in cystic fibrosis (CF) (Laki et al., 2006). We showed that RNF5 expression was negatively associated with eight-digit HLA types for all six genes composing the 8.1AH. Interestingly, RNF5 was recently described as a potential novel drug target in CF, because, when RNF5 activity is decreased, the mutant CFTR protein becomes stabilized and its function as a cAMP-dependent chloride channel in the lung epithelium is partially rescued (Sondo et al., 2018; Tomati et al., 2015). Impaired ion transport mediated by CFTR is associated with reduced airway hydration and decreased mucociliary clearance, which leads to increased vulnerability to bacterial infection (Dechecchi et al., 2018; Munder and Tümmler, 2015). Hence, we postulate that downregulation of RNF5 expression due to the presence of specific regulatory variants in 8.1AH is the likely mechanism underlying delayed colonization by S. aureus and P. aeruginosa. While our findings do not preclude that HLA types comprising 8.1AH have a direct role in the protection against bacterial colonization in CF, they suggest that reduced RNF5 expression likely plays a major factor in the association between 8.1AH and colonization.

Taken together, our data provide insights into the genetic architecture of the MHC region and show that disease associations with HLA types cannot only be due to epitope differences but differential expression of HLA genes and associated non-HLA genes. Further, these data are a valuable resource for both the genetics and disease-focus communities that can be exploited to identify causal mechanisms underlying disease associations in the MHC region, and potentially novel drug targets.

Materials and methods

Key resources table
Reagent type
(species) or resource
DesignationSource or referenceIdentifiersAdditional
information
Software, algorithmGATK(Van der Auwera et al., 2013; Zhang et al., 2017)
Software, algorithmHLA-VBSeq(Nariai et al., 2015)
Software, algorithmSAMtools version 1.2(Li et al., 2009)
Software, algorithmSTAR (2.5.0a)(Dobin et al., 2013; Harrow et al., 2012)
Software, algorithmRSEM (1.2.20)(Li and Dewey, 2011)
Software, algorithmbcftools(Li et al., 2009)
Software, algorithmlmekin (coxme package)R 3.5.1
Software, algorithmMBASED(Mayba et al., 2014)
Software, algorithmcoloc.abf (coloc package)R 3.5.1 (Giambartolomei et al., 2014)

Whole-genome sequencing data

Request a detailed protocol

iPSCORE Resource: Whole genome sequencing (WGS) was performed for 273 individuals (152 females and 121 males) at Human Longevity, Inc (HLI) as described in DeBoever et al. (2017). Briefly, DNA was isolated from blood or skin fibroblasts (254 and 19 samples, respectively) using DNEasy Blood and Tissue Kit. A Covaris LE220 instrument was used to quantify, normalize and shear the DNA samples, which were then normalized to 1 µg, DNA libraries were made using the Illumina TruSeq Nano DNA HT kit and their size and concentration determined using LabChip DX Touch (Perkin Elmer) and Quant-iT (Life Technologies), respectively. We normalized all concentrations to 2–3.5 nM, prepared combined pools of six samples and sequenced them at 150 bp paired-end on Illumina HiSeqX. WGS data were generated at an average of 1.1 billion paired-end reads (depth of coverage: 52X) per sample (range 720 million to 3.1 billion). WGS reads were aligned to the reference genome (GRCh37/hg19) including decoy sequences (hs37d5) (Auton et al., 2015) using BWA-MEM version 0.7.12 (Li and Durbin, 2010). Fastq file quality was estimated using FASTQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). All BAM files were processed using GATK best practices (Van der Auwera et al., 2013; Zhang et al., 2017), as described in DeBoever et al. (2017) to detect single nucleotide polymorphisms (SNPs) genome-wide. For chromosome 6, the chromosome-wide SNP density (79 SNPs per 10 kb) was calculated by dividing the total number of SNPs in the 273 iPSCORE individuals by the length of the chromosome, excluding undefined ‘N’ nucleotides.

HipSci Resource: WGS was performed for 146 individuals (85 females and 61 males) and 231 associated iPSC lines. Fibroblasts obtained from each HipSci individual were reprogrammed and between one to two iPSC lines were derived. We downloaded CRAM files for the 377 samples (146 fibroblast samples and 231 iPSCs) from the European Nucleotide Archive (ENA, https://www.ebi.ac.uk/ena, study PRJEB15299) and transformed to BAM files using Samtools (Li et al., 2009). The WGS data were an average of 835 million paired-end reads (depth of coverage: 38X) per sample (range 450 million to 1.8 billion). BAM files were processed using GATK to detect SNPs on chromosome 6.

WGS coverage and distribution of variants in MHC region

Request a detailed protocol

We assessed the alignment and variant calling of the WGS in the MHC region (chr6:29640168–33115544) by determining SNP density in consecutive 10 kb windows across the region. The SNP density in the WGS data varied widely: 1) for the 273 iPSCORE genomes from 1 SNP/10 kb to 900 SNPs/10 kb (mean = 149.2 SNPs/10 kb) compared with the average of 79.6 SNPs/10 kb across chromosome 6 (Figure 1A); and 2) for the 377 HipSci genomes (146 fibroblast samples and 231 iPSCs) from 1 SNP/10 kb to 995 SNPs/10 kb (mean = 155.4 SNPs/10 kb) compared with the average of 58.1 SNPs/10 kb across chromosome 6. The region with the highest SNP density contained the HLA-DQA1 and HLA-DQB1 genes, which is consistent with previous studies (Norman et al., 2017). To determine whether the high SNP density interfered with read alignment, we analyzed the read depth across the region. We observed that the read depth (mean 49.8 ± 12.1 for iPSCORE and 37.9 ± 7.1 for HipSci), overall, was comparable with the genome average (52X and 38X, respectively, Figure 1B). We found, on average, 98.1% and 96.5% of the MHC region had a read depth >20X, and for each iPSCORE and HipSci sample only 284 ± 122 bp and 252 ± 115 bp respectively were not covered by any reads.

HLA typing from whole-genome sequencing data with HLA-VBSeq

Request a detailed protocol

HLA-VBSeq (Nariai et al., 2015) estimates the most probable HLA types from WGS data at eight-digit resolution by simultaneously optimizing read alignments to a database of HLA type sequences and the abundance of reads on the HLA types by variational Bayesian inference. HLA typing was carried out as follows. First, we identified 30 target genes that had more than one allele at eight-digit resolution in the IMGT/HLA database (release 3.30.0) (Robinson et al., 2016), which included six HLA class I genes (HLA-A, -B, C, -E, -F, -G), eight HLA Class I pseudogenes (HLA-H, -J, -K, -L, -P, -T, -V, -W), twelve HLA class II genes (HLA-DMA, -DMB, -DOA, -DOB, -DPA1, -DPA2, -DPB1, -DBP2, -DQA1, -DQB1, -DRA, -DRB1) and four non-HLA genes (MICA, MICB, TAP1, TAP2). Second, for each of the 650 samples (273 iPSCORE and 377 HipSci), reads which aligned to the target genes were extracted from each BAM file using coordinates from Gencode v19 (https://www.gencodegenes.org/releases/19.html) and SAMtools version 1.2 (Li et al., 2009). Third, for each sample the extracted reads were re-aligned to the collection of all genomic HLA sequences in the IPD-IMGT/HLA database (release 3.30.0) (Robinson et al., 2016), allowing each read to be aligned to multiple reference sequences using the ‘-a’ option in BWA-MEM. The expected read counts for each HLA type were obtained with HLA-VBSeq software (http://nagasakilab.csml.org/hla/). For each HLA gene, only the alleles with mean coverage ≥20% of the average coverage calculated over the whole genome were considered, and a target HLA genotype was determined as follows: 1) If no allele passed the threshold, then there were not enough reads aligned to correctly identify an HLA type, and hence no alleles were called; 2) If there was only one allele that passed the threshold, and the depth of coverage was two or more times greater than that of the threshold, then the HLA locus was considered to be homozygous for that HLA allele; 3) If there was only one allele that passed the threshold, and the depth of coverage was less than twice that of the threshold, then the allele was called heterozygous with the second allele not determined; 4) If there were two or more alleles that passed the threshold, the alleles were sorted based on the depth of coverage (from high to low), if the depth of coverage of the top allele was more than twice as that of the second one, then the HLA locus was called homozygous for the top allele; 5) If there were two or more alleles that passed the threshold, the alleles were sorted based on the depth of coverage (from high to low), if the depth of coverage of the top allele was less than twice that of the second one, then the two alleles with the highest coverage were selected as the HLA type. For HLA types at four-digit resolution (used in Figure 2—figure supplement 1D) the eight-digit resolution was converted by removing high digit values.

Determining recall, twin concordance, Mendelian inheritance concordance and fibroblast-iPSC concordance of HLA type detection

Request a detailed protocol

Reproducibility of HLA typing for each gene at eight-digit resolution was assessed by calculating: 1) recall, the fraction of determined HLA types for each gene; 2) HLA type concordance in 25 monozygotic twin pairs in the iPSCORE resource; 3) Mendelian inheritance concordance across 17 trios in the iPSCORE resource; and 4) HLA type concordance in 231 fibroblast-iPSC pairs. The recall was calculated independently for the 273 iPSCORE samples and the 377 HipSci samples as the number of determined HLA types divided by the total number of samples. For each gene, we calculated concordance as the fraction of HLA types that matched in the 25 monozygotic twin pairs in iPSCORE or the 231 fibroblast-iPSC paired genomes from the same individual in the HipSci resource. We calculated Mendelian inheritance concordance for each gene as the percentage of HLA types segregating in non-overlapping trios (i.e. in families with multiple trios, each individual was only used once for this analysis). If one or more HLA types were undetermined for a given pair/trio, we excluded the pair/trio during the concordance calculation.

Hardy-Weinberg equilibrium (HWE)

Request a detailed protocol

To investigate HWE, we performed a likelihood ratio test on the allele frequency of each HLA allele independently in 275 unrelated Caucasians in the iPSCORE and HipSci cohorts. We tested each HLA allele independently for several reasons: 1) calculating HWE on sites with tens or hundreds of alleles is highly computationally intensive and complicated (Eguchi and Matsuura, 1990; Graffelman et al., 2017; Guo and Thompson, 1992; Li et al., 2014); and 2) the MHC region is under high selective pressure, therefore the proportions between different alleles do not satisfy the assumptions for HWE (Hardy, 1908). All associations with p-values<1×10−6 were considered as significantly deviating from HWE and flagged in the HLA-type QTL analysis.

Comparison of observed HLA type frequencies with other populations

Request a detailed protocol

In the Allele Frequency Net Database (AFND) (González-Galarza et al., 2015), 3,556,301 people were genotyped for 12 out of the 30 genes we examined in the MHC region at two-digit resolution, 3,469,268 people (17 genes) at four-digit, 124,721 people (14 genes) at six-digit, and 10,212 people (seven genes) at eight-digit. This would result in testing 226 alleles at two-digit resolution (115 iPSCORE and 101 HipSci), 395 alleles at four-digit resolution (222 iPSCORE and 173 HipSci), 321 alleles at six-digit resolution (174 iPSCORE and 147 HipSci), and 42 alleles at eight-digit resolution (21 iPSCORE and 21 HipSci). We conducted the allele frequency comparative analysis using four-digit resolution (identified by collapsing eight-digit HLA types) to maximize the number of individuals in AFND, the number of genes and the number of alleles. The allele frequency of each HLA type in iPSCORE, HipSci and in the AFND was calculated by dividing the total number of individuals containing the given HLA type by the total number of people in each cohort and a correlation value was calculated after fitting a linear model.

Human iPSCs used in expression studies

Request a detailed protocol

For the expression studies, we used RNA-seq data for 446 iPSC samples from 361 individuals (215 iPSCORE and 146 HipSci, Supplementary file 2). For the 215 iPSCORE individuals, fibroblasts were reprogrammed, an iPSC clone obtained and RNA-seq data generated and processed as previously described (D'Antonio-Chronowska et al., 2019; D'Antonio et al., 2018; DeBoever et al., 2017; Panopoulos et al., 2017). Briefly, sequenced RNA-seq reads were aligned to the Gencode V.19 transcriptome using STAR (2.5.0a) (Dobin et al., 2013; Harrow et al., 2012). Gene expression levels (TPM) were calculated using RSEM (1.2.20) (Li and Dewey, 2011). For the 231 iPSC samples in HipSci (corresponding to 146 individuals), RNA-seq data were downloaded from the European Nucleotide Archive (ENA, https://www.ebi.ac.uk/ena, study PRJEB15299) and processed using the same pipeline used for the iPSCORE iPSC RNA-seq data.

Estimating HLA type-specific gene expression levels

Request a detailed protocol

The IMGT-IPD HLA database release 3.30.0 (Robinson et al., 2016), RSEM version 1.2.20 (Li and Dewey, 2011) contained cDNA sequences corresponding to each eight-digit resolution allele of the 24 HLA genes expressed in iPSCs. To minimize the number of variant differences between the RNA-seq reads mapping to the HLA genes, for each individual, we estimated HLA type-specific expression by replacing the canonical cDNA reference sequences of the 24 genes in Gencode v19 with cDNA sequences corresponding to the HLA types detected in the individual using HLA-VBSeq (Nariai et al., 2015). Because individuals with heterozygous HLA types had two cDNA sequences associated with the same HLA gene, while homozygous individuals had only one, we determined HLA type-specific expression in two distinct ways. For heterozygous individuals, the TPM associated with each cDNA sequence was considered as HLA type-specific expression, while for homozygous individuals, the TPM was divided by two to obtain HLA type-specific expression levels.

HLA types associated with expression levels of cognate HLA gene

Request a detailed protocol

We compared the expression level distributions of each HLA type against all other alleles from the cognate HLA gene (Figure 2; Supplementary file 5). Nominal p-values were calculated using Mann-Whitney U test. To FDR-correct, we used the p.adjust function in R with method = ‘BH’ (Benjamini-Hochberg). All tests with adjusted p-value<0.05 were considered significant.

eQTL detection

Request a detailed protocol

For eQTL analyses, we used the WGS from 361 individuals (215 iPSCORE and 146 HipSci, Supplementary file 2) for which iPSC RNA-seq data were available. We filtered variants that: 1) had a call rate <95%; 2) had minor allele frequencies < 1%; 3) deviated from HWE in the 275 unrelated Caucasian individuals (p<1.0×10−6). We obtained 45,245 filtered variants (40,668 SNPs and 4577 small indels) in the MHC region (plus 1 Mb upstream and downstream; exact coordinates: chr6:29640168–33115544) using bcftools (Li et al., 2009) and decomposed multiallelic genotypes using vt (Tan et al., 2015). Out of the 383 Gencode v19 genes in the MHC locus (chr6:29640168–33115544), 146 were expressed genes (TPM ≥2 in at least 10 samples), including 24 HLA genes. eQTL analysis requires expression levels associated with each gene; therefore, for the 24 HLA genes, we transformed the estimated HLA type-specific expression levels to gene-level expression by summing the expression of both HLA types in each individual. Expression levels of all 146 genes were quantile-normalized across all individuals. Normalized gene expression levels were adjusted for sex, age, batch (iPSCORE and HipSci), iPSC passage, read depth, read length, ten PEER factors (Stegle et al., 2010) and, to account for ethnicity, ten principal components calculated on the genotypes of 90,081 common SNPs in linkage equilibrium that we previously used to estimate the ethnicity of iPSCORE individuals (Panopoulos et al., 2017). To adjust for potential unknown sources of noise in gene expression, the ten PEER factors were calculated on genome-wide expression levels of 20,595 genes (with TPM ≥2 in at least 10 samples), rather than only on the expression levels of genes in the MHC region. We performed eQTL analysis for the 146 expressed genes and the 45,245 variants using a linear mixed model (lmekin from the coxme package in R) with a kinship matrix (to take population structure into account, including multiple ethnicities, members of the same family and twin pairs) calculated on 90,081 common SNPs distributed across the genome and in linkage equilibrium.

In order to detect significant eQTLs for each eGene, and then detect eGenes at genome-wide significance levels, we employed a three-step procedure to perform hierarchical multiple testing correction (Huang et al., 2018): 1) locally adjusted p-value at gene-level: nominal p-values were Bonferroni-corrected for the number of tests performed for each gene (45,245, Supplementary file 6A); the variant with the minimum adjusted p-value was considered as lead eQTL (Supplementary file 6B); 2) genome-wide adjusted p-value of lead eQTLs: to FDR-correct across all expressed genes, we further adjusted the p-values of the lead eQTLs using the p.adjust function in R with method = ‘BH’ (Benjamini-Hochberg) and n = 20,595 (corresponding to the number of expressed genes genome-wide) with a threshold of 0.05; and 3) eQTLs were identified for each significant eGene as variants with a locally adjusted p-value from step 1 < 0.1.

To detect QTLs conditional on the lead eQTL, for each gene we repeated the eQTL adding the genotype from the lead eQTL as a covariate (Conditional 1). To detect eQTLs conditional on the top two independent eQTLs (Conditional 2), we repeated eQTL detection using the genotypes from the lead eQTL and the genotype from the lead Conditional 1 eQTL as covariates.

Enrichment of eQTLs for regulatory elements

Request a detailed protocol

To test if eQTLs were enriched for overlapping iPSC regulatory elements, we used the following approach: 1) For all 45,245 tested variants, we determined the overlapping chromatin mark in each of five Roadmap iPSCs; 2) We defined effect size deciles for all tested eQTLs (each eQTL-eGene pair identified by testing 45,245 variants in 146 expressed genes) and binned the eQTLs based on their effect size (Figure 4—figure supplement 3, Supplementary file 6); 3) For each Roadmap iPSC sample, we counted the number of eQTLs associated with each mark in each bin, creating a 10-by-15 matrix (10 bins and 15 chromatin marks); 4) For normalization, we divided each column by its average across all bins and log2-transformed each value (each chromatin mark was tested individually), in order to obtain enrichments centered around zero; and then 5) we plotted these enrichment values for each chromatin mark (Figure 4—figure supplement 3).

Associations between eQTLs and allele-specific expression in eGenes

Request a detailed protocol

We obtained allele-specific expression (ASE) data from our previous genome-wide eQTL analysis on 215 iPSCORE iPSCs (DeBoever et al., 2017) and defined the strength of ASE for each gene as the fraction of RNA transcripts that were estimated to originate from the allele with higher expression (hereto referred to as ‘allelic imbalance fraction’, AIF). We tested the distributions of AIF calculated using MBASED (Mayba et al., 2014) on all samples and performed two distinct analyses. First, we tested if the 56 eGenes were more likely to have ASE by comparing their AIF in the samples where their associated lead eQTL was heterozygous with the 90 expressed MHC genes that did not have eQTLs using Mann-Whitney U test (Figure 4—figure supplement 4). Next, we tested if eGenes in samples with a heterozygous lead eQTL had stronger AIF than samples homozygous for the lead eQTL using Mann-Whitney U test (Figure 4—figure supplement 5A). Finally, we examined all primary and conditional eQTLs (39,891 eQTL-eGene pairs tested, corresponding to 21,739 variants) to determine if eGenes in samples with heterozygous primary eQTLs had stronger AIF than samples homozygous for primary eQTLs (Figure 4—figure supplement 5B).

Associations between eQTLs and GWAS

Request a detailed protocol

To test the associations between eQTL signal and complex traits and disease, we downloaded the summary statistics of 4,083 GWAS experiments from the UK BioBank (http://www.nealelab.is/uk-biobank). We identified 40,998 variants that were both tested for eQTLs and tested in the UK BioBank. For each of the 56 eGenes and each GWAS experiment, we extracted all p-values associated with these variants and performed colocalization using the coloc.abf from the coloc package in R (Giambartolomei et al., 2014). Using this tool, for each tested eGene-GWAS trait pair, we obtained posterior probabilities (PP) for five hypotheses: 1) hypothesis 0 (H0): neither the eGene nor the GWAS have a genetic association in the tested region; 2) H1: only the eGene has a significant association in the tested region; 3) H2: only the GWAS trait has a significant association in the region; 4) H3: both the eGene and the GWAS trait have significant associations, but their causal variants are different; and 5) H4: both the eGene and the GWAS trait and share the same causal variant. The sum of all five posterior probabilities is one for each tested eGene-GWAS trait pair. All pairs with PP H4 >0.2 are reported in Supplementary file 8.

HLA-type eQTL detection

Request a detailed protocol

For each sample, the HLA types were assigned dosage values and converted to VCF file format. Dosage was assigned as follows: 0 = the sample did not harbor the analyzed allele; 0.5 = the sample was heterozygous for the analyzed allele; 1 = the sample was homozygous for the analyzed allele. We adjusted the expression levels of all 146 expressed using the same covariates described for the eQTL detection (sex, age, batch, iPSC passage, read depth, read length, ten PEER factors and ten principal components calculated on the genotypes of the 90,081 common SNPs in linkage equilibrium) and performed HLA-type eQTL analysis investigating associations between each single HLA type with allele frequency >1% (283 in total) and expression levels of each of the 146 expressed genes using a linear mixed model (lmekin from the coxme package in R) with a kinship matrix calculated on the 90,081 common SNPs in linkage equilibrium. P-values were adjusted for FDR by using the p.adjust function in R with method = ‘BH’ (Benjamini-Hochberg) on all combinations between the 146 expressed genes by 283 HLA types present in the 361 individuals (corresponding to 41,318 tests) with a threshold of 0.05. The 15 HLA types that were not in HWE (p<1×10−6, Figure 2—figure supplement 1) are flagged in Supplementary file 9. We initially analyzed self-associations (defined as HLA-type eQTL between an HLA type and the expression of its associated gene, for example HLA*01:01:01:01 and the expression of HLA-A) independently to examine associations between expression levels of each HLA gene and its HLA types. To determine if HLA-type eQTLs were independent from single-variant eQTLs, we repeated the HLA-type eQTL detection by conditioning on: 1) the lead primary eQTL (defined as ‘conditional 1’); 2) the lead primary and lead conditional 1 eQTLs (‘conditional 2’); and 3) on the lead primary, conditional 1 and conditional two lead eQTLs (‘conditional 3’).

Coexpression analysis

Request a detailed protocol

To test if the four clusters of HLA-type eQTLs (Figure 6) identified coexpressed genes associated with the same regulatory variants, we calculated the correlation between the expression levels of each pair of the 146 genes in the MHC region (Figure 7—figure supplement 1).

To determine coexpression between each pair of genes in the MHC region, we calculated a linear model as:

ei~ej+covariates

Where ei and ej represent the expression of the i-th and j-th genes in the MHC region and covariates are the same used for the eQTL analysis. Beta was used as measure of coexpression.

Accession numbers

Request a detailed protocol

Whole-genome sequencing data of 273 individuals in the iPSCORE cohort (Panopoulos et al., 2017) is publicly available through dbGaP: phs001325. RNA-seq data of 215 iPSCs in the iPSCORE cohort (DeBoever et al., 2017) is publicly available through dbGaP:phs000924. Whole-genome sequencing data of 377 individuals in the HipSci cohort is publicly available through EGA: PRJEB15299. RNA-seq data of 231 iPSCs in the HipSci cohort is publicly available through ENA: PRJEB7388. eQTL and HLA-type eQTL summary statistics are available at Figshare (https://doi.org/10.6084/m9.figshare.10084733.v1 and https://doi.org/10.6084/m9.figshare.10084736.v1).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
    The quest for better understanding of HLA-disease association: scenes from a road less travelled by
    1. J Holoshitz
    (2013)
    Discovery Medicine 16:93–101.
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77

Decision letter

  1. Mark I McCarthy
    Senior Editor; Genentech, United States
  2. Calliope Dendrou
    Reviewing Editor; University of Oxford, United Kingdom
  3. Vivian G Cheung
    Reviewer; HHMI, University of Michigan, United States
  4. Heather Cordell
    Reviewer; Newcastle University, United Kingdom
  5. Diogo Meyer
    Reviewer; Institute of Biosciences of the University of São Paulo, Brazil

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The MHC is the region of the human genome associated with the greatest number of different human diseases. Despite this, unravelling the relationship between genetic variation in the region and the expression of HLA and non-HLA genes with relevance to disease is a significant challenge. Here the authors have used deep sequencing of stem cells from 419 individuals to create an extensive map of regulatory variation in the MHC region, exploring genetic associations for over 4,000 different disease and traits, and uncovering links between 180 such associations with variation in the expression of specific genes. This study provides insights into the genetic architecture of the MHC region and sheds new light on the relatively underappreciated concept that non-coding, regulatory genetic variation in the MHC region can also contribute to disease risk. Therefore, this work has important implications for understanding the role of the MHC in human disease, and will be of interest to geneticists, immunologists, molecular biologists and clinical scientists.

Decision letter after peer review:

Thank you for submitting your article "In-depth genetic analysis of 6p21.3 reveals insights into associations between HLA types and complex traits and disease" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Mark McCarthy as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Vivian G Cheung (Reviewer #1); Heather Cordell (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

In this manuscript whole genome sequencing is employed to obtain HLA genotypes, which are then used to in a gene expression analysis to identify the regulatory variants affecting HLA and non-HLA genes in the major histocompatibility complex region of human chromosome 6. The study implies that HLA types can predict the expression of non-HLA genes, and identifies potential clusters of co-expressed genes. As the HLA region is highly polymorphic and many of its variants have been associated with immune-based and other human diseases, an accurate approach for determining HLA genotypes from whole genome sequencing data, and an improved understanding of the impact of variants within the region on modulating gene expression, are important.

The essential revisions detailed below provide the collated comments and queries from the reviewers that need to be addressed. The reviewers were in agreement with respect to the need to give the paper a more cohesive structure, defining one or a few key questions and focusing on the most novel findings. We think that the suggestions regarding the statistical methods and the employment of simulation-based approaches are helpful and should be carefully considered. All reviewers also commented on the need to provide further support for the eQTL findings. Although we do not consider additional experimental validation to be necessary, we would like an indication of whether the eQTLs found were shared with those reported in previous studies or other available data sets. Other types of evidence that can support the eQTLs, such as the localisation of relevant variants in regulatory elements or GWAS enrichment of the variants, should also be investigated and reported.

Essential revisions:

1) Ideas presented in the Introduction are not taken up in the conclusion. The Introduction refers to "altered self" and "mistaken identity" hypotheses. Although the latter is taken up in the proof of principle, these terms are not used again, making the Introduction and Discussion somewhat disjointed. More generally, the paper covers a broad array of issues, but lacks focus: it was unclear exactly what the main idea is. If the focus is on demonstrating the enrichment of sites with regulatory role, a concern is that with respect to comparisons across studies with power that is non-comparable. If the idea is to defend the "mistaken identity" hypothesis, this needs to be taken up more thoroughly in the Discussion, given that the idea is now built upon, almost exclusively, a case study of one locus and its putative HLA-eQTL.

2) The history of HLA typing need to be amended. The text states that "the development of array-based technologies provided methods to estimate HLA types at four-digit resolution". Four-digit resolution greatly predates array technologies, which were brought in late in the history of HLA typing as a means to impute four-digit alleles, given appropriate reference samples and SNP density. Similarly, the possibility of "incorporating synonymous variants and non-coding regulatory variants" is attributed to NGS-based methods (supported by a string of citations between 2012 and 2017) but (Sanger) sequencing and PCR-based with allele-specific primers describing these levels of variation have long been around (in fact, the data in the IMGT resource was generated by Sanger sequencing of cloned sequences). The critical step of the cited papers was in fact adapting NGS pipelines for HLA typing.

3) In Figure 4, there are 3 groups of genes (described in the subsection “Genes in the MHC region are enriched for SNP-eQTLs”): those who expression levels show association with 283 SNP alleles, those with 10 or fewer and those with more than 1,000 associations. The data shown in Figure 4 do not show the distributions of expression levels for each allele so it is not possible to tell whether the associations are driven by just a few samples. Can the authors replot some of the results with graphs where expression levels are on the y-axis, and the three genotype groups on the x-axis?

4) An important aspect of the study is the identification of "four groups of genes that have alleles with highly correlated expression". This is largely supported by the results in Figure 5, which show sets of genes with significant correlations with HLA genotype (i.e., HLA-eQTL in the paper's terminology). This is a potentially interesting finding. Why was this analysis carried out in the first place? What are prior results that motivate the expectation that HLA allele identity will be predictive of expression of neighboring genes? Is this a causal mechanism or are HLA genes marking a haplotype with other genetic signatures? If documenting co-expression is of interest, why not examine it directly, rather than mediated through associations to an HLA type?

5) In the treatment of HLA-eQTLs, an interesting idea that yielded useful results, the authors chose not to examine the role of HLA genotypes on the expression of other HLA genes in the MHC region. This examination should either be expanded upon or a justification for the lack of a more extensive analysis needs to be provided.

6) The entire discussion of how levels of variability at 2 and 8 digit relate occupies a fair portion of the paper. Although essentially correct, much of the pattern seen is naturally expected given the nesting of 8 digit definitions within 2 digits. Thus, large sections of the paper relating diversity, using correlation analyses (Figure 2) and discussions of findings could be omitted. The very definition of these levels makes the patterns found and described expected. This shows up again in the discussion of concordance, where the fact that the values are higher for 2 digits is presented a result. This is unavoidable, given the nesting, which constrain the results so that the 2-digit concordance can only be the same as or greater than that at 8-digits. This is a feature of the data, not a result.

7) The paper dedicates considerable space to discussing variation among alleles. The usage of "allele-specific expression" as employed here does not match the usual definition, which refers to the "comparing the expression of the alleles in a heterozygous individual (allele specific expression or ASE analysis)" (https://doi.org/10.1186/s12864-018-5181-0). What is being done is in fact a test for differences in expression among alleles, a different question. Figure 3 treats the data as "each point reporting one allele". However, in the text referring to this same analysis we read "We (…) compared the TPM distributions between all samples carrying a specific HLA allele and all the other samples." This is a different analysis, since the two alleles of heterozygous individuals are confounded, both contributing to the locus-level expression. Although sometimes employed, extrapolating an individual's expression to an allele which is carried is an approach that introduces an additional source of uncertainty, which can be dealt with by analyzing the allele-level expression (since the allelic sequence in present in the index, allele-level expression could be directly compared). The paper needs to clarify what was done and justify the choices.

8) Throughout the manuscript, the authors refer to the personalized component in the index used to estimate expression levels as "personalized transcripts". However, according to the Methods, cDNA sequences from the IMGT database were used instead of personalized isoform sequences. While this is a personalized index at the HLA level, it is not a personalized transcriptome-based index since transcriptomic diversity is not being assayed as in Gensterblum-Miller et al. (https://doi. org/10.4049/jimmunol.1701061).

9) The paper fails to acknowledge previous work. The statement "The inability to genetically interrogate the MHC region" fails to acknowledge work such as Lee et al. (https://doi.org/10.1093/bioinformatics/bty), Jensen et al. (10.1101/gr.218891.116), Aguiar et al. (10.1371/journal.pgen.1008091). Similarly, the paper states that high resolution data for high depth data has not been provided across the region, but the above papers do address this issue. Later, in Results, when discussing the use of personalized transcripts to determine gene expression levels the authors argue that "the expression levels of genes in this region measured by RNA-seq may be inaccurate due to large numbers of mismatches between sequence reads and reference sequence" but do not have results on this, and fail to cite papers such as Aguiar et al., 2019, which explores exactly the relationship between mismatches and accuracy of expression through simulations or more general treatments such as Panousis et al., 2014.

10) The paper presents several analyses to support the reliability of the HLA genotype calls. Initially, the study examines the SNP density interferes with alignment coverage. It then argues that values which are sequencing depths (not coverage) are comparable to other genome regions and thus claims this implies "high quality". While the info on coverage and depth is useful, it is not clear that these values alone can be treated as proxies of high quality.

11) The paper also uses concordance among twins to "examine the accuracy of the high-resolution HLA types". While concordance is useful to document if there are no experimental errors, it does not guarantee accuracy. For example, if there is a systematic error in the alignment procedure, the same and incorrect sample could be obtained in each case. Thus, this may be indicated as an assay to measure experimental replicability, but not accuracy.

12) The argument for an improvement in expression estimates using the personalized pipeline (subsection “Using personalized transcripts to determine gene expression levels”) is that class I estimates are higher and class II are lower, when the pipeline is used. It is unclear why this is support for greater accuracy. It would imply that the higher expression of class II with the standard reference index results in overestimates, but mismatches between that index and the individual's alleles don't make such a prediction. The fact that HLA class II genes have lower expression doesn't naturally lead to the prediction that the effect of a better index would be to lower expression. Accuracy would require some external measure of expression (which isn't available) or a simulation-based approach (as implemented in some previously cited studies).

13) The paper spends some time correlating genotype and allele frequencies, including a report of negative correlation between number of alleles and number of homozygous genotypes. The standard way to explore this is to test for deviations from Hardy-Weinberg proportions, something that can be done using various existing software packages.

14) The study also shows a correlation between allele frequencies in the two datasets. While this is comforting, it is a very coarse indicator, which at most could be reported as a supplementary result, but doesn't warrant inclusion in the body of the paper.

15) For the eQTL mapping, the authors have adjusted the expression levels for different batches, ethnicities etc, but have not controlled for the effects of population structure on the genotype data (see for example the use of PCA in Lappalainen et al., Nature, 2013). It would also be helpful to provide more details on how the samples were used in the analyses. More details on how the adjustments were made to ensure that ethnicity did not add unexpected structure to the data would be helpful. Additionally, only one of the MZ twins should be used – it is not clear that the twins were removed.

16) The authors should justify the use of the overly conservative Bonferroni correction in the eQTL mapping instead of the widely-used FDR correction in genome-wide studies. By querying the MHC region alone, the definition of significance becomes different from genome-wide studies, making contrasts between these (including the proportion of eGenes) not directly comparable. Wouldn't it be preferable to carry out a genome-wide analysis, and then use the distribution of p-values to estimate corrected values under a chosen FDR? This may give readers a better ability to place these findings in the context of previous eQTL studies.

17) The authors claim that there is an enrichment of eGenes on the MHC region by comparing with a previous genome-wide work which analyzed only a subset of the samples (215 of the 446 samples). Since sample size affects power to detect eGenes, the direct comparison is not appropriate. More generally, the notion that "in the MHC region a large fraction of expressed genes are eGenes" should be stated in that qualifies the identify of an eGene as contingent on the dataset. A gene can be an eGene in some studies and not others, as a function of sampling. Thus, the statement that "genes in the MHC region are significantly enriched for having their expression levels modulated by regulatory variants" is not supported by the analysis.

18) Further support for the eQTL findings. For some of the regulatory variants, are there any samples that are heterozygous and show that the two alleles are indeed expressed at the different levels? Were the eQTL variants enriched for localisation in regulatory sites? Were they enriched for GWAS signals and have the eGenes in question been previous linked to disease?

https://doi.org/10.7554/eLife.48476.sa1

Author response

Essential revisions:

1) Ideas presented in the Introduction are not taken up in the conclusion. The Introduction refers to "altered self" and "mistaken identity" hypotheses. Although the latter is taken up in the proof of principle, these terms are not used again, making the Introduction and Discussion somewhat disjointed. More generally, the paper covers a broad array of issues, but lacks focus: it was unclear exactly what the main idea is. If the focus is on demonstrating the enrichment of sites with regulatory role, a concern is that with respect to comparisons across studies with power that is non-comparable. If the idea is to defend the "mistaken identity" hypothesis, this needs to be taken up more thoroughly in the Discussion, given that the idea is now built upon, almost exclusively, a case study of one locus and its putative HLA-eQTL.

We thank the reviewers for their thoughtful and insightful comments and agree that the manuscript lacked focus. For these reasons, we have substantially re-written the manuscript and have focused on the following topics:

• Using whole genome sequence from 419 individuals to call eight-digit HLA types and RNA-seq data from 361 matched induced pluripotent stem cells (iPSCs) to create a comprehensive map of regulatory variants in the MHC region.

• Building on this regulatory map, we explore GWAS signals for 4,083 complex traits, and finding strong evidence of colocalization for 180 disease loci with expression quantitative trait loci (eQTLs).

• Showing that taking eight-digit HLA type haplotypes into account substantially increases power to identify eGenes over single-variant eQTL analysis.

• Examining the association between the 8.1 ancestral haplotype and delayed colonization in Cystic Fibrosis, and postulating that downregulation of RNF5 expression is the likely causal mechanism.

2) The history of HLA typing need to be amended. The text states that "the development of array-based technologies provided methods to estimate HLA types at four-digit resolution". Four-digit resolution greatly predates array technologies, which were brought in late in the history of HLA typing as a means to impute four-digit alleles, given appropriate reference samples and SNP density. Similarly, the possibility of "incorporating synonymous variants and non-coding regulatory variants" is attributed to NGS-based methods (supported by a string of citations between 2012 and 2017) but (Sanger) sequencing and PCR-based with allele-specific primers describing these levels of variation have long been around (in fact, the data in the IMGT resource was generated by Sanger sequencing of cloned sequences). The critical step of the cited papers was in fact adapting NGS pipelines for HLA typing.

To focus the manuscript and Introduction, we have removed this entire section describing the different resolutions of HLA typing. Instead we have added a new figure (Figure 1) that describes current standard HLA nomenclature

3) In Figure 4, there are 3 groups of genes (described in the subsection “Genes in the MHC region are enriched for SNP-eQTLs”): those who expression levels show association with 283 SNP alleles, those with 10 or fewer and those with more than 1,000 associations. The data shown in Figure 4 do not show the distributions of expression levels for each allele so it is not possible to tell whether the associations are driven by just a few samples. Can the authors replot some of the results with graphs where expression levels are on the y-axis, and the three genotype groups on the x-axis?

We thank the reviewers for this comment: to show that the eQTL associations are not driven by just a few samples we added Figure 4—figure supplement 2, which shows boxplots of the expression levels for each of the 56 eGenes divided by the three genotypes (homozygous reference, heterozygous and homozygous alternative) of the lead eQTL. These plots show that the eQTL signals are supported by many samples. We have also made the eQTL summary statistics publicly available

(https://figshare.com/s/c8533530204e1822c6f4 and https://figshare.com/s/eab712af0e4baeb99251).

4) An important aspect of the study is the identification of "four groups of genes that have alleles with highly correlated expression". This is largely supported by the results in Figure 5, which show sets of genes with significant correlations with HLA genotype (i.e., HLA-eQTL in the paper's terminology). This is a potentially interesting finding. Why was this analysis carried out in the first place? What are prior results that motivate the expectation that HLA allele identity will be predictive of expression of neighboring genes? Is this a causal mechanism or are HLA genes marking a haplotype with other genetic signatures? If documenting co-expression is of interest, why not examine it directly, rather than mediated through associations to an HLA type?

We thank the reviewers for this observations. We agree that the rationale for analyzing the four sets of genes with significant associations with HLA genotypes was not clearly explained.

Why was this analysis carried out in the first place?

We hypothesized that HLA types, which represent a haplotype of many regulatory and coding SNPs, could exert a greater effect on gene expression than single variants. This is now clearly stated in the Discussion and alluded to in the Results.

What are prior results that motivate the expectation that HLA allele identity will be predictive of expression of neighboring genes?

There were no prior results that motived us to hypothesize that this would be the case. However, we have since identified a previous study that conducted a similar analysis but on a considerably smaller scale. We have cited a recent paper by Lam et al.(Lam et al., 2017), which describes eQTLs in the MHC region in 31 LCL samples, in the Discussion.

Is this a causal mechanism or are HLA genes marking a haplotype with other genetic signatures?

We believe that these findings suggest that the HLA-type eQTLs capture regulatory variants associated with the expression of one or a few neighboring genes. Furthermore, unlike the eQTLs from the single-variant analysis that tend to span long distances most likely due to the strong LD in the MHC region, the HLA-type eQTLs capture regulatory variants that act over relatively short distances.

We have expanded the Results and Discussion to describe these associations in detail.

If documenting co-expression is of interest, why not examine it directly, rather than mediated through associations to an HLA type?

In response to this comment, we have added a new analysis section in the Results “HLA-type eQTLs capture regulatory variation associated with expression of neighboring genes” and in the Materials and methods section “Coexpression analysis”.

“We investigated the molecular underpinnings of the clustering of the most significant

HLA-type eQTLs within the four distinct genomic regions (Figure 6A): 1)

chr6:29640168-30152231, including HLA class I genes (HLA-A, HLA-F and HLA-G) and HLA class I pseudogenes; 2) chr6:31171745-31512238, including HLA class I genes

(HLA-B, HLA-C), MICA and MICB; 3) chr6: 32361463-32789609, including HLA class

II genes (HLA-DRA, HLA-DRB1, HLA-DQA1, HLA-DQB1 and HLA-DOB), TAP1 and

TAP2; and 4) chr6:32916389-33115544, including HLA class II genes (HLA-DPA1, HLA-DPB1, and HLA-DPB2). We observed that different eGenes in each cluster were typically associated with the genotypes of different HLA types, suggesting that the observed clustering was not due to co-expression (i.e., correlated expression levels) of neighboring HLA genes (Figure 7). There was one exception in cluster 3, with the expression of both HLA-DRB5 and HLA-DRB1 associated with many of the same HLA types (Figure 7C). For this reason, we examined the co-expression of the eGenes in each of the four clusters. Only five eGenes in cluster 3 (HLA-DRA, HLA-DRB5, HLA-DQA1, HLA-DQB1 and HLA-DRB1) displayed moderate correlated expression (Figure 7—figure supplement 1), despite having their expression levels largely associated with different HLA types of neighboring genes (Figure 7). These findings suggest that the HLA-type eQTLs capture regulatory variants associated with the expression of one or a few neighboring genes. Furthermore, unlike the eQTLs from the single-variant analysis that tend to span long distances most likely due to the strong LD in the MHC region (Figure 4), the HLA-type eQTLs capture regulatory variants that act over relatively short distances.”

5) In the treatment of HLA-eQTLs, an interesting idea that yielded useful results, the authors chose not to examine the role of HLA genotypes on the expression of other HLA genes in the MHC region. This examination should either be expanded upon or a justification for the lack of a more extensive analysis needs to be provided.

We apologize for the confusion: in our original submission, we only removed “self” associations, i.e., we did not analyze the role of HLA genotypes on the same HLA gene (for example HLAA:01:01:01:01 and HLA-A expression), but we investigated the associations between the HLA type of one HLA gene and the expression of other HLA genes (for example HLA-A:01:01:01:01 and HLA-B expression). In the new version of the manuscript, we have added the self-association analysis (for example HLA-A:01:01:01:01 and HLA-A expression) in the Results section “HLAtype eQTLs associated with the expression of 114 eGenes in the MHC region” and in Supplementary file 9. We have tried to make the text in the Results and Materials and methods describing this analysis clearer.

6) The entire discussion of how levels of variability at 2 and 8 digit relate occupies a fair portion of the paper. Although essentially correct, much of the pattern seen is naturally expected given the nesting of 8 digit definitions within 2 digits. Thus, large sections of the paper relating diversity, using correlation analyses (Figure 2) and discussions of findings could be omitted. The very definition of these levels makes the patterns found and described expected. This shows up again in the discussion of concordance, where the fact that the values are higher for 2 digits is presented a result. This is unavoidable, given the nesting, which constrain the results so that the 2-digit concordance can only be the same as or greater than that at 8-digits. This is a feature of the data, not a result.

We agree that the discussion about the differences between two-digit and eight-digit resolutions is superfluous and contributes to the lack of focus that the reviewers pointed out in Comment 1. Therefore, we have removed this analysis from the manuscript including the Results and Materials and methods, and all related figures.

7) The paper dedicates considerable space to discussing variation among alleles. The usage of "allele-specific expression" as employed here does not match the usual definition, which refers to the "comparing the expression of the alleles in a heterozygous individual (allele specific expression or ASE analysis)" (https://doi.org/10.1186/s12864-018-5181-0). What is being done is in fact a test for differences in expression among alleles, a different question. Figure 3 treats the data as "each point reporting one allele". However, in the text referring to this same analysis we read "We (…) compared the TPM distributions between all samples carrying a specific HLA allele and all the other samples." This is a different analysis, since the two alleles of heterozygous individuals are confounded, both contributing to the locus-level expression. Although sometimes employed, extrapolating an individual's expression to an allele which is carried is an approach that introduces an additional source of uncertainty, which can be dealt with by analyzing the allele-level expression (since the allelic sequence in present in the index, allele-level expression could be directly compared). The paper needs to clarify what was done and justify the choices.

We have greatly reduced the space in the text devoted to discussing expression differences among HLA types (alleles). We have eliminated all discussion of comparing alignments to reference versus cDNA sequences specific for the HLA types from the Results. We moved some of this to the supplementary figures.

We apologize for the confusion in the description of the analysis about obtaining expression levels from alignments to cDNA sequences specific for the HLA types. We have edited the Materials and methods to clarify that we treated homozygous and heterozygous individuals in different ways:

“The IMGT-IPD HLA database release 3.30.0(Robinson et al., 2016), RSEM version 1.2.20(Li and Dewey, 2011) contained cDNA sequences corresponding to each eight-digit resolution allele of the 24 HLA genes expressed in iPSCs. […] For heterozygous individuals, the TPM associated with each cDNA sequence was considered as HLA type-specific expression, while for homozygous individuals, the TPM was divided by two to obtain HLA type-specific expression levels.”

8) Throughout the manuscript, the authors refer to the personalized component in the index used to estimate expression levels as "personalized transcripts". However, according to the Materials and methods, cDNA sequences from the IMGT database were used instead of personalized isoform sequences. While this is a personalized index at the HLA level, it is not a personalized transcriptome-based index since transcriptomic diversity is not being assayed as in Gensterblum-Miller et al. (https://doi. org/10.4049/jimmunol.1701061).

We agree that “personalized transcripts” is not the correct term to use. As discussed in response to comment 7, we have eliminated all discussion of comparing alignments to reference versus cDNA sequences specific for the HLA types from the Results. We referenced the paper by Gensterblum-Miller et al. in the Results section “Eight-digit HLA types associated with expression of cognate HLA gene”

9) The paper fails to acknowledge previous work. The statement "The inability to genetically interrogate the MHC region" fails to acknowledge work such as Lee et al. (https://doi.org/10.1093/bioinformatics/bty), Jensen et al. (10.1101/gr.218891.116), Aguiar et al. (10.1371/journal.pgen.1008091). Similarly, the paper states that high resolution data for high depth data has not been provided across the region, but the above papers do address this issue. Later, in Results, when discussing the use of personalized transcripts to determine gene expression levels the authors argue that "the expression levels of genes in this region measured by RNA-seq may be inaccurate due to large numbers of mismatches between sequence reads and reference sequence" but do not have results on this, and fail to cite papers such as Aguiar et al., 2019, which explores exactly the relationship between mismatches and accuracy of expression through simulations or more general treatments such as Panousis et al., 2014.

We thank the reviewers for these comments. As part of focusing the paper we have removed all sections to which comment 9 is concerning. Nevertheless we have included references for: Jensen et al.(Jensen et al., 2017), Lee et al.(Lee et al., 2018), Panousis et al.(Panousis et al., 2014) and Aguiar et al.(Aguiar et al., 2019).

10) The paper presents several analyses to support the reliability of the HLA genotype calls. Initially, the study examines the SNP density interferes with alignment coverage. It then argues that values which are sequencing depths (not coverage) are comparable to other genome regions and thus claims this implies "high quality". While the info on coverage and depth is useful, it is not clear that these values alone can be treated as proxies of high quality.

In our original manuscript, we had concluded that WGS samples were high-quality because the read depth in the MHC region is comparable to the rest of the genome and that the SNP density is comparable to previous studies. This was based on the assumptions that, if the alignment were low-quality, the read depth would be lower than the rest of the genome and the SNP density distribution would be noisy and would not have the same shape as was previously observed. We understand that these assumptions are not sufficient to claim that WGS had high quality.

To address this comment, we have changed the title of the first section of the Results to “Eight-digit HLA typing” and we removed any reference to quality. We have modified the take-home message of this section to the fact that the alignment of all WGS samples is comparable to the rest of the genome and that the SNP density is similar to previous studies(Norman et al., 2017).

11) The paper also uses concordance among twins to "examine the accuracy of the high-resolution HLA types". While concordance is useful to document if there are no experimental errors, it does not guarantee accuracy. For example, if there is a systematic error in the alignment procedure, the same and incorrect sample could be obtained in each case. Thus, this may be indicated as an assay to measure experimental replicability, but not accuracy.

We agree that concordance is not a measure of accuracy. When we describe high concordance we no longer refer to this as “accuracy” but rather “replicability” and have changed the section title in the Results section to “HLA types have high recall rates and reproducibility”.

12) The argument for an improvement in expression estimates using the personalized pipeline (subsection “Using personalized transcripts to determine gene expression levels”) is that class I estimates are higher and class II are lower, when the pipeline is used. It is unclear why this is support for greater accuracy. It would imply that the higher expression of class II with the standard reference index results in overestimates, but mismatches between that index and the individual's alleles don't make such a prediction. The fact that HLA class II genes have lower expression doesn't naturally lead to the prediction that the effect of a better index would be to lower expression. Accuracy would require some external measure of expression (which isn't available) or a simulation-based approach (as implemented in some previously cited studies).

We have eliminated most of the discussion on personalized expression from the manuscript. See answers to comments 7 and 8 for a description as to how we modified the manuscript.

13) The paper spends some time correlating genotype and allele frequencies, including a report of negative correlation between number of alleles and number of homozygous genotypes. The standard way to explore this is to test for deviations from Hardy-Weinberg proportions, something that can be done using various existing software packages.

We agree with the reviewers and have removed all the analyses describing correlations between genotypes and allele frequencies.

We have added the following to the Results:

“To estimate the accuracy of the 526 called HLA alleles we tested for Hardy-Weinberg equilibrium (HWE) in the unrelated individuals and found that the vast majority (511, 97.1%) were in HWE (p > 1 x 10-6, Figure 2—figure supplement 1A-C), consistent with the reproducibility values in twin pairs (96.7%) and in fibroblast-iPSC pairs (95.7%). […] Overall, HLA-VBSeq(Nariai et al., 2015) was able to accurately detect eight-digit resolution HLA alleles in the 650 WGS samples as evidenced by high recall rates, reproducibility and allele frequencies consistent with Hardy-Weinberg law.”

And added the following to the Materials and methods:

“Hardy-Weinberg equilibrium (HWE)

To investigate HWE, we performed a likelihood ratio test on the allele frequency of each HLA allele independently in 275 unrelated Caucasians in the iPSCORE and HipSci cohorts. […] All associations with p-values < 1 x 10-6 were considered as significantly deviating from HWE and flagged in the HLA-type QTL analysis.”

14) The study also shows a correlation between allele frequencies in the two datasets. While this is comforting, it is a very coarse indicator, which at most could be reported as a supplementary result, but doesn't warrant inclusion in the body of the paper.

We agree with the reviewers and have removed the figure showing the correlation between allele frequencies in the two datasets from the manuscript.

15) For the eQTL mapping, the authors have adjusted the expression levels for different batches, ethnicities etc, but have not controlled for the effects of population structure on the genotype data (see for example the use of PCA in Lappalainen et al., Nature, 2013). It would also be helpful to provide more details on how the samples were used in the analyses. More details on how the adjustments were made to ensure that ethnicity did not add unexpected structure to the data would be helpful. Additionally, only one of the MZ twins should be used – it is not clear that the twins were removed.

We agree that the use of a linear model (such as MatrixEQTL) to test for eQTLs is not appropriate when tested individuals are related. Therefore, we repeated the eQTL analysis (both on single variants and HLA types) using a linear mixed model that includes a kinship matrix to take population structure into account.

We have substantially re-written the Materials and methods to describe these analyses:

“eQTL detection

For eQTL analyses, we used the WGS from 361 individuals (215 iPSCORE and 146 HipSci, Supplementary file 2) for which iPSC RNA-seq data were available. […] We performed eQTL analysis for the 146 expressed genes and the 45,245 variants using a linear mixed model (lmekin from the coxme package in R) with a kinship matrix (to take population structure into account, including multiple ethnicities, members of the same family and twin pairs) calculated on 90,081 common SNPs distributed across the genome and in linkage equilibrium.”

We have changed Supplementary files 6, 7 and 9 with the new QTL results.

16) The authors should justify the use of the overly conservative Bonferroni correction in the eQTL mapping instead of the widely-used FDR correction in genome-wide studies. By querying the MHC region alone, the definition of significance becomes different from genome-wide studies, making contrasts between these (including the proportion of eGenes) not directly comparable. Wouldn't it be preferable to carry out a genome-wide analysis, and then use the distribution of p-values to estimate corrected values under a chosen FDR? This may give readers a better ability to place these findings in the context of previous eQTL studies.

We thank the reviewers for suggesting how we could improve the eQTL analysis. In response to Comment 15, we have changed the method to perform the eQTL analysis from a linear model to a linear mixed model to account for random genetic effects, which has allowed us to use all the iPSC samples (including multiple samples from the same individual and from twin pairs) and therefore has improved our eQTL detection power.

To address the reviewers’ comments, we improved the explanation of the eQTL detection method in the Materials and methods section “eQTL detection”.

“In order to detect significant eQTLs for each eGene, and then detect eGenes at genomewide significance levels, we employed a three-step procedure to perform hierarchical multiple testing correction (Huang et al., 2018): 1) locally adjusted p-value at gene-level: nominal p-values were Bonferroni-corrected for the number of tests performed for each gene (45,245, Supplementary file 7A); the variant with the minimum adjusted p-value was considered as lead eQTL (Supplementary file 7B); 2) genome-wide adjusted p-value of lead eQTLs: to FDRcorrect across all expressed genes, we further adjusted the p-values of the lead eQTLs using the p.adjust function in R with method = “BH” (Benjamini-Hochberg) and n = 20,595 (corresponding to the number of expressed genes genome-wide) with a threshold of 0.05; and 3) eQTLs were identified for each significant eGene as variants with a locally adjusted p-value from step 1 < 0.1.”

17) The authors claim that there is an enrichment of eGenes on the MHC region by comparing with a previous genome-wide work which analyzed only a subset of the samples (215 of the 446 samples). Since sample size affects power to detect eGenes, the direct comparison is not appropriate. More generally, the notion that "in the MHC region a large fraction of expressed genes are eGenes" should be stated in that qualifies the identify of an eGene as contingent on the dataset. A gene can be an eGene in some studies and not others, as a function of sampling. Thus, the statement that "genes in the MHC region are significantly enriched for having their expression levels modulated by regulatory variants" is not supported by the analysis.

We agree that it would be misleading to compare the number of eQTLs across studies, where many other factors (most importantly, sample size) play important roles. Therefore, we have removed this statement from the text (Introduction, Results and Discussion).

18) Further support for the eQTL findings. For some of the regulatory variants, are there any samples that are heterozygous and show that the two alleles are indeed expressed at the different levels? Were the eQTL variants enriched for localisation in regulatory sites? Were they enriched for GWAS signals and have the eGenes in question been previous linked to disease?

We thank the reviewers for these suggestions. To respond to this comment we performed three additional analyses (explained below).

For some of the regulatory variants, are there any samples that are heterozygous and show that the two alleles are indeed expressed at the different levels?

To respond to this comment, we performed enrichment analyses to determine if eQTLs were associated with allele-specific expression (ASE):

1) We tested if the 56 eGenes were more likely to have ASE than the other expressed genes in the MHC region. We obtained ASE data from our previous eQTL analysis on 215 iPSCORE iPSCs(DeBoever et al., 2017) and tested the distributions of major allele frequencies calculated using MBASED(Mayba et al., 2014) on all samples between the 56 eGenes and the other 90 genes expressed in the MHC region. We found that eGenes were significantly more likely to have ASE [mean allelic imbalance fraction (AIF) = 67.3%, compared with 61.9% in the other genes, p = 4.38 x 10-64, Mann-Whitney U test, Figure 4—figure supplement 4, Supplementary file 7].

2) Next, we investigated if eGenes in samples with a heterozygous lead eQTL had stronger AIF than samples homozygous for the eQTL using Mann-Whitney U test. We observed that eGenes in samples with heterozygous eQTLs (primary and conditional) have stronger ASE compared with eGenes in samples homozygous for eQTLs (39,891 eQTLs in total examined, corresponding to 21,739 variants; Figure 4—figure supplement 5A, B).

We describe these analyses in the Results (“eQTLs associated with the expression of 56 eGenes in the MHC region”), in the Materials and methods (“Associations between eQTLs and allele-specific expression in eGenes”), Figure 4—figure supplements 4 and 5.

Were the eQTL variants enriched for localisation in regulatory sites?

We obtained 15 ChromHMM(Ernst and Kellis, 2012) chromatin states in five iPSC lines from the Roadmap Epigenomics Consortium(Roadmap Epigenomics et al., 2015) and tested if the effect size of eQTLs changed between states. Intuitively, variants in regulatory elements are more likely to be associated with gene expression, therefore have stronger effect size. To perform this analysis, we used the following approach:

1) For all 45,245 tested variants, we determined the overlapping chromatin mark in each of five Roadmap iPSCs;

2) We defined effect size deciles for all tested eQTLs (each eQTL-eGene pair identified by testing 45,245 variants in 146 expressed genes) and binned the eQTLs based on their effect size (Figure 4—figure supplement 3, Supplementary file 6);

3) For each Roadmap iPSC sample, we counted the number of eQTLs associated with each mark in each bin, creating a 10-by-15 matrix (10 bins and 15 chromatin marks);

4) For normalization, we divided each column by its average across all bins and log2-transformed each value (each chromatin mark was tested individually), in order to obtain enrichments centered around zero;

5) We plotted these enrichment values for each chromatin mark (Figure 4—figure supplement 3).

Figure 4—figure supplement 3 shows that, in all Roadmap iPSC samples, the eQTLs in the MHC region with the largest effect sizes were enriched for overlapping regulatory elements active in iPSCs including active transcription start sites (TSS), regions flanking active TSS, UTRs (defined as regions transcribed at 5’ and 3’) and transcribed regions.

We added a description of this analysis in the Results section “eQTLs associated with the expression of 56 eGenes in the MHC region”, in the Materials and methods section “Enrichment of eQTLs for regulatory elements” and in Figure 4—figure supplement 3.

Were they enriched for GWAS signals and have the eGenes in question been previous linked to disease?

We changed our approach to determine associations between eQTLs and complex traits to determine whether eQTL signals in the MHC region, rather than single variants, were associated with any of the 4,083 traits described in the UK BioBank. Using this approach, we used a colocalization approach(Giambartolomei et al., 2014) to test all the eQTL signals with each trait in the UK BioBank and observed 180 associations with PPA > 0.8 (Supplementary file 8, Figure 5). These results demonstrate the utility of our data to serve as a resource to narrow down the location of causal regulatory variants underlying known associations between genes in the MHC region and complex traits, as well as to identify biologically meaningful novel associations.

We have added a section in the Results (“Regulatory variants in the MHC region play important roles in complex traits”), in the Materials and methods (“Associations between eQTLs and GWAS”), Supplementary file 8 and Figure 5.

https://doi.org/10.7554/eLife.48476.sa2

Article and author information

Author details

  1. Matteo D'Antonio

    1. Institute for Genomic Medicine, University of California, San Diego, San Diego, United States
    2. Department of Pediatrics, Rady Children’s Hospital, University of California, San Diego, San Diego, United States
    Contribution
    Formal analysis, Investigation, Methodology, Writing—original draft
    Contributed equally with
    Joaquin Reyna
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5844-6433
  2. Joaquin Reyna

    1. Department of Pediatrics, Rady Children’s Hospital, University of California, San Diego, San Diego, United States
    2. Biomedical Sciences Graduate Program, University of California, San Diego, La Jolla, United States
    Contribution
    Data curation, Software, Methodology, Writing—original draft
    Contributed equally with
    Matteo D'Antonio
    Competing interests
    No competing interests declared
  3. David Jakubosky

    1. Biomedical Sciences Graduate Program, University of California, San Diego, La Jolla, United States
    2. Bioinformatics and Systems Biology Graduate Program, University of California, San Diego, San Diego, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  4. Margaret KR Donovan

    1. Bioinformatics and Systems Biology Graduate Program, University of California, San Diego, San Diego, United States
    2. Department of Biomedical Informatics, University of California, San Diego, San Diego, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  5. Marc-Jan Bonder

    European Molecular Biology Laboratory, European Bioinformatics Institute, Cambridge, United Kingdom
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  6. Hiroko Matsui

    Institute for Genomic Medicine, University of California, San Diego, San Diego, United States
    Contribution
    Data curation
    Competing interests
    No competing interests declared
  7. Oliver Stegle

    European Molecular Biology Laboratory, European Bioinformatics Institute, Cambridge, United Kingdom
    Contribution
    Supervision
    Competing interests
    No competing interests declared
  8. Naoki Nariai

    Department of Pediatrics, Rady Children’s Hospital, University of California, San Diego, San Diego, United States
    Present address
    Illumina, Inc, San Diego, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Methodology
    Competing interests
    No competing interests declared
  9. Agnieszka D'Antonio-Chronowska

    1. Institute for Genomic Medicine, University of California, San Diego, San Diego, United States
    2. Department of Pediatrics, Rady Children’s Hospital, University of California, San Diego, San Diego, United States
    Contribution
    Resources
    Competing interests
    No competing interests declared
  10. Kelly A Frazer

    1. Institute for Genomic Medicine, University of California, San Diego, San Diego, United States
    2. Department of Pediatrics, Rady Children’s Hospital, University of California, San Diego, San Diego, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    kafrazer@ucsd.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6060-8902

Funding

California Institute for Regenerative Medicine (GC1R-06673)

  • Kelly A Frazer

National Institutes of Health (HG008118)

  • Kelly A Frazer

National Institutes of Health (HL107442)

  • Kelly A Frazer

National Institutes of Health (DK105541)

  • Kelly A Frazer

National Institutes of Health (DK112155)

  • Kelly A Frazer

National Science Foundation (CMMI division award 1728497)

  • Kelly A Frazer

National Institutes of Health (T15LM011271)

  • Margaret KR Donovan

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Ivan Carcamo-Orive and the other members of the i2QTL Consortium for their comments. This work was supported in part by a California Institute for Regenerative Medicine grant GC1R-06673-B, National Science Foundation CMMI division award 1728497, and NIH grants HG008118, DK105541 and DK112155. MKRD was supported by T15LM011271. These funding agencies played no role in the design or conclusions of this study. The iPSCORE iPSC lines are available; please contact Dr. Kelly A Frazer.

Senior Editor

  1. Mark I McCarthy, Genentech, United States

Reviewing Editor

  1. Calliope Dendrou, University of Oxford, United Kingdom

Reviewers

  1. Vivian G Cheung, HHMI, University of Michigan, United States
  2. Heather Cordell, Newcastle University, United Kingdom
  3. Diogo Meyer, Institute of Biosciences of the University of São Paulo, Brazil

Publication history

  1. Received: May 15, 2019
  2. Accepted: November 19, 2019
  3. Accepted Manuscript published: November 20, 2019 (version 1)
  4. Version of Record published: December 10, 2019 (version 2)

Copyright

© 2019, D'Antonio et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 856
    Page views
  • 133
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

  1. Further reading

Further reading

    1. Computational and Systems Biology
    2. Epidemiology and Global Health
    Isa Kristina Kirk et al.
    Research Article
    1. Cancer Biology
    2. Computational and Systems Biology
    Adam C Palmer et al.
    Research Article Updated