Immune diversity sheds light on missing variation in worldwide genetic diversity panels.

Defining worldwide human genetic variation is a critical step to reveal how genome plasticity contributes to disease. Yet, there is currently no metric to assess the representativeness and completeness of current and widely used data on genetic variation. We show here that Human Leukocyte Antigen (HLA) genes can serve as such metric as they are both the most polymorphic and the most studied genetic system. As a test case, we investigated the 1,000 Genomes Project panel. Using high-accuracy in silico HLA typing, we find that over 20% of the common HLA variants and over 70% of the rare HLA variants are missing in this reference panel for worldwide genetic variation, due to undersampling and incomplete geographical coverage, in particular in Oceania and West Asia. Because common and rare variants both contribute to disease, this study thus illustrates how HLA diversity can detect and help fix incomplete sampling and hence accelerate efforts to draw a comprehensive overview of the genetic variation that is relevant to health and disease.

Introduction

Defining worldwide human genetic variation at a genomewide level promises breakthroughs to reconstruct human evolution [1] and to define how genetic variation contributes to disease, hence paving the way towards precision medicine [2]. By decreasing sequencing costs by orders of magnitude [3], Next-Generation sequencing (NGS) facilitated the development of large scale genome sequencing projects [4-6] aimed at fulfilling this goal. Notably, the 1,000 Genomes Project is currently the reference panel for human genetic diversity [7], having targeted 2,693 individuals from 26 worldwide populations (S1 Fig), with the aim to define common genetic variants whose frequency is higher than 1% in populations. Although completion of the 1,000 Genomes Project represents a milestone for human genomics, there is currently no metric to assess how representative such reference panels for genomewide variation are and sampling bias in particular remains a concern [8]. Here we propose to use Human Leukocyte Antigen (HLA) diversity as a benchmark to perform this assessment. Indeed, because of their critical role in immunity and reproduction, HLA genes evolved to become the most polymorphic human genes, with over 11,000 distinct protein variants encoded by the five most variable HLA loci (HLA-A, -B, -C, -DRB1 and -DQB1) [9], hence making them superb markers of human diversity [10]. Their impact on health and disease [11] and in transplantation medicine [12-14] in particular also ensured that HLA sampling is second to none, with hundreds of populations and millions of individuals studied worldwide [15]. The HLA system is thus both the most polymorphic and the most studied genetic system (Fig 1A), and HLA sampling is second to none and considerably more widespread than that of the 1,000 Genomes Project panel for example (Fig 1B). HLA variation can therefore provide a high resolution picture of the representativeness of reference panels for genomewide diversity [15].

Fig 1

HLA sampling (A) is much more widespread than that of the 1,000 Genomes Project (B), which represents the current reference panel of human genomewide diversity.

A. Circles represent the populations that were HLA typed. Colors and size of the circles correspond to the number of common alleles observed in each population (scales in the bottom left corner). B. Circles represent the region of origin for each of the 26 populations of the 1,000 Genomes Project. Each population has a color that corresponds to the associated geographical region: Europe (blue), Africa (yellow), Americas (red), South Asia (purple), and East Asia (green). A-B. Maps were generated using the ggplot2 package in R [16] and the world database.

Materials and methods

1,000 Genomes Project sequence data

Exome sequence data for all 2,693 individuals forming the 1,000 Genomes Project was obtained (S1 Fig) and then spilt into two groups: a first set of 992 exomes for which there was already available HLA typing data [17] and a second set of 1,701 exomes for which there was no HLA typing data. HLA typing was then performed with the Polymorphisms to Phenotypes for Medicine (PolyPheMe) software, an in silico solution to perform highly-accurate HLA typing from any type of NGS data (see S1 File). In this study, the first set of data was used to evaluate PolyPheMe’s performance and the second set was used for de novo, in silico HLA typing.

PolyPheMe HLA typing tool

The PolyPheMe software (Xegen, France) can perform high accuracy HLA typing for five genes (HLA-A, -B, -C, -DRB1, and -DQB1) for a wide range of sequence data, including whole exome or whole genome data (see S1 File for details about PolyPheMe). All analyses were performed with PolyPheMe v1.2 on exome sequences using the IMGT 3.28 database [9] as reference.

Evaluation of PolyPheMe performance

To investigate the performance of PolyPheMe, we first analyzed the subset of the 1,000 Genomes Project for which there was already available HLA typing data [17]. In silico typing was thus performed for the HLA-A, -B, -C, -DRB1, and -DQB1 genes for 992 individuals. 294 individuals were excluded for the HLA-DQB1 analysis due to lack of data likely linked to enrichment issues and a total of 9,332 types were generated.

The types obtained with PolyPheMe were then compared to those already established [17], and three categories were defined: 1/‘concordant’ when the result with PolyPheMe is precise (only one type isolated) and compatible with previously established types, 2/‘imprecise’ when the result with PolyPheMe includes multiple possible types but at least one of them is contained in the list of previously established types, and 3/‘discordant’ when the results are not compatible (S2 Fig).

Results initially labelled as‘imprecise’ or‘discordant’ were further investigated through manual inspection of the sequence reads. For 21 individuals, genomic DNA was also obtained from Coriell Institute and a new HLA typing was performed with a commercial kit (Holotype, Omixon) using NGS. Results were directly analyzed both with commercial Omixon software and PolyPheMe software. This verification step confirmed that 105 typing results initially identified as‘discordant’ and four typing results initially identified as imprecise were in fact correctly typed by PolyPheMe (S3 Fig).

After verification of the typing results, the performance of PolyPheMe reaches a precision of 99.7% (9,300/9,332) and only 0.3% of the types are erroneous (32/9,332) (S2 Fig). Analysis of these 32 cases shows that 29 of them represent unique cases and all 29 of these unique variants were observed in the panel in other individuals (S1 and S2 Tables). This shows that these errors are not biased towards particular variants. Instead, analysis of these cases indicate that the problem stems from under or overrepresentation of one of the two alleles of the locus.

After validating the performance of PolyPheMe, the complete panel of the 1,000 Genome Project was investigated (S4 Fig).

Defining common HLA variants

For each of the five target HLA genes, we established a list of the variants having an allele frequency >1% in at least one population in the Allele Frequency Net Database [15]. Populations with a low sample size can produce incorrect estimates of the real allele frequencies and, conversely, setting too high a limit for the minimum population size can eliminate populations and lead to a lack of coverage for certain geographical regions. Thus a threshold must be selected to reduce the noise produced by populations with low sampling and to maximize geographical coverage. Using this approach, we selected a minimum population size of 150 individuals (S5 Fig) and 340 distinct alleles were defined as common (S3 Table). As an alternative, we also performed an analysis where the number of times an allele was observed in populations is used as a measure of how common it is (S6 Fig).

Definition of the common HLA haplotypes

Common five-genes haploypes (HLA-A, -B, -C, -DRB1, and -DQB1) were defined as those having a frequency >1% in at least one population in the National Bone Marrow Program reference panel [18] (S7 Fig), the largest source of HLA haplotype data, with a total of 74,425 different haplotypes already observed in one or more population.

HLA haplotype prediction for the 1,000 Genomes Project panel

The PolyPheMe software integrates a module for haplotype prediction based on a previously described methodology [19]. Briefly, using the precise, five-locus HLA genotype, all haplotype combinations are first generated (S4 Table). Then, in selecting the most likely haplotype structures, a priority is given to combinations maximizing the frequency observed in a single population. If a genotype cannot be explained by a haplotype combination in a unique population, the combination that maximizes frequency across all populations is selected. If no haplotype combination can explain the observed genotype, the haplotype with the best frequency is selected and the other haplotype is predicted.

Results

As a first step to compare the HLA diversity observed in the 1,000 Genomes Project to that commonly observed in worldwide populations, we set to define the HLA types of all 2,693 individuals forming the 1,000 Genomes Project panel (S1 Fig) at the five most polymorphic HLA loci: HLA-A, -B, -C, -DRB1 and -DQB1. This analysis used the Polymorphisms to Phenotypes for Medicine (PolyPheMe) software, an in silico solution to perform highly-accurate HLA Typing from any type of NGS data (see Methods and S1 File). As an initial test, PolyPheMe’s performance was first assessed on the 992 individuals of the panel that had already been HLA typed using standard methodology [17]. This comparison revealed a 98.5% concordance between PolyPheMe’s results and previous HLA typing (S2 Fig), but after careful sequence inspection and targeted resequencing (S3 Fig), the accuracy was estimated to be higher, reaching 99.7% (S2 Fig, S1 Table). Following this validation step, PolyPheMe was then used to determine the HLA types for the remaining 1,701 individuals for whom there was no prior information (S4 Fig, S2 Table). Out of the 26,342 cases of in silico typing performed, only 132 types (0.5%) were not precisely resolved (i.e. defined without ambiguities at a two-field level) and for an additional 32 cases (0.1%), no alleles were identified (S4 Fig). Thus 26,178 precise types (99.4%) were obtained and represent 414 distinct variants: 84 for HLA-A, 160 for HLA-B, 67 for HLA-C, 31 for HLA-DQB1, and 72 for HLA-DRB1.

To assess how well the variants observed in the 1,000 Genomes Project panel represent common HLA variation, we used data from the Allele Frequency Net Database to define a common set of variants corresponding to those whose frequency is higher than 1% in populations with sufficient sampling (n>150 individuals, see Methods and S5 Fig). From the analysis of 256 such populations, 340 HLA variants were observed (S3 Table): comparison of these 340 variants to the 414 variants observed in the 1,000 Genomes Project panel shows that the two sets only share 266 variants in common (Fig 2). Thus 74 of the 340 common HLA variants (22%) are missing in the 1,000 Genomes Project panel. That the missing fraction is homogenously distributed across the five HLA genes investigated, ranging between 17 and 28%, shows that this number is not due to a sampling issue with a given gene but is rather a general representation of the common HLA variation missing in the 1,000 Genomes Project panel.

Fig 2

The HLA diversity in the 1000 Genomes Project panel only represents 78% of the expected diversity for the alleles with a frequency >1%.

The top part of the figure shows the % match (Y axis) between the expected HLA diversity at different frequency cutoffs (X axis) and the HLA diversity observed in the 1,000 Genomes Project panel. For each frequency cutoff, the number of expected alleles is displayed at the top of each histogram. The bottom part of the figure displays the same information on a locus by locus basis for the 1% cutoff.

To understand why 74 common HLA variants are not present in the 1,000 Genomes Project panel, we investigated their geographical and population distribution (Fig 3A). This analysis illustrates two limitations of the 1,000 Genomes Project sampling: missing populations and undersampled regions. The former refers to populations in regions of the world that are not represented in the 1,000 Genomes Project sampling. Indeed, ten and eleven of the missing variants originate for example from populations in Oceania and West Asia, respectively, two geographical regions that are not represented in the 1,000 Genomes Project sampling. Importantly, in those regions, the missing variants can be extremely common, with frequencies higher than 10% (Fig 3B). The latter, ‘undersampled regions’, corresponds to regions where the sampling does not fully grasp the population diversity of the region. This category includes both a lack of depth in the sampling in terms of the number of individuals studied and an insufficient geographical sampling within a region. For example, while five European populations are included in the 1,000 Genomes Project, fifteen common variants were missed: ten whose frequency is 1–2.5% and that might have been characterized with more sampling depth, and four whose frequency is >2.5% (2.5–10%) that were likely missed due to the limited geographical sampling within the region.

Fig 3

The common HLA alleles that are missing in the 1,000 Genome Project panel define a worldwide distribution with variable frequencies that range from low (1–2.5%) to high (>10%).

A. Worldwide distribution of the populations harboring alleles that are missing in the 1,000 Genomes Project. Colors and size of the circles correspond to the number of common alleles observed in each population (scale on the bottom left corner). The map was generated using the ggplot2 package in R [16] and the world database. B. Geographical distribution of the alleles that are missing in the 1,000 Genomes Project. For each region, the number of distinct alleles is given together with their maximum frequencies.

Thus, while a large fraction (78%) of the common HLA alleles are observed in the 1,000 Genomes Project panel, the HLA benchmark underlines two important areas of improvement. First, a significant fraction (7%) of the very common variants (frequency >10%) are missing in the 1,000 Genomes Project panel. Second, the representativeness of the panel drops dramatically when less common HLA variants are considered: to 33% for example if a larger set of expected alleles with frequency >0.1% is considered (Fig 2). A similar result is obtained if, instead of using a fixed frequency cutoff, the reference set of expected variants is set to include all the variants observed at least five times in populations, as only ~29% of them are observed in the 1,000 Genomes Project panel (S6 Fig). Likewise, assessing the 1,000 Genomes Project in terms of HLA haplotype diversity, a relevant unit for disease association studies, shows that only 37 of the 58 haplotypes (63.8%) whose frequency is >1% in the largest HLA haplotype database (National bone marrow program database [18]) are observed in the 1,000 Genomes Project panel (S7 Fig). That 20 of the 21 missing haplotypes are associated with native populations from the USA or from Oceanian populations, two groups already identified as ‘missing’ (Fig 3), shows that ‘missing populations’ are also the main source of missing haplotypes.

Discussion

Worldwide HLA diversity is thus a helpful metric to assess the representativeness of panels of human genetic diversity, as it can define both the common and rare genetic variation that is missing, with both types of variation having potential to impact expression and contribute to individual disease risk [20-22]. The approach isolated coverage issues in a reference panel as widely used as the 1,000 Genomes Project panel. While some of the results are intuitive, like the missing variation linked to regions of the world that were not sampled, HLA variation can help assess how much is missing and isolate missing information that is more difficult to predict, like undersampling in particular geographical regions. This information can in turn be used to improve existing panels and assess how relevant are new panels, a critical step to produce an exhaustive catalogue of human variation [23] and develop new therapeutic solutions [24]. This is important because while genomics is a critical step to understand human diversity [25], it is still failing due to historical bias in sampling that focused too much on individuals of European ancestry [8]. The 1,000 Genomes Project was the first to attempt to provide genomewide data on global human diversity but complementary projects that attempt to increase sample size [26] or focus on specific populations [6, 27] are also developing, together with approaches that gather and compile data from multiple sources [28]. Yet all these approaches to generate a better overview of human genetic variation require a metric to assess how successful they are and we show here that HLA variation is a compelling system to perform that assessment.

1,000 Genomes Project sampling.

The 1,000 Genomes Project sampling is given on a regional basis (top part) or on a population basis (bottom part). In both cases the number of individuals used to evaluate PolyPheMe’s performance is also given (performance test column).

(PDF)

Click here for additional data file.

PolyPheMe’s performance before and after the verification steps.

PolyPheMe’s performance was assessed on a subset of 992 individuals of the 1,000 Genomes Project panel who were already HLA typed using standard methodology. This figure summarizes the concordant/correct, imprecise, and discordant/erroneous results obtained before and after the validation step (S3 Fig).

(PDF)

Click here for additional data file.

Validation of PolyPheMe results.

PolyPheMe’s performance was assessed on a subset of 992 individuals of the 1,000 Genomes Project panel who were already HLA typed using standard methodology. 192 HLA types out of 9,332 were initially characterized either as imprecise or discordant comparing to existing results (S2 Fig). For 109 of these 192 cases, we could show that the PolyPheMe results were correct: this figure summarizes the validation steps used for those 109 cases.

(PDF)

Click here for additional data file.

In silico HLA typing of the 1,000 Genomes Project panel by PolyPheMe.

This figure summarizes the total number of in silico HLA types realized by PolyPheMe on a locus-by-locus basis. For each locus, the number of precise and imprecise types is given, together with the number of cases for which no type was obtained.

(PDF)

Click here for additional data file.

Impact of the population size threshold (horizontal axis) on the number of populations used (vertical axis on the left, blue line), on the number of common HLA variants defined (vertical axis on the left, orange line), and on the fraction of common HLA variants observed in the 1,000 Genomes Project panel (vertical axis on the right, grey line).

A threshold of at least 150 individuals per population was used for the final analysis (dotted line).

(PDF)

Click here for additional data file.

More than 70% of the HLA variants observed at least five times in populations are missing in the 1,000 Genomes Project panel.

For each occurrence threshold (minimum number of observations of a given allele in populations), the number of alleles of HLA-A, -B, -C,–DQB1, and–DRB1 observed with this minimum occurrence is given (line ‘1.’). The second line (‘2.’) indicates how many of the alleles listed in 1. are also observed in the 1,000 Genomes Project panel. Finally the third line corresponds to the fraction ‘2./1’. For example, 1,202 HLA-A, -B, -C,–DQB1, and–DRB1 alleles are observed at least five times in populations. Out of those 1,202 alleles, 354 are also observed in the 1,000 Genomes Project panel, which corresponds to a fraction of 0.29.

(PDF)

Click here for additional data file.

Common HLA haplotypes.

This figure summarizes the common HLA haplotypes and the populations that harbor them with frequency > 1%. The third column gives the number of observations of the haplotypes in the 1,000 Genomes Project panel. Population codes correspond to those described in S8 Fig. Two haplotypes that are not observed in the 1,000 Genomes Project are in bold because they carry common HLA variants (HLA-C*08:06 and HLA-DRB1*14:08).

(PDF)

Click here for additional data file.

Population codes for the haplotype analysis.

This figure provides the population codes used in the haplotype analysis (S7 Fig).

(PDF)

Click here for additional data file.

Evaluation of PolyPheMe performance.

This file provides the results of the analysis with PolyPheMe for the subset of the 1,000 Genomes Project for which there was already available HLA typing data. For 928 of the 992 individuals, the PolyPheMe results were ‘correct’ for all loci (first tab ‘Correct typing’). In some cases (marked by asterisks), the PolyPheMe typing was validated through manual verifications (see Methods). For 64 of the 992 individuals (tab ‘Others’), at least one HLA type was either imprecise (marked in orange) or erroneous (marked in red). The correct type was obtained through manual verifications (see Methods) and is given in green below the imprecise or erroneous type.

(XLSX)

Click here for additional data file.

Final HLA types for the complete 1,000 Genome Project panel.

This table gives the final HLA types in our analysis for the 2,693 donors of the 1,000 Genome Project panel.

(XLSX)

Click here for additional data file.

Common HLA variants.

This Table lists the 340 distinct HLA variants defined as ‘common’, i.e. whose frequency is higher than 1% in populations with sufficient sampling (see Methods).

(XLSX)

Click here for additional data file.

Five-locus HLA haplotypes predicted for the donors of the 1,000 Genome Project panel.

This Table provides the predicted five-locus HLA haplotypes (HLA-A, -B, -C, -DRB1, -DQB1) for the donors of the 1,000 Genome Project. For each donor, the two haplotypes are given, together with the population where each was observed with the highest frequency (population codes are given in S8 Fig).

(XLSX)

Click here for additional data file.

PolyPheMe HLA typing tool.

A description of PolyPheMe software is provided: general description, typing strategy, typing steps, results, and computing power requirements and availability.

(PDF)

Click here for additional data file.