PyHLA: tests for the association between HLA alleles and diseases
 Yanhui Fan^{1, 2, 3} and
 YouQiang Song^{1, 2}Email author
DOI: 10.1186/s1285901714960
© The Author(s). 2017
Received: 26 July 2016
Accepted: 27 January 2017
Published: 6 February 2017
Abstract
Background
Recently, several tools have been designed for human leukocyte antigen (HLA) typing using single nucleotide polymorphism (SNP) array and nextgeneration sequencing (NGS) data. These tools provide highthroughput and costeffective approaches for identifying HLA types. Therefore, tools for downstream association analysis are highly desirable. Although several tools have been designed for multiallelic marker association analysis, they were designed only for microsatellite markers and do not scale well with increasing data volumes, or they were designed for largescale data but provided a limited number of tests.
Results
We have developed a Python package called PyHLA, which implements several methods for HLA association analysis, to fill the gap. PyHLA is a tailormade, easy to use, and flexible tool designed specifically for the association analysis of the HLA types imputed from genomewide genotyping and NGS data. PyHLA provides functions for association analysis, zygosity tests, and interaction tests between HLA alleles and diseases. Monte Carlo permutation and several methods for multiple testing corrections have also been implemented.
Conclusions
PyHLA provides a convenient and powerful tool for HLA analysis. Existing methods have been integrated and desired methods have been added in PyHLA. Furthermore, PyHLA is applicable to small and large sample sizes and can finish the analysis in a timely manner on a personal computer with different platforms. PyHLA is implemented in Python. PyHLA is a free, open source software distributed under the GPLv2 license. The source code, tutorial, and examples are available at https://github.com/felixfan/PyHLA.
Keywords
HLA Association Interaction MultiallelicBackground
The human leukocyte antigen (HLA) loci on chromosome 6 (6p21.3) are the most polymorphic and genedense region of the human genome. HLA proteins play an important role in transplant rejection. Association of variants in the HLA region and infectious, autoimmune diseases and cancers has been established. Directly typing HLA using experiments is still laborious, expensive, and timeconsuming [1]. Several algorithms and pipelines, such as HLA*IMP:02 [2] and MGAPrediction [3] have been developed for HLA imputation using data from genomewide association studies (GWAS), whereas OptiType [4], HLAVBSeq [5] and HLAreporter [6] have been developed for HLA typing using data from nextgeneration sequencing (NGS) studies. All tools use HLA allele sequences from the IMGT/HLA database [7] as reference. These tools have provided us a costefficient, highthroughput approach for HLA typing by using the already available GWAS and NGS data.
Given the continuously increasing amounts of HLA types being generated, integrating the workflow for their downstream association analysis is highly desirable. Several existing tools, such as CLUMP [8], PyPop [9] and SKDM [10], can be used to analyze HLA types. These tools are not ideal for association analysis of HLA types inferred from GWAS and NGS data as they were designed for analyzing microsatellite markers or provided limited functions. In this study, we present PyHLA, a Pythonbased standalone tool, for the association analysis between diseases and HLA types inferred from GWAS and NGS data.
Implementation
Data summary (module 1)
Gene, allele and population level summary of the frequency can be produced in the case and control populations.
Association analysis (module 2)
It is a simple and easy way to implement methods for localization of susceptibility genes by comparing the allele frequencies between cases and controls from the same population. Usually, Pearson’s chisquared test or Fisher’s exact test is performed on a 2 × 2 contingency table, which contains the counts of minor and major alleles for a single locus in cases and controls. As the most polymorphic part of the human genome, HLA genes, such as HLAA, HLAB and HLAC, have several thousand known alleles [7]. PyHLA performs Pearson’s chisquared test or Fisher’s exact test on the 2 × 2 contingency table, which compares one allele with the other alleles grouped together.
If the HLAA gene has n common alleles in cases and controls, then n tests are performed. In each test, one allele is compared with the other n − 1 alleles grouped together. The allelic 2 × 2 contingency table for a specific HLA allele contains the counts of this allele and the counts of other n − 1 alleles in cases and controls. The dominant and recessive models assume that each allele is dominant and recessive to the other n − 1 alleles, respectively. The dominant 2 × 2 contingency table for a specific HLA allele contains the counts of individuals with and without the allele in cases and controls. The recessive 2 × 2 contingency table for a specific HLA allele contains the counts of individuals with and without two copies of this allele in cases and controls.
Pearson’s Chisquared test
where χ ^{ 2 } is the chisquare critical value with a degree of freedom equals to 1. O _{ i,j } and E _{ i,j } are the observed and expected frequencies of the cell in row i and column j, respectively. r is the number of rows and c is the number of columns; both are 2 for the 2 × 2 contingency table.
Fisher’s exact test
where O _{ i,j } is the observed frequency of the cell in row i and column j. r _{ i } and c _{ i } are the rows and columns of marginal totals, respectively. N is the grand total. P _{ cutoff } is the exact probability of obtaining such set of observed values. Then, the probability for all possible tables with the same marginal totals is calculated. The twosided p value for the Fisher’s exact test is calculated by summing all probabilities less than or equal to P _{ cutoff }.
Logistic and linear regression
Genotype coding for additive, dominant, and recessive models, with D being the risk allele
Genotype  Code  

Additive  Dominant^{a}  Recessive^{b}  
DD  2  1  1 
Dd  1  1  0 
dd  0  0  0 
Multiple testing correction
The p values can be adjusted by using the Bonferroni correction or false discovery rate (FDR) correction. The empirical pvalues can also be calculated using a permutation test, which randomly shuffles the phenotypes for individuals, while keeping the HLA alleles unchanged.
Amino acid association analysis
PyHLA can perform not only allele level association analysis but also the AA level association analysis. The aligned AA sequences were retrieved from the IMGT/HLA database [7]. Fisher’s exact test or Pearson’s chisquared test can be conducted to investigate AA occurrence that are significantly associated with a disease.
Zygosity test (module 3)
Three tests were performed to investigate homozygous, heterozygous, and zygosity associations. These three tests evaluate the frequency difference of subjects carrying the homozygous and heterozygous alleles and the absence of a particular allele/AA in cases and controls. An individual carrying two same alleles is considered homozygous in the allele level test. An individual carrying two identical alleles or an individual carrying two different alleles that code for the same AA residue is considered homozygous in the AA level test. Fisher’s exact test or Pearson’s chisquared test for a 2 × 2 contingency table can be used for the zygosity test.
Interaction test (module 4)
Interaction test performs eight tests for detecting the strongest association. These tests involve tests for independence, differential association, combined association, linkage disequilibrium, and interaction [10, 11]. Each of the eight tests is based on a 2 × 2 contingency table. Fisher’s exact test or Pearson’s chisquared test can be used for the interaction test.
Results
Since it is hard to find a publicly available real dataset. A simulated data set with 1000 cases and 1000 controls was used to demonstrate the usage of PyHLA. Detailed commands, inputs and outputs are available on https://github.com/felixfan/PyHLA/tree/master/demo.
Association test suggested that the two most significant alleles are HLAA*01:01 (P = 4.03E24, OR = 2.15) and HLADQB1*05:02 (P = 3.32E11, OR = 1.58). Zygosity test further showed that the susceptibility to disease between homozygote and heterozygote of these two alleles are different (P = 2.46E11 and P = 1.10E7 for HLAA*01:01 and HLADQB1*05:02, respectively.). The heterozygotes are individually associated with the disease (P = 1.21E19 and P = 6.60E8 for HLAA*01:01 and HLADQB1*05:02, respectively.). Finally, the interaction test suggested that HLAA*01:01 and HLADQB1*05:02 are in linkage disequilibrium in cases; their combined action is contributory to disease susceptibility.
Discussion
PyHLA provides an integrated pipeline for detecting HLA association in antigen (twodigit allele level), protein (fourdigit allele level) and AA levels. Zygosity tests will examine the homozygous, heterozygous, and zygosity associations once the associated alleles and AAs are identified. In addition, interaction test examines the independence, differential association, combined association, interaction, and linkage disequilibrium between two factors.
In addition to identifying alleles and AA residues that are significantly associated with the disease, PyHLA also tests whether the increased HLA homozygosity or heterozygosity contributes to the increased susceptibility to a disease. When several factors are associated with the disease, the interaction test identifies the strongest one between each pair of the two factors. The factor with the strongest association is more likely to be the causative factor that truly contributed to the disease [11].
In this work, Pearson’s chisquared test and Fisher’s exact test performed on a 2 × 2 contingency table were implemented in PyHLA. Linear and logistic regressions were also included to consider multiple covariates simultaneously.
Bonferroni adjustment and correction via FDR estimation are widely used for multiple testing corrections. Bonferroni correction assumes that all tests are independent and is conservative in genetic association analysis, whereas FDR is less stringent [12–15]. In addition, the empirical p values can also be calculated using the permutation test, which randomly shuffle the phenotypes for individuals while keeping the HLA alleles unchanged. The permutation test preserves the correlation structure among HLA alleles but requires a large number of random shuffles. Given that the number of HLA alleles is relatively smaller than the number of SNPs in the genome, the computing time and resources needed for the permutation test are significantly less. PyHLA can perform these analyses on a single modern personal computer in a timely manner.
Four chisquared tests were implemented in CLUMP [8] to test the association between disease and alleles at highly polymorphic loci, and Monte Carlo imputation was performed to estimate the significance level. CLUMP is mainly designed for analyzing microsatellite markers in qualitative trait studies (casecontrol study), but not in quantitative trait studies. CLUMP cannot perform residual level tests as well. SKDM [10] is specialized in casecontrol HLA analysis through the identification and subsequent dissection of AA association; it is not designed for quantitative studies. Only the Fisher’s exact test is available for association test, and only Bonferroni correction is available for multiple testing adjustment. PyPop [9] is designed to handle large sample sizes for population statistics, haplotype frequency estimation and linkage disequilibrium significance testing. PyHLA is designed to supplement and extend these existing software. PyHLA can handle both qualitative and quantitative trait studies in both amino acid level and different resolutions of allele levels. Both chisquared test and Fisher's exact test are implemented to test the association, and both Bonferroni correction and FDR are available for multiple testing adjustment. Monte Carlo imputation is also implemented to estimate the significance level. Moreover, logistic regression and linear regression implemented in PyHLA can also include covariates in the association analysis.
Conclusions
In summary, PyHLA is a userfriendly tool for HLA association analysis. Existing methods are integrated and additional desired methods are included in PyHLA. PyHLA is applicable to small and large sample sizes and can complete the analysis in a timely manner on a personal computer. PyHLA is designed for casecontrol studies. PyHLA is currently unable to analyze familybase datasets.
Abbreviations
 AA:

Amino acid
 FDR:

False discovery rate
 GWAS:

Genomewide association study
 HLA:

Human leukocyte antigen
 NGS:

Nextgeneration sequencing
Declarations
Acknowledgements
Not applicable.
Funding
This work was supported by grants from the National Natural Science Foundation of China (81271226), the RGC General Research Fund (17117715), and the Health and Medical Research Fund of Hong Kong Government (01121726, RRG08, and 04152256).
Availability of data and materials
Project name: PyHLA
Project homepage: https://github.com/felixfan/PyHLA
Operating systems: Windows, Linux, Mac OS X
Programming language: Python
License: GPLv2
Authors’ contributions
YQS and YHF conceived the idea and designed the research. YHF wrote and tested the software. YHF and YQS wrote the manuscript. Both authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Hsieh AR, Chang SW, Chen PL, Chu CC, Hsiao CL, Yang WS, Chang CC, Wu JY, Chen YT, Chang TC, et al. Predicting HLA genotypes using unphased and flanking singlenucleotide polymorphisms in Han Chinese population. BMC Genomics. 2014;15:81.View ArticlePubMedPubMed CentralGoogle Scholar
 Dilthey A, Leslie S, Moutsianas L, Shen J, Cox C, Nelson MR, McVean G. Multipopulation classical HLA type imputation. PLoS Comput Biol. 2013;9:e1002877.View ArticlePubMedPubMed CentralGoogle Scholar
 Li S, Wang H, Smith A, Zhang B, Zhang X, Schoch G, Geraghty D, Hansen JA, Zhao LP. Predicting multiallelic genes using unphased and flanking single nucleotide polymorphisms. Genet Epidemiol. 2011;35(2):85–92.View ArticlePubMedGoogle Scholar
 Szolek A, Schubert B, Mohr C, Sturm M, Feldhahn M, Kohlbacher O. OptiType: precision HLA typing from nextgeneration sequencing data. Bioinformatics. 2014;30:3310–6.View ArticlePubMedPubMed CentralGoogle Scholar
 Nariai N, Kojima K, Saito S, Mimori T, Sato Y, Kawai Y, YamaguchiKabata Y, Yasuda J, Nagasaki M. HLAVBSeq: accurate HLA typing at full resolution from wholegenome sequencing data. BMC Genomics. 2015;16 Suppl 2:S7.View ArticlePubMedPubMed CentralGoogle Scholar
 Huang Y, Yang J, Ying D, Zhang Y, Shotelersuk V, Hirankarn N, Sham PC, Lau YL, Yang W. HLAreporter: a tool for HLA typing from next generation sequencing data. Genome Med. 2015;7(1):1–12.View ArticleGoogle Scholar
 Robinson J, Halliwell JA, McWilliam H, Lopez R, Parham P, Marsh SGE. The IMGT/HLA database. Nucl Acids Res. 2013;41:D1222–1227.View ArticlePubMedGoogle Scholar
 Sham PC, Curtis D. Monte Carlo tests for associations between disease and alleles at highly polymorphic loci. Ann Hum Genet. 1995;59:97–105.View ArticlePubMedGoogle Scholar
 Lancaster AK, Single RM, Solberg OD, Nelson MP, Thomson G. PyPop update – a software pipeline for largescale multilocus population genomics. Tissue Antigens. 2007;69:192–7.View ArticlePubMedPubMed CentralGoogle Scholar
 Kanterakis S, Magira E, Rosenman KD, Rossman M, Talsania K, Monos DS. SKDM human leukocyte antigen (HLA) tool: a comprehensive HLA and disease associations analysis software. Hum Immunol. 2008;69(8):522–5.View ArticlePubMedGoogle Scholar
 Svejgaard A, Ryder LP. HLA and disease associations: detecting the strongest association. Tissue Antigens. 1994;43(1):18–27.View ArticlePubMedGoogle Scholar
 Duggal P, Gillanders EM, Holmes TN, BaileyWilson JE. Establishing an adjusted pvalue threshold to control the familywide type 1 error in genome wide association studies. BMC Genomics. 2008;9:516–6.View ArticlePubMedPubMed CentralGoogle Scholar
 Gao X, Becker LC, Becker DM, Starmer JD, Province MA. Avoiding the high Bonferroni penalty in genomewide association studies. Genet Epidemiol. 2010;34:100–5.PubMedPubMed CentralGoogle Scholar
 Johnson R, Nelson G, Troyer J, Lautenberger J, Kessing B, Winkler C, O’Brien S. Accounting for multiple comparisons in a genomewide association study (GWAS). BMC Genomics. 2010;11:724.View ArticlePubMedPubMed CentralGoogle Scholar
 Perneger TV. What’s wrong with Bonferroni adjustments. Brit Med J. 1998;316:1236–8.View ArticlePubMedPubMed CentralGoogle Scholar