ORIGINAL RESEARCH published: 21 August 2020 doi: 10.3389/fgene.2020.00951 High-Throughput Sequencing to Investigate Associations Between HLA Genes and Metamizole-Induced Agranulocytosis Anca Liliana Cismaru 1,2† , Livia Grimm 1† , Deborah Rudin 3,4 , Luisa Ibañez 5 , Evangelia Liakoni 6,7 , Nicolas Bonadies 8 , Reinhold Kreutz 9 , Pär Hallberg 10 , Edited by: Mia Wadelius 10 , EuDAC Collaborators, Manuel Haschke 6,7 , Carlo R. Largiadèr 1 and Ulrich M. Zanger, Ursula Amstutz 1* Dr. Margarete Fischer-Bosch Institut 1 Department of Clinical Chemistry, Inselspital, Bern University Hospital, University of Bern, Bern, Switzerland, 2 Graduate für Klinische Pharmakologie (IKP), School for Cellular and Biomedical Sciences, University of Bern, Bern, Switzerland, 3 Division of Clinical Pharmacology & Germany Toxicology, University Hospital Basel, Basel, Switzerland, 4 Department of Biomedicine, University of Basel, Basel, Reviewed by: Switzerland, 5 Clinical Pharmacology Service, Hospital Universitari Vall d’Hebron, Department of Pharmacology, Therapeutics Leonid Padyukov, and Toxicology, Fundació Institut Català de Farmacologia, Autonomous University of Barcelona, Barcelona, Spain, 6 Clinical Karolinska Institutet (KI), Sweden Pharmacology and Toxicology, Department of General Internal Medicine, Inselspital, Bern University Hospital, University Jonathan Marchini, of Bern, Bern, Switzerland, 7 Institute of Pharmacology, University of Bern, Bern, Switzerland, 8 Department of Hematology Regeneron Pharmaceuticals, Inc., and Central Hematology Laboratory, Inselspital, Bern University Hospital, University of Bern, Bern, Switzerland, 9 Charité - United States Universitätsmedizin Berlin, Corporate Member of Freie Universität Berlin, Humboldt-Universität zu Berlin, and Berlin Institute Loukas Moutsianas, of Health, Institut für Klinische Pharmakologie und Toxikologie, Berlin, Germany, 10 Department of Medical Sciences, Clinical Wellcome Sanger Institute (WT), Pharmacology and Science for Life Laboratory, Uppsala University, Uppsala, Sweden United Kingdom *Correspondence: Ursula Amstutz Background and Objective: Agranulocytosis is a rare and potentially life-threatening ursula.amstutz@insel.ch complication of metamizole (dipyrone) intake that is characterized by a loss of circulating † These authors share first authorship neutrophil granulocytes. While the mechanism underlying this adverse drug reaction is not well understood, involvement of the immune system has been suggested. In Specialty section: This article was submitted to addition, associations between genetic variants in the Human Leukocyte Antigen (HLA) Pharmacogenetics region and agranulocytosis induced by other drugs have been reported. The aim of and Pharmacogenomics, the present study was to assess whether genetic variants in classical HLA genes are a section of the journal Frontiers in Genetics associated with the susceptibility to metamizole-induced agranulocytosis (MIA) in a Received: 27 March 2020 European population by targeted resequencing of eight HLA genes. Accepted: 29 July 2020 Published: 21 August 2020 Design: A case-control cohort of Swiss patients with a history of neutropenia or Citation: agranulocytosis associated with metamizole exposure (n = 53), metamizole-tolerant Cismaru AL, Grimm L, Rudin D, (n = 39) and unexposed controls (n = 161) was recruited for this study. A high-throughput Ibañez L, Liakoni E, Bonadies N, Kreutz R, Hallberg P, Wadelius M, resequencing (HTS) and high-resolution typing method was used to sequence and EuDAC Collaborators, Haschke M, analyze eight HLA loci in a discovery subset of this cohort (n = 31 cases, n = 38 controls). Largiadèr CR and Amstutz U (2020) Identified candidate alleles were investigated in the full Swiss cohort as well as in two High-Throughput Sequencing to Investigate Associations Between independent cohorts from Germany and Spain using HLA imputation from genome- HLA Genes and Metamizole-Induced wide SNP array data. In addition, variant calling based on HTS data was performed in Agranulocytosis. Front. Genet. 11:951. the discovery subset for the class I genes HLA-A, -B, and -C using the HLA-specific doi: 10.3389/fgene.2020.00951 mapper hla-mapper. Frontiers in Genetics | www.frontiersin.org 1 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis Results: Eight candidate alleles (p < 0.05) were identified in the discovery subset, of which HLA-C∗ 04:01 was associated with MIA in the full Swiss cohort (p < 0.01) restricted to agranulocytosis (ANC < 0.5 × 109 /L) cases. However, no candidate allele showed a consistent association in the Swiss, German and Spanish cohorts. Analysis of individual sequence variants in class I genes produced consistent results with HLA typing but did not reveal additional small nucleotide variants associated with MIA. Conclusion: Our results do not support an HLA-restricted T cell-mediated immune mechanism for MIA. However, we established an efficient high-resolution (three- field) eight-locus HTS HLA resequencing method to interrogate the HLA region and demonstrated the feasibility of its application to pharmacogenetic studies. Keywords: drug-induced agranulocytosis, metamizole (dipyrone), HLA, pharmacogenetics, high-throughput sequencing, next generation sequencing, genetic association studies INTRODUCTION such as HLA-B∗ 27:05, HLA-B∗ 08:01, HLA-A∗ 31:01 and HLA- DQB1 alleles, respectively (Hallberg et al., 2016; Legge et al., Metamizole (dipyrone) is a non-opioid pro-drug with analgesic 2017; Wadelius et al., 2018). Furthermore, the frequencies of HLA and antipyretic properties that has been used in clinical practice alleles can vary between populations and consequently so could since the early 1930s (Hinz et al., 2016). Despite an overall the susceptibility to associated adverse drug reactions (Lucena good efficacy and low gastrointestinal toxicity compared to other et al., 2011). A recent commentary suggested the possibility medications (Arellano and Sacristán, 1990; Andrade et al., 1998; of an increased susceptibility to MIA in the British, Irish and Sánchez et al., 2002; Blaser et al., 2015), exposure to this drug Scandinavians compared to other populations, which could be has been associated with the risk of agranulocytosis, a rare mediated by specific HLA alleles with a higher frequency in and potentially fatal blood dyscrasia characterized by a decrease these populations (Shah, 2018). So far, larger studies investigating in circulating neutrophil granulocytes (Blaser et al., 2017). possible associations between genetic variation in the HLA Previous studies have reported widely varying risk estimates region and MIA are lacking and are sorely needed to improve for metamizole-induced agranulocytosis (MIA) such as 1:1500 our knowledge of the mechanism underlying the hematological prescriptions in Sweden and 1:1’000’000 inhabitants and year complications of pharmacotherapy with metamizole and to in Spain (Hedenmalm and Spigset, 2002; Ibáñez et al., 2005). attempt to minimize the risk of MIA in the near future. Currently, national regulations worldwide differ concerning The human leukocyte antigen (HLA) region stands out from metamizole use (Van der Klauw et al., 1999; Ibáñez et al., 2005; the rest of the human genome by virtue of its exceptional level Basak et al., 2010) with market withdrawal in countries such as of polymorphism and high gene density (Shiina et al., 2009; the United States, the United Kingdom, Denmark, and Sweden, Robinson et al., 2015). A wealth of associations between the HLA whereas it remains a frequently used prescription drug in other region and adverse drug reactions (ADRs) have been reported, countries like Germany and Switzerland, or even readily available yet the identification of underlying causal variants has been over-the-counter in countries like China, Mexico, and Brazil handicapped by the strong linkage disequilibrium (LD) in this (Chan and Chan, 1996; Hedenmalm and Spigset, 2002; Blaser region (Shiina et al., 2009; McCormack et al., 2011; Alfirevic et al., 2015; Stammschulte et al., 2015; Miljkovic et al., 2018). et al., 2012; Hosomichi et al., 2013). With the advent of high- The underlying mechanisms by which metamizole induces throughput sequencing (HTS), it is now possible to profile HLA agranulocytosis have not yet been fully elucidated. Similarly, alleles at unprecedented resolution and reduced ambiguity and risk factors that may render some patients more susceptible analyze full-length genes including intronic regions, potentially to this adverse drug reaction are poorly understood and there enabling deeper insights into the etiopathogenesis of HLA-ADR are no effective strategies to prevent the occurrence of MIA. associations (Kishore and Petrek, 2018). Previous studies have suggested that the mechanism of toxicity While several commercial HTS-based HLA typing kits are may involve an immunological component, including T-cell available, their cost may be prohibitive for population-based mediated elimination of neutrophils (Claas, 1996; Uetrecht, studies in a research context. Alternatively, non-commercial 1996; García-Martínez et al., 2003; Andres et al., 2009). Given primers to amplify full-length HLA genes have been described their role in cell-mediated innate and adaptive immunity, an (Shiina et al., 2012; Hosomichi et al., 2013; Ehrenberg et al., association of specific human leukocyte antigen (HLA alleles), 2017), as well as several freely available bioinformatics software as indicated by a small previous study (Vlahov et al., 1996), tools to achieve HLA typing from HTS reads (Warren et al., may provide evidence of such an immune-mediated mechanism. 2012; Bai et al., 2014; Nariai et al., 2015; Kawaguchi et al., 2017; Indeed, agranulocytosis induced by several other drugs such Lee and Kingsford, 2018). Previously published non-commercial as thiamazole (methimazole), clozapine and sulfasalazine has HLA HTS methods (Hosomichi et al., 2013; Ehrenberg et al., recently been found to be associated with specific HLA alleles 2017) used the tagmentation approach during sequencing library Frontiers in Genetics | www.frontiersin.org 2 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis preparation, which may lead to a biased pattern of coverage, such the absolute neutrophil count (ANC) but on the lowest number as an underrepresentation of GC-rich segments in introns 2 and 3 of potentially confounding factors such as comorbidities (i.e., of HLA-A (Lima et al., 2019). Potentially, this bias may be reduced ongoing infectious diseases such as HIV), chemotherapy, and with alternative fragmentation approaches. immunosuppressive therapy with cytotoxic drugs that may have Overall, while the literature is highly enthusiastic about HTS- contributed to agranulocytosis or neutropenia (Hallberg et al., based HLA typing (Carapito et al., 2016; Kishore and Petrek, 2016). Due to limits of sequencing library multiplexing, two 2018), less attention has been brought to novel challenges that of the originally selected cases could not be included for HTS go hand in hand with these new methods. The vast degree of analysis. In addition, the HTS discovery subset included all polymorphism of HLA genes has spawned new bioinformatics metamizole-tolerant individuals of European ancestry (n = 38). challenges associated with profiling HLA HTS data, such as This discovery subset was used to obtain high-resolution HLA identifying the best match of sequence reads to highly similar typing for a cost-effective discovery of potential candidate alleles reference sequences (e.g., different alleles of the same gene as well associated with MIA/MIN. as conserved regions between genes) and the lack of reference For follow-up of candidate alleles, HLA imputation data from sequences encompassing the full genomic sequence for many the full Swiss cohort, including the 14 cases not selected for alleles, which have recently been summarized (Klasberg et al., HTS analysis and 153 population controls was used, as well 2019). In this study, we assessed the feasibility of high-throughput as data from 286 independent subjects recruited in Germany sequencing of eight full-length classical HLA genes using non- (12 cases and 96 controls) and Spain (31 cases and 147 commercial primers and evaluated freely available bioinformatics controls) by investigators from the European Drug-induced tools to study classical HLA genes at the level of individual Agranulocytosis Consortium (EuDAC) as described previously alleles as well as individual small nucleotide variants (SNVs) to (Hallberg et al., 2016). investigate potential associations between variation in these genes and metamizole-induced agranulocytosis. DNA and Template Preparation Genomic DNA from whole blood or buffy coat was extracted and purified using the QIAamp DNA Blood Maxi and Midi MATERIALS AND METHODS Kits (Qiagen) according to the manufacturer’s instructions. DNA samples from subjects of the discovery cohort were Participants and Setting amplified at eight loci (HLA-A, -B, -C, -DRB1, -DQA1, Patients for this retrospective observational case-control study -DQB1, -DPA1 and -DPB1) by full gene-length long-range were recruited at two Swiss University Hospital (Basel and PCR (Polymerase Chain Reaction). Novel primers for HLA- the Inselspital Bern University Hospital). Patients diagnosed DQA1, -DPA1 (reverse primer), and -DPB1 (reverse primer) were with new-onset neutropenia or agranulocytosis while under designed in-house, while a combination of previously described metamizole therapy between 2005 and 2017, metamizole- locus-specific primers (Shiina et al., 2012; Hosomichi et al., tolerant control patients, as well as healthy population controls 2013; Ehrenberg et al., 2017) was used to amplify the other five were included. Collection of biological materials and clinical genes (Supplementary Table S1). Each locus was amplified in information was undertaken with a written informed consent 25 µL reaction volumes consisting of 50–250 ng of gDNA, 1x and following the principles of good clinical practice according PrimeSTAR GXL Buffer, 0.625 U PrimeSTAR GXL Polymerase to the Declaration of Helsinki. This study was approved by (Takara Bio Inc), 0.2 µM dNTPs and 0.08 mM of the respective the local ethics committees (“Ethikkommission Nordwest- und forward and reverse primer mix using the PCR cycling conditions Zentralschweiz” and “Kantonale Ethikkommission Bern,” BASEC described in Supplementary Table S2. To confirm amplification, protocol number 2015–00231). Clinical data was collected as 2 µL of each PCR product were analyzed on 0.8% agarose described previously (Rudin et al., 2019c). gels. The remaining product was purified by the QIAquick Cases were restricted to ≥18 years of age and development PCR Purification Kit (Qiagen) and quantified with a Qubit 4 of neutropenia (absolute neutrophil count < 1.5 × 109 /L) Fluorometric system (Life Technologies). within at least 1 day after the first intake or at the latest 2 weeks after discontinuing metamizole. Neutrophil counts prior to metamizole intake were not available for the majority Multiplex Library Preparation and of cases. Metamizole-tolerant control patients had received Sequencing of the Discovery Subset metamizole treatment with a daily dose of at least 500 mg for For library preparation with the QIAseq FX DNA Library Kit at least 28 consecutive days without any reported drug-related (Qiagen), all eight HLA amplicons were pooled in equimolar hematological complications. Population controls with no prior amounts (except for HLA-DRB1, where the total input amount exposure to metamizole were recruited via advertisement in was doubled due to lower coverage observed in preliminary regional print or online media. experiments) in two amplicon pools according to their sizes (see Out of the total of 45 cases having met the abovementioned Supplementary Table S1) resulting in two libraries per subject. inclusion criteria, a discovery subset of 31 metamizole-induced Amplicon pools were enzymatically fragmented (3 min for the agranulocytosis/neutropenia (MIA/MIN) patients was selected shorter pool “S” and 4 min for the longer pool “L”), then subjected for HTS analysis based on clinical factors supporting metamizole to size-selection clean-up (GeneRead Size Selection Kit, Qiagen), as the most likely causing agent. Selection was not based on uniquely dual-indexed at the subject level, purified with magnetic Frontiers in Genetics | www.frontiersin.org 3 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis beads (CleanNGS Kit, GC biotech) and the resulting libraries Hardy-Weinberg equilibrium (HWE) using Fisher’s exact test of “S” and “L” amplicons of all subjects were pooled separately implemented in Genepop on the web v.4.2 (Raymond and in equimolar amounts. For sequencing, one part of pooled “S” Rousset, 1995; Rousset, 2008). Depth of coverage was assessed libraries (3 nM) and three parts (volume) of pooled “L” libraries using the mean depth of coverage in exon 2 (or exon 3 in the (3 nM) were combined and the final library was diluted to an case of HLA-DRB2, which lacks exon 2) and refers to the sum of approximate concentration of 2 nM. the depth of coverage of the two alleles of the respective gene and The final library pool was diluted in HT-1/5% Tris-HCL, sample as reported by HLA-HD. pH8.5 with 5% PhiX spiked in and sequenced on a MiSeq instrument using the paired-end 500 cycle (2x 250bp paired- Genotyping Using Genome-Wide SNP end reads) MiSeq Reagent v2 Kit (both Illumina). Libraries for Arrays two samples were sequenced as part of the protocol optimization Genome-wide genotyping was performed using the Infinium using a MiSeq Nano v2 (2 × 250 bp paired-end reads) Kit. CoreExome-24 V.1.2 BeadChip (Illumina) in the 253 Swiss The PCR amplification and library preparation protocol used individuals recruited for the present study (53 cases, 39 tolerant, in this study was adapted based on preliminary experiments and 161 unexposed controls). Thirty-one agranulocytosis cases using anonymized samples previously HLA typed by validated from Spain and 147 Spanish controls had previously been SSP-PCR based methods (LinkSeq for HLA, Linkage Biosciences; genotyped with the Illumina HumanOmni 2.5M chip (Illumina), HLA-FluoGene ABC, Innotrain). Using this optimized protocol, while an additional 36 Spanish controls had been genotyped 100% concordant results with SSP-PCR typing were obtained with the Illumina HumanOmni1-Quad 1M chip. Cases (n = 12) for HLA-A, -B, -C, -DQA1, -DQB1, -DPA1, -DPB1 and 95% and controls (n = 96) from Germany had been genotyped concordant results in HLA-DRB1 in 10 investigated subjects (20 with the Illumina HumanOmniExpress 700K. Genotype calls alleles total; data not shown). were generated with the GenomeStudio Software (Illumina) and coordinates were based on the Genome Reference Consortium Quality Control and HLA Typing of the GRCh37 (hg19) build. Discovery Subset Standard quality control (Turner et al., 2011; Guo et al., 2014; Quality of the obtained raw sequence reads was assessed using Alonso et al., 2015) and data management was performed with FastQC v.0.11.71 and MultiQC v.1.5 (Ewels et al., 2016) and PLINK v.1.7 and v.1.9 (Purcell et al., 2007). Specifically, data was trimmed as described in more detail in the Supplementary filtered individually for each cohort (Swiss, Spanish, German) Appendix using the BiomedicalWorkbench v.4.1.1 (QIAgen). on the basis of single nucleotide polymorphisms (SNP) genotype Briefly, trimming and read filtering consisted of removal of call rates (>98%), sample genotype call rates (>98%), deviation Universal Illumina Adapters, removal of PCR primer sequences, from HWE in controls (p < 0.001) and minor allele frequency and removal of reads with length < 100 bp. For each individual, (MAF > 1%). paired reads from both sequencing libraries were combined for We performed multi-dimensional scaling (MDS) analysis subsequent analyses. including the combined data from the Swiss and EuDAC cohorts HLA typing of the trimmed and paired sequencing reads together with data from HapMap phase III samples to identify was performed with both the freely available tool HLA-HD nine and five ethnic outliers in the Swiss and the EuDAC v.1.2.0.1 (Kawaguchi et al., 2017) using read mapping with cohorts, respectively. Based on pairwise identity-by-descent Bowtie 2 (v.2.3.4.1) (Langmead and Salzberg, 2012) and the (IBD) estimation, we confirmed that three pairs of individuals IPD-IMGT/HLA reference dictionary v.3.34.0 (Robinson et al., in the EuDAC cohorts were closely related. One individual from 2015), as well as with the commercial software NGSengine each pair was randomly selected to be excluded from subsequent v.2.8.0.9796 (GenDx, IPD-IMGT/HLA database v.3.30.0). Details analyses. MDS was also applied to variants in the HLA region are provided in the Supplementary Appendix. to investigate potential global genetic differences between the Typing results from HLA-HD and NGSengine were compared cohorts in this region (see more detailed description in the at three-field resolution, with the exception of HLA-DRB1. If Supplementary Appendix). typing results were discordant, allele frequency information as reported in the Allele Frequency Net Database (González- Imputation of SNPs and HLA Types From Galarza et al., 2015), as well as reported DRB1-DQB1 haplotype Genome-Wide SNP Array Data frequencies (Maiers et al., 2007) were incorporated to infer the SNP imputation was performed by minimac3 on the Michigan more likely typing result. The plausibility of HLA-DRB1 typing Imputation Server (Das et al., 2016) using the Haplotype results was assessed based on the presence of related genes Reference Consortium (HRC 1.1 2016) reference panel and pseudogenes and information on known DRB haplotypes (McCarthy et al., 2015) on the basis of pre-imputation SNP-level (Kotsch and Blasczyk, 2000). Additional details on how the validated data2 from each cohort (see Supplementary Table S3). comparisons were performed are provided in the Supplementary HLA types of the five genes with at least one identified Methods. Genotype frequencies of the accepted three-field candidate allele (HLA-B, -C, -DQA1, -DQB1 and -DRB1) resolution HLA typing results were tested for deviations from were imputed using HIBAG (Zheng et al., 2014) with the 1 https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (last accessed 2 04.12.2019) https://www.well.ox.ac.uk/$\sim$wrayner/tools/ Frontiers in Genetics | www.frontiersin.org 4 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis published parameter estimates from the European reference Variants were converted to PLINK BED files using PLINK population based on hg19 (Zheng et al., 2014) in the three v.1.9 retaining only biallelic sites. Biallelic variants were filtered independent cohorts. HLA imputation of the Swiss cohort was on the basis of variant call rate (removing variants missing in > 1 performed using imputed SNPs with R2 ≥ 0.9 to increase the sample), minor allele frequency (removing variants with < 4 number of SNPs used in the imputation model. Unimputed minor variant allele count), and deviations from HWE in tolerant SNPs were used to impute HLA types of the German and controls (removing variants with HWE exact test p < 0.001). Spanish cohorts, since a higher proportion of SNPs of the imputation model were directly genotyped. HLA imputation accuracy was estimated using the samples from the discovery Statistical Analysis subset that had HLA typing results available from both Statistical analyses of the HTS typing results in the discovery HTS and HLA imputation assuming results from HTS as subset were performed using the R language for statistical the correct typing. computing v.3.5.3 (R Development Core Team,, 2019). Individual alleles at both three- and two-field resolution were tested for association with metamizole-induced Investigation of Individual Variants in neutropenia/agranulocytosis in a two-sided allelic Fisher’s exact Classical HLA Class I Genes From HTS test. Only alleles with at least four observations were included Data in the Discovery Subset in the association analyses. Given the exploratory nature of this Obtaining reliable variant calls in HLA genes has proved to discovery analysis, unadjusted p-values were reported and alleles be challenging with mapping bias impeding the accuracy of with p < 0.05 were selected as candidate alleles. genotype calls from short reads (Brandt et al., 2015). To Candidate alleles were subsequently tested in the full Swiss address this challenge, hla-mapper (Castelli et al., 2018) [v.2.3, and in the Spanish and German EuDAC cohorts to assess database v002.1, default parameters plus the flags –noclean and – replication of these candidate associations. Statistical analyses multiple_hits_MQ0; bwa v.0.7.17 (Li and Durbin, 2010) and were performed as described above but only at two-field SAMtools v.1.8 (Li et al., 2009)] was used to map the trimmed resolution available from HLA imputation using a two-sided reads to the respective HLA class I genes of the human reference allelic Fisher’s exact test and using the R language for statistical genome build hg19. Per-sample variants from the respective computing v.3.6.1 (R Development Core Team,, 2019). In the BAM files of HLA-A, -B, and -C were called using Genome Swiss cohort, candidate alleles were investigated in the full cohort, Analysis Toolkit (GATK) v.4.1.3.0 (Poplin et al., 2017) in GVCF as well as with agranulocytosis cases only (ANC < 0.5 × 109 /L; mode with default parameters considering chromosome 6 of n = 30). Association results were visualized using the R package hg19 as reference. Per-sample GVCFs were consolidated using forestplot (Gordon and Lumley, 2019). For replication analyses, GATK GenomicsDBImport and jointly genotyped using GATK the significance threshold was adjusted to the number of GenotypeGVCFs producing a set of jointly called SNPs and investigated candidate genes, i.e., p < 0.01 was considered Indels. Variants were annotated using bcftools (SAMtools) with statistically significant, corresponding to a Bonferroni correction the pattern ‘CHR:POS:REF:ALT.’ for 5 independent tests. TABLE 1 | Patient characteristics of the three independent cohorts. Cohort Swiss EUDAC EUDAC NGS discovery German Spanish Cases Controls Cases Controls Cases Controls N = 31(45) N = 38‡ (191) N = 12 N = 92 N = 29 N = 145 ANC < 500/µL [%] 25[87] NA 12 NA 29 NA Gender, male [%] 13[42] 17[45] 4[33] NA 6[21] NA Age, years [%] <25 7[23] 1[3] 2[16.6] NA 3[10.3] NA 25-44 11[35] 6[16] 6[50] NA 6[21] NA 45-64 10[32] 15[39] 2[16.6] NA 12[41.4] NA 65-74 3[10] 9[24] 2[16.6] NA 5[17] NA >74 – 7[18] – NA 3[10.3] NA BMI, median [range] 24 28 NA NA NA NA [19-47] [16-39] Latency time*a /treatment durationb , days 17 25 33.5 NA 11.5 NA [1-204] [1-5297] [4-9855] [1-235] The HTS discovery subset is shown for the Swiss cohort. a Cases, b controls. *Latency missing for 3 agranulocytosis cases and 1 neutropenia case. ‡ Metamizole-tolerant controls. NA, not available or not applicable. Frontiers in Genetics | www.frontiersin.org 5 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis Allelic associations of the filtered biallelic variants identified in the discovery subset as described in Section “Investigation of individual variants in classical HLA class I genes from HTS data in the discovery subset” were tested in PLINK v.1.9 using Fisher’s exact test. Finally, a discovery meta-analysis of all alleles was conducted for HLA-A, -B, -C, -DRB1, -DQA1, and -DQB1 including two-field and one-field association results from a two- sided allelic Fisher’s exact test for the three independent cohorts (71 agranulocytosis cases and 428 controls). The meta-analysis was performed using a sample-size-weighted fixed effect scheme with METAL (Willer et al., 2010). The analysis was based on p-values, taking both sample size and direction of effect relative to the effect allele into account, and testing for heterogeneity of effects between the cohorts. For this discovery meta-analysis, p < 1 × 10−3 was considered statistically significant. Given the linkage between HLA genes and particularly among class I and class II genes, this threshold was selected using the maximum FIGURE 1 | Mean depth of coverage in exon 2 as reported by HLA-HD of the number of alleles tested in any cohort for the most diverse class eight sequenced genes of interest. Note that two cases were sequenced as I and II genes, respectively, HLA-B (29 alleles in the Spanish part of preliminary experiments for protocol optimization and are thus not cohort) and HLA-DRB1 (20 alleles in the Spanish cohort) as an included in this figure. Additionally, samples that were sequenced but estimate for the number of independent tests. In each cohort, excluded for analysis are also excluded in this figure. only alleles with at least four observations were included in the initial meta-analysis. For two alleles with the lowest meta-analysis p-values, for which the initial meta-analysis did not include all samples (Figure 1). A considerably higher variation in depth of three cohorts due to this frequency cutoff, the meta-analysis was coverage between samples was observed in HLA-DRB1, with a repeated including data from all three cohorts. mean depth as low as 24x (Figure 1). Two cases sequenced as part of preliminary method optimization with a different pooling protocol were not included in these summary statistics. Mean RESULTS depth of coverage in exon 2 of HLA-DRB1 was dependent on the allele length with lower depth observed for longer alleles Participants (Supplementary Figure S1). In particular, alleles of the DR53 Of the 96 cases of metamizole-induced agranulocytosis recruited, haplotype, which is characterized by the presence of HLA-DRB4, ten were excluded as follows: eight cases from the Swiss cohort and to a lesser extent those of the DR52 haplotype, characterized (two as being ethnic outliers, two due to ongoing infectious by the presence of HLA-DRB3, showed lower depth of coverage diseases, four due to immunosuppressive therapy with cytotoxic (Supplementary Figure S1). In addition to the genes of interest, drugs) and two cases from the Spanish cohort (suggested ethnic additional HLA loci were amplified in some samples, specifically outliers). Out of the total 443 controls available, fifteen were the two HLA-A-related pseudogenes HLA-H and HLA-Y, the excluded as follows: nine from the Swiss cohort (one was HLA-DRB1-related genes HLA-DRB3, -DRB4, and -DRB5, and recruited twice as determined by IBD analysis, seven ethnic pseudogenes HLA-DRB2, -DRB6, and –DRB7 (Supplementary outliers and one with low sample call rate), two from the Figure S2). The presence of HLA-Y, of which only HLA-Y∗ 02:01 Spanish cohort (due to cryptic relatedness) and four from the was typed, appeared to be linked to specific HLA-A alleles German cohort (three ethnic outliers and one due to cryptic (Supplementary Table S4). Five out of eight genes of interest relatedness). As determined by the MDS analysis, the proportion displayed balanced allele ratios > 0.7 in every sample; whereas of participants of European ancestry among the total number HLA-DQB1 and HLA-DRB1 showed a wider range of observed recruited for the MIA-CH cohort was thus 96.5%, of which allele ratios (Supplementary Figure S3). 75.5, 74, and 74% in the case, tolerant control and unexposed In most instances, HLA typing results were accepted as control groups, respectively, reported “Swiss” as their ethnicity. reported by HLA-HD. Details of discordant typing results that Characteristics of patients and control subjects are summarized were investigated in more detail during quality control are shown in Tables 1, 2. The discovery subset included 25 out of 30 (full in Supplementary Table S5. In summary, 547 of 552 results Swiss cohort of MIA cases with agranulocytosis, ANC < 500/µL) were accepted from HLA-HD, of which 477 were identical at and 6 out of 15 cases of neutropenia (ANC < 1500/µL), while all three-field resolution to the result from NGSengine (total of cases of the German and Spanish cohorts had agranulocytosis. 483 comparisons, excluding HLA-DRB1), and three results were accepted from NGSengine. In addition, quality control identified Sequencing and HTS-Based Typing of allelic dropout of HLA-DRB1∗ 07 in two samples, which was the Discovery Subset confirmed in both instances using PCR SSP-based HLA typing The current HTS method resulted in a mean depth of coverage by the Transplantation Immunology Laboratory of the Center in exon 2 of >250x in seven out of eight target genes for all for Laboratory Medicine (Inselspital Bern University Hospital). Frontiers in Genetics | www.frontiersin.org 6 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis FIGURE 2 | Candidate alleles with evidence for association with MIA/MIN (p < 0.05) in the discovery subset at three-field resolution. Odds ratios (OR), 95% confidence intervals, and unadjusted p-values are shown. No deviation from HWE was detected for any of the eight allele, which was overrepresented in the control group, showed genes of interest. a stronger association with respect to MIA in the full cohort (Figure 4). HLA–C∗ 04:01 was strongly linked to HLA-B∗ 35:01 Association Analyses and Identification with 31 out of 32 carriers of the less common HLA-B∗ 35:01 allele of Candidate Alleles in the Discovery also carrying HLA-C∗ 04:01. When expanding the analyses to the EuDAC cohorts, Subset the general trend of the effect of these candidate alleles Association analysis of individual HLA alleles obtained with HTS typing in the discovery subset identified eight candidates for suggestive association with MIN/MIA in five different genes (Figure 2 and Supplementary Table S6). The same alleles were TABLE 2 | Patient characteristics of the full Swiss cohort. identified at three-field (Figure 2) and at two-field resolution Full Swiss cohort (Supplementary Figure S4). In total, 73 (three-field) and 75 (two-field) alleles were tested. Allele frequencies of all alleles Cases Tolerant Controls Unexposed Controls at three-field resolution in cases and tolerant controls of the N = 45 N = 38 N = 153 discovery subset are summarized in Supplementary Table S6. ANC < 500/µL 30 NA NA Gender, male [%] 23[51] 17[45] 72[47] Replication of Candidate Associations Age, years [%] Using HLA Imputation <25 10[22] 1[2.5] 27[17.5] Accuracy of HLA imputation from genome-wide SNP array data 25-44 14[31] 6[16] 73[48] in the discovery subset was >94% for four out of five imputed 45-64 12[27] 15[39.5] 53[34.5] genes (Table 3). Moreover, the identified candidate alleles were 65-74 7[15.5] 9[24] – imputed with high sensitivity and specificity in the discovery >74 2[4.5] 7[18] – subset (Table 4). BMI, median [range] 24[19-47] 28[16-39] 23[18-40] Overall, candidate alleles showed a similar trend in the full Latency time*a / 15[1-204] 25 [5197] NA Swiss cohort as in the discovery subset, albeit with a smaller effect treatment durationb , size (Figure 3). However, none of the candidate alleles showed days p < 0.01 in the full cohort, whereas a single allele (HLA–C∗ 04:01) a Cases, b metamizole-tolerantcontrols. *Latency missing for 3 agranulocytosis did so when restricting the analysis to cases with MIA only. This cases and 1 neutropenia case. NA, not applicable. Frontiers in Genetics | www.frontiersin.org 7 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis TABLE 3 | Accuracy of HLA imputation of the five genes where a candidate allele variation in the HLA region from genome-wide SNP was identified by HTS in the discovery cohort. array data was observed in Spanish versus German and Gene Imputation accuracy Swiss controls based on MDS analysis of HLA region SNPs (Supplementary Figure 5). HLA-B 96.4% (133/138) HLA-C 94.9% (131/138) HLA-DQA1 97.1% (134/138) Associations of Individual Sequence HLA-DQB1 98.6% (136/138) Variants in Full-Length Genes HLA-DRB1 88.4% (122/138) In total, 13 out of 725 tested individual variants obtained by HTS in the three full-length classical HLA class I Accuracies were determined by comparing the typing results of the discovery subset, i.e., the samples for which a typing result by HTS was available, to those genes were identified as candidates for an association with obtained by imputation. Number in brackets indicate the number of correctly MIA/MIN (p < 0.05) (Table 5). Over 50% of carriers of imputed alleles and the total number of alleles imputed for each gene. the alternate allele for 12 out of 13 candidate variants also carried the candidate allele HLA-C∗ 04:01 and 100% TABLE 4 | Sensitivity and specificity of HLA imputation for the eight candidate of carriers of eight of these variants also carried HLA- alleles identified by HTS in the discovery subset compared to HLA typing results of C∗ 04:01. The remaining candidate variant was carried together the HTS method. with the candidate allele HLA-C∗ 07:04 in all carriers of the Allele Sensitivity Specificity alternate allele. B*35:01 9/11 (82%) 127/127 (100%) C*04:01 13/13 (100%) 125/125 (100%) Discovery Meta-Analysis of Imputed HLA C*07:04 4/4 (100%) 134/134 (100%) Alleles Across Cohorts DQA1*01:02 24/24 (100%) 114/114 (100%) In the meta-analysis of all two-field and one-field resolution DQB1*05:01 18/18 (100%) 120/120 (100%) imputed HLA typing across all three cohorts, no alleles were DQB1*06:04 7/7 (100%) 131/131 (100%) significantly associated with MIA with a p-value of ≤1 × 10−3 DRB1*04:01 11/11 (100%) 127/127 (100%) (Supplementary Table S7). In the analysis of alleles at two-field DRB1*13:02 8/8 (100%) 130/130 (100%) resolution, the top allele with the smallest p-value was HLA- A∗ 68:02 (p = 3 × 10−3 ). This result was derived only from the Spanish cohort as this allele had less than four observations on MIA was not consistent in all three independent in the other two cohorts. Further exploration of the data for cohorts (Figure 5). The two most promising, linked this allele revealed that HLA-A∗ 68:02 was not observed in any candidate alleles HLA-B∗ 35:01 and –C∗ 04:01 showed a of the German or Swiss cases, and a meta-analysis including similarly protective trend in the Swiss and German cohorts the data from these two cohorts yielded a weaker association including only cases of agranulocytosis, but not in the result (p = 0.08, Supplementary Table S7). In the analysis Spanish cohort (Figure 4). No global difference of genetic of candidate alleles with one-field resolution, the allele with FIGURE 3 | Allele frequencies of the two linked candidate alleles, HLA-B∗ 35:01 and –C∗ 04:01, that showed a stronger association in the full cohort when restricting the analysis to agranulocytosis cases only, in all analyzed cohorts and subsets thereof. Frontiers in Genetics | www.frontiersin.org 8 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis FIGURE 4 | Associations of candidate alleles from the discovery subset at two-field resolution in the full Swiss cohort and in an analysis including only cases with MIA. HLA typing results obtained by HTS were used if available, else typing results obtained by HLA imputation were used. Alleles in strong linkage disequilibrium are highlighted in one contiguous gray block. P-values < 0.01 in the full cohort are highlighted with an asterisk. the smallest p-value was HLA-A∗ 31 (p = 4.5 × 10−3 ). This MIA/MIN. Candidate alleles were further investigated in an result did not consider the German cohort where this allele extended Swiss cohort as well as two independent European was not observed in any of the cases and only in two controls. replication cohorts. However, none of these alleles were Also for this allele, inclusion of all three cohorts resulted in consistently associated in all three cohorts and no significant an overall weaker association (p = 0.01, Supplementary Table associations of other alleles were identified in a discovery S7). When considering these results based on all data from all meta-analysis of all cohorts, thus challenging the hypothesis three cohorts, the candidate allele HLA-C∗ 04:01 was thus the of an HLA-restricted T cell-mediated mechanism behind top association signal (p = 0.0059) also in the meta-analysis metamizole-induced agranulocytosis. (Supplementary Table S7). Metamizole-Induced Agranulocytosis DISCUSSION and HLA Gene Variation The possibility of certain HLA alleles conferring an increased In this study, we investigated the association between susceptibility to MIA, which could put certain populations at a genetic variation in HLA genes and metamizole-induced greater risk, has previously been suggested (Vlahov et al., 1996; agranulocytosis/neutropenia (MIA/MIN). To our knowledge, Shah, 2018). However, in the current study, no major HLA this is the first pharmacogenetic association analysis of risk allele with a strongly increased frequency among patients full-length HLA genes using an HTS approach based with MIA or MIN was detected, thus making a T cell-mediated on non-commercial primers and utilizing freely available immune mechanism restricted by a specific HLA allele unlikely. software tools. This exploratory analysis in a Swiss discovery Interestingly, previous studies investigating agranulocytosis subset identified eight candidate alleles in five different related to other drugs found strong associations with HLA risk classical class I and II genes to be potentially associated with alleles (Yunis et al., 1995; Mallal et al., 2002, 2008; Illing et al., Frontiers in Genetics | www.frontiersin.org 9 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis FIGURE 5 | Associations of candidate alleles from the discovery subset at two-field resolution in the three independent cohorts, where cases are restricted to those meeting the criterion for agranulocytosis. HLA typing results obtained by HTS were used if available, else typing results obtained by HLA imputation were used. Alleles in strong linkage disequilibrium are highlighted in one contiguous gray block. All p-values < 0.01 are highlighted with an asterisk. TABLE 5 | Results from Fisher’s exact test for association of individual variants in the classical class I genes HLA-A, -B, and –C obtained by variant calling of the HTS data. CHR Position Ref Alt Cases Tolerants P OR Gene Exon/Intron Comment 6 31237779 C T 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Exon 5 Co-occured with C*04:01:01 6 31238027 C A 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Exon 4 6 31238267 C T 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Intron 3 6 31238538 G C 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Intron 3 6 31239082 G C 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Exon 3 6 31239501 G T 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Exon 2 6 31239607 G A 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Exon 2 6 31239742 G C 1 (2%) 12 (16%) 0.0063 0.087 HLA-C Intron 1 6 31239652 G A 3 (5%) 15 (20%) 0.011 0.21 HLA-C Intron 1 6 31239660 G T 3 (5%) 15 (20%) 0.011 0.21 HLA-C Intron 1 6 31237987 C T 2 (3%) 13 (17%) 0.012 0.16 HLA-C Exon 4 6 31239057 C T 5 (8%) 18 (24%) 0.020 0.28 HLA-C Exon 3 6 31238928 TCA T 4 (6%) 0 (0%) 0.039 NA HLA-C Exon 3 Co-occured with C*07:04:01 Only variants with p < 0.05 are shown and sorted according to p-value. Number of alleles and allele frequency (%) of the alternate allele (Alt) of cases and tolerant controls are indicated in the columns “Cases” and “Tolerants” respectively. Blocks of variants which are present in the same samples are highlighted in one block of the same shade of gray. If a shown variant co-occurred with a candidate HLA allele in the discovery subset, i.e., 100% of Alt allele carriers also carried the respective candidate allele, it is indicated in the column “Comment.” Genomic positions are based on hg19. Frontiers in Genetics | www.frontiersin.org 10 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis 2012; Chen et al., 2015; Hallberg et al., 2016; Saito et al., 2016; Daly et al., 2009) and others with somewhat weaker effect sizes He et al., 2017; Thao et al., 2018; Wadelius et al., 2018), some of still reported ORs in the range of 7–15 (McCormack et al., which have been replicated in independent cohorts. 2011; Wadelius et al., 2018). Given the sample size of the Our findings for MIA do not align with findings for other discovery subset, our study had sufficient power (80%) to detect agranulocytosis-inducing drugs, raising the question whether the associations of rarer (OR ≥ 8.5 for allele frequencies of ≥2%) mechanism of inducing agranulocytosis differs between causal and more common alleles with clinically relevant effect sizes [OR drugs. Indeed, our findings are consistent with recent functional of ≥5.0 for allele frequency of ≥5% and OR of ≥3.7 for allele studies indicating a direct toxic effect of the main metabolite of frequency of ≥10%; calculated using G∗ Power v3.1 (Faul et al., metamizole in the presence of hemin on granulocyte precursors 2007; Erdfelder et al., 2009)]. These minimum detectable ORs are (Rudin et al., 2019a,b). Moreover, a recent in vitro investigation thus lower than most previously reported associations of HLA of patient-derived peripheral blood cells did not detect any risk alleles mentioned above. activation of the adaptive immune system by either metamizole Another challenging aspect to be noted is the ascertainment or its main metabolite (Dina Manuel, 2019). To our knowledge, of metamizole as the causative agent for agranulocytosis and only one earlier study investigated HLA variation in the context the lack of data on neutrophil counts prior to metamizole of metamizole-induced agranulocytosis (Vlahov et al., 1996), administration. Multiple medications that are commonly reporting some associated alleles. However, this study included used concurrently with metamizole such as beta-lactam only five cases of MIA, and the reported alleles did not show the antibiotics, have also been associated with lowered neutrophil same effect in the cohorts analyzed here, which may be explained counts. However, restricting analyses to patients without by the very limited sample size of the previous investigation. such concomitant medications would result in very limited Taken together, these results thus do not support an HLA- sample sizes and an associated loss in statistical power. Finally, restricted T cell-mediated mechanism of MIA. while HLA alleles in the discovery subset were determined Of note, the candidate allele HLA-C∗ 04:01:01 with the at higher resolution using HTS, the typing resolution in the strongest evidence for association showed a putative protective replication cohorts based on HLA imputation was limited. effect, i.e., it was observed at very low frequency or was absent While concordance of HLA imputation at two-field resolution in MIA cases. Interestingly, a similar trend was observed for with HTS-based typing was generally high, particularly for this candidate allele and the linked allele HLA-B∗ 35:01 in the one of the candidate alleles (HLA-B∗ 35:01), a slightly reduced German, but only to a much weaker extent in the Spanish cohort. sensitivity was observed. Importantly, HLA imputation accuracy While no global difference was observed between Spanish, is limited by the resolution of the reference data available, which German and Swiss individuals by MDS analysis of SNPs in is currently primarily based on traditional HLA typing methods. the HLA region, we cannot exclude that individual alleles may differ in their frequency between these populations, resulting in differences in the observed associations with MIA as reported HTS-Based HLA Sequencing for for other idiosyncratic ADRs (Lucena et al., 2011). Nevertheless, Pharmacogenetic Association Studies given the exploratory nature of our analyses in the discovery With the present method, we obtained three-field resolution HLA subset, the identified candidates are expected to include false typing results, achieving higher resolution and less ambiguity positive results. Further validation in additional independent compared to traditional HLA typing methods based on core study populations would thus be required to assess whether HLA- exons alone or with HLA imputation. As a result, using HTS- C∗ 04:01, which also yielded the top association signal in the based direct analysis of the HLA region may have benefits for meta-analysis of all alleles in all cohorts, or any of the other future genetic association studies. In this context, our results putative candidates confer true associations. However, even if demonstrate the feasibility of HTS-based HLA analyses in such a protective effect of HLA-C∗ 04:01 were indeed present, pharmacogenetic studies using non-commercial methods and its underlying mechanisms would differ from a T cell mediated data analysis tools. mechanism mediated by a specific HLA allele. Importantly, an The development of HTS-based HLA typing approaches has association detected in classical HLA genes may also result from been met with great enthusiasm (Carapito et al., 2016; Kishore a mechanism mediated by other genes or regulatory variants and Petrek, 2018) but currently faces certain limitations. Similar in linkage disequilibrium (LD) with HLA alleles as it is often to other HTS-methods, our analysis included full-length genes, challenging to differentiate the effect of classical HLA genes and enabling, in principle, typing at four-field resolution. To date, other non-HLA genes or regulatory variants within this region however, this resolution was limited to three fields due to the (Shiina et al., 2009; McCormack et al., 2011; Alfirevic et al., 2012; limited characterization of intronic and UTR regions in the HLA Hosomichi et al., 2013). reference database (Kawaguchi et al., 2017). A more complete As MIA is a rare occurrence, identifying a large number of dictionary of full-length reference sequences for HLA alleles cases for genetic association studies is challenging. Nevertheless, is imperative to support full four-field resolution typing and due to the large effect sizes reported for HLA associations with maximize the potential of HTS-based HLA analysis. other idiosyncratic ADRs, our analyses in the Swiss discovery In contrast to another recent pharmacogenetic association subset were sufficiently powered to detect associations of similar study using HTS-based HLA typing (Nakatani et al., 2019), magnitude. Specifically, several such previous studies reported our analyses not only included HLA typing, but also assessed odds ratios (ORs) > 50 (Mallal et al., 2002; Hung et al., 2005; sequence variation at the level of SNVs in full-length class Frontiers in Genetics | www.frontiersin.org 11 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis I genes. Such analyses at the level of individual SNVs may HLA genes without amplifying closely related pseudogenes. complement traditional haplotype-level allele typing approaches Specifically, in addition to co-amplification of HLA-DRB1 related in association studies by enabling the detection of associations for genes, we also detected co-amplification of the HLA-A related SNVs contained on multiple HLA alleles, for which analyses at pseudogene HLA-Y (Kawaguchi et al., 2017). Of note, this the allele-level may be underpowered. Similarly, it may uncover pseudogene is not included in any of the human reference causal variants underlying known haplotype-level associations genome sequences (primary assembly or alternate haplotype with ADRs or diseases by providing effect sizes for individual thereof), making it difficult to account for during primer design. sequence changes. Given that such co-amplification of pseudogene sequences may We restricted our present analysis of individual HLA SNVs not be fully avoided for HLA genes, it is thus crucial that to class I genes as hla-mapper, the HLA-dedicated mapping their presence be considered at the data analysis stage. For tool, was only available for these genes. Due to the extensive example, HLA-HD detected reads originating from HLA-Y while sequence variation in HLA genes, variant calls in the HLA the version of NGSengine used in the presented analyses did region obtained from standard read mapping pipelines have been not. Given these numerous challenges associated with PCR reported to show a higher error rate compared to other parts amplification of full-length HLA genes, alternative PCR-free of the human genome with a reference allele bias (Brandt et al., approaches for isolating genomic target regions for subsequent 2015). The congruence of our variant calling results with the HTS, such as the use of CRISPR/Cas9 (Nachmanson et al., results of HLA typing demonstrates the feasibility of achieving 2018), may prove advantageous for future analysis of this sequence variation data at the level of individual SNVs with the genomic region. help of publicly available HLA-specific software. However, such software is currently scarce and not available for HLA class II genes, which display a more challenging pattern of variation, CONCLUSION including structural variants encompassing larger gene regions In this study, we did not find any consistent associations between and pseudogenes with highly similar sequences (Klasberg et al., HLA alleles and MIA across different cohorts by investigating 2019). Thus, for a more comprehensive assessment of HLA variation in HLA genes in detail. In particular, no strong variation at the level of SNVs, additional or updated HLA-specific enrichment of distinct HLA alleles in patients with MIA was mapping tools that include class II genes are needed. found. Our results, therefore, do not support an HLA-restricted The present method differs from other non-commercial HLA T cell-mediated immune mechanism of MIA. Nevertheless, our HTS methods (Hosomichi et al., 2013; Ehrenberg et al., 2017) study presents the feasibility of HTS-based analysis of HLA genes in that it relies on an enzymatic fragmentation of amplicons in the context of a pharmacogenetic association study, without and instead of a tagmentation approach in order to minimize the need to rely on costly commercial solutions, and highlights coverage bias (Lima et al., 2019). This alternative library the need for further development of dedicated bioinformatics preparation approach presented here may thus be applied to tools for HTS-based HLA data analysis. long-range PCR products generated with any available primers to amplify HLA genes. Furthermore, in contrast to the method used here, some commercially available HTS-based HLA typing kits EUDAC COLLABORATIVE AUTHORS do not amplify exon 1 and parts of intron 1 of most HLA class II genes, thus leading to ambiguous three-field typing results. This Niclas Eriksson (Uppsala Clinical Research Center and omission of exon 1 likely reflects the difficulty in amplifying full- Department of Medical Sciences, Uppsala University, Uppsala, length class II genes as single amplicons with even amplification Sweden); Mariam Molokhia (Department of Population Health efficacy across alleles. In particular, the genomic length of HLA- Sciences, King’s College London, London, United Kingdom); DRB1 varies considerably between alleles, primarily due to large Alfonso Carvajal (Centro de Estudios sobre la Seguridad de los structural variants in intron 1 (Von Salomé et al., 2007). Indeed, Medicamentos, Universidad de Valladolid, Valladolid, Spain); M. we observed a lower amplification efficiency of longer HLA-DRB1 Isabel Lucena (S Farmacologia Clinica, IBIMA, H Universitario alleles leading to more pronounced allelic imbalance and some Virgin de la Victoria, Universidad de Malaga, CIBERehd, allelic dropout. Amplification efficiency of some HLA-DRB1 Malaga, Spain); Javier Martin (Instituto de Parasitología y alleles may be further impaired by co-amplification of HLA-DRB4 Biomedicina Lopez-Neyra, CSIC, Granada, Spain); Esther or -DRB2. Here, we detected suspected allelic dropout for several Sancho Ponce (Servei d’Hematologia I Banc de Sang, Hospital samples using our quality control procedure. Nevertheless, given Genral de Catalunya, Sant Cugat del Vallès, Spain); Lourdes the low allele ratios for HLA-DQB1 observed in a few samples, Vendrell (Fundació Institut Català de Farmacologia, Hospital and it is possible that some additional allelic dropout was not Universitari Vall d’Hebron, Barcelona, Spain); Ramon Puig detected. This allelic imbalance could potentially be reduced by Treserra (Fundació Institut Català de Farmacologia, Barcelona, using higher PCR primer concentrations as previously described Spain); Jose Luis Caro (Banc de Sang i Teixits, Barcelona, (Hosomichi et al., 2013; Ehrenberg et al., 2017). Alternatively, for Spain); Eduard Palou (Banc de Sang i Teixits, Barcelona, future studies, it may be advantageous to amplify these class II Spain); María José Herrero (Banc de Sang i Teixits, Barcelona, genes in two separate amplicons (Shiina et al., 2012). Spain); Esmeralda de la Banda Ledrado (Hospital Universitario Finally, our study highlights additional inherent difficulties Bellvitge, Barcelona, Spain); Eva Montané (Hospital Universitari in identifying primers that amplify all alleles of full-length Germans Trias i Pujol, Universitat Autònoma de Barcelona); Frontiers in Genetics | www.frontiersin.org 12 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis Elisa Orna Montero (Institut Català d’Oncologia-Hospital FUNDING Germans Trias I Pujol, Badalona, Spain); José Tomás Navarro Ferrando (Catalan Institute of Oncology; Francesc This work was supported by funding from the Swiss National Salvador Rudilla, Banc de Sang i Teixits, Barcelona, Spain); Science Foundation (SNF grant no. 31003A_160206), the Consuelo Pedrós Cholvi (Hospital Universitario Bellvitge, Carlos III Spanish Health Institute (FIS10/02632), the European Barcelona, Spain). Regional Development Fund FEDER, the Swedish Research Council (Medicine 521-2011-2440, 521-2014-3370, and 2018- 03307), the Swedish Heart and Lung Foundation (20120557, DATA AVAILABILITY STATEMENT 20140291, and 20170711) and a grant from the Federal Institute for Drugs and Medical Devices (Bonn, Germany). The The datasets for this article are not publicly available due EUDRAGENE collaboration has received support from the EC to insufficient consent to share individual-level data publicly. 5th Framework program (QLRI-CR-2002-02757) and the Serious Aggregate data (allele frequencies for HTS based typing) are Adverse Events Consortium, SAEC, a collaboration for academia included as Supplementary Material to this article for all and industry. Mariam Molokhia is supported by the National investigated HLA genes. Requests to access individual-level Institute for Health Research Biomedical Research Center at datasets should be directed to the corresponding author. Guy’s and St Thomas’ National Health Service Foundation Trust and King’s College London. ETHICS STATEMENT ACKNOWLEDGMENTS The studies involving human participants were reviewed and We would like to cordially thank the laboratory technicians of approved by Ethikkommission Nordwest- und Zentralschweiz the University Institute of Clinical Chemistry (Inselspital Bern and Kantonale Ethikkommission Bern, Switzerland. The University Hospital), Gisela Andrey, Daniel Schärer, and Daniela patients/participants provided their written informed consent Sommer, for their support during the DNA extraction and library to participate in this study. preparation phases. We would like to thank the Transplantation Immunology Laboratory of the Center for Laboratory Medicine (Inselspital Bern University Hospital) for their support with AUTHOR CONTRIBUTIONS SSP-PCR based HLA typing. We further thank the Interfaculty Bioinformatics Unit of the University of Bern and its head, AC and LG contributed to the sample preparation, data Rémy Bruggmann, for providing access to their computational acquisition, and analysis of the data, as well as manuscript infrastructure to run our HTS data analyses. Last but not least, preparation. MH and UA were key senior contributors and made we would like to extend our thanks to Dr. Shuji Kawaguchi for substantial contributions to the study conception and design as his kind help with troubleshooting for the HLA-HD tool. Finally, well as interpretation of the data. DR, EL, LI, RK, and EuDAC we gratefully thank all study participants. collaborators recruited the study participants and collected the clinical data and biological material. CL provided input on the statistical analysis and data interpretation. MW and PH SUPPLEMENTARY MATERIAL contributed to the conception and design of the EuDAC study. All authors and EuDAC collaborators provided critical feedback The Supplementary Material for this article can be found to help shape the manuscript. All authors read and approved the online at: https://www.frontiersin.org/articles/10.3389/fgene. final manuscript. 2020.00951/full#supplementary-material REFERENCES Arellano, F., and Sacristán, J. A. (1990). Metamizole: reassessment of its therapeutic role. Eur. J. Clin. Pharmacol. 38, 617–619. doi: 10.1007/BF00278592 Alfirevic, A., Gonzalez-Galarza, F., Bell, C., Martinsson, K., Platt, V., Bretland, G., Bai, Y., Ni, M., Cooper, B., Wei, Y., and Fury, W. (2014). Inference of high et al. (2012). In silico analysis of HLA associations with drug-induced liver resolution HLA types using genome-wide RNA or DNA sequencing reads. BMC injury: use of a HLA-genotyped DNA archive from healthy volunteers. Genome Genomics. 15:325. doi: 10.1186/1471-2164-15-325 Med. 4:51. doi: 10.1186/gm350 Basak, G. W., Drozd-Sokołowska, J., and Wiktor-Jedrzejczak, W. (2010). Alonso, N., Lucas, G., and Hysi, P. (2015). Big data challenges in bone research: Update on the incidence of metamizole sodium-induced blood dyscrasias genome-wide association studies and next-generation sequencing. Bonekey Rep. in Poland. J. Int. Med. Res. 38, 1374–1380. doi: 10.1177/14732300100380 4:635. doi: 10.1038/bonekey.2015.2 0419 Andrade, S. E., Martinez, C., and Walker, A. M. (1998). Comparative safety Blaser, L., Hassna, H., Hofmann, S., Holbro, A., Haschke, M., Rätz Bravo, A., evaluation of non-narcotic analgesics. J. Clin. Epidemiol. 51, 1357–1365. doi: et al. (2017). Leucopenia associated with metamizole: a case-control study. Swiss 10.1016/S0895-4356(98)00076-6 Med. Wkly. 147, w14438. doi: 10.4414/smw.2017.14438 Andres, E., Dali-Youcef, N., Serraj, K., and Zimmer, J. (2009). Recognition and Blaser, L. S., Tramonti, A., Egger, P., Haschke, M., Krähenbühl, S., and Rätz management of drug-induced cytopenias: the example of idiosyncratic drug- Bravo, A. E. (2015). Hematological safety of metamizole: retrospective analysis induced thrombocytopenia. Expert Opin. Drug Saf. 8, 183–190. doi: 10.1517/ of WHO and Swiss spontaneous safety reports. Eur. J. Clin. Pharmacol. 71, 14740330902784162 209–217. doi: 10.1007/s00228-014-1781-z Frontiers in Genetics | www.frontiersin.org 13 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis Brandt, D. Y. C., Aguiar, V. R. C., Bitarello, B. D., Nunes, K., Goudet, J., and agranulocytosis in a Han population from northern China. Sci. Rep. 7:11950. Meyer, D. (2015). Mapping bias overestimates reference allele frequencies at doi: 10.1038/s41598-017-12350-2 the HLA genes in the 1000 genomes project phase I data. G3 5, 931–941. Hedenmalm, K., and Spigset, O. (2002). Agranulocytosis and other blood doi: 10.1534/g3.114.015784 dyscrasias associated with dipyrone (metamizole). Eur. J. Clin. Pharmacol. 58, Carapito, R., Radosavljevic, M., and Bahram, S. (2016). Next-generation 265–274. doi: 10.1007/s00228-002-0465-2 sequencing of the HLA locus: Methods and impacts on HLA typing, population Hinz, B., Cheremina, O., Bachmakov, J., Renner, B., Zolk, O., Fromm, M. F., et al. genetics and disease association studies. Hum. Immunol. 77, 1016–1023. doi: (2016). Dipyrone elicits substantial inhibition of peripheral cyclooxygenases in 10.1016/j.humimm.2016.04.002 humans new insights into the pharmacology of an old analgesic. FASEB J. 21, Castelli, E. C., Paz, M. A., Souza, A. S., Ramalho, J., and Mendes-Junior, C. T. 2343–2351. doi: 10.1096/fj.06-8061com (2018). Hla-mapper: an application to optimize the mapping of HLA sequences Hosomichi, K., Jinam, T. A., Mitsunaga, S., Nakaoka, H., and Inoue, I. (2013). produced by massively parallel sequencing procedures. Hum. Immunol. 79, Phase-defined complete sequencing of the HLA genes by next-generation 678–684. doi: 10.1016/j.humimm.2018.06.010 sequencing. BMC Genomics 14:355. doi: 10.1186/1471-2164-14-355 Chan, T. Y. K., and Chan, A. W. K. (1996). Aminopyrine-induced blood dyscrasias Hung, S.-I., Chung, W.-H., Liou, L.-B., Chu, C.-C., Lin, M., Huang, H.-P., et al. - still a problem in many parts of the world. Pharmacoepidemiol. Drug (2005). HLA-B∗ 5801 allele as a genetic marker for severe cutaneous adverse Saf. 5, 215–219. doi: 10.1002/(sici)1099-1557(199607)5:4<215::aid-pds208>3.0. reactions caused by allopurinol. Proc. Natl. Acad. Sci. U.S.A. 102, 4134–4139. co;2-5 doi: 10.1073/pnas.0409500102 Chen, P. L., Shih, S. R., Wang, P. W., Lin, Y. C., Chu, C. C., Lin, J. H., et al. Ibáñez, L., Vidal, X., Ballarín, E., and Laporte, J. R. (2005). Agranulocytosis (2015). Genetic determinants of antithyroid drug-induced agranulocytosis by associated with dipyrone (metamizol). Eur. J. Clin. Pharmacol. 60, 821–829. human leukocyte antigen genotyping and genome-wide association study. Nat. doi: 10.1007/s00228-004-0836-y Commun. 6:7633. doi: 10.1038/ncomms8633 Illing, P. T., Vivian, J. P., Dudek, N. L., Kostenko, L., Chen, Z., Bharadwaj, M., Claas, F. H. J. (1996). Immune mechanisms leading to drug-induced blood et al. (2012). Immune self-reactivity triggered by drug-modified HLA-peptide dyscrasias. Eur. J. Haematol. Suppl. 60, 64–68. doi: 10.1111/j.1600-0609.1996. repertoire. Nature 486, 554–558. doi: 10.1038/nature11147 tb01648.x Kawaguchi, S., Higasa, K., Shimizu, M., Yamada, R., and Matsuda, F. (2017). HLA- Daly, A. K., Donaldson, P. T., Bhatnagar, P., Shen, Y., Pe’er, I., Floratos, A., et al. HD: an accurate HLA typing algorithm for next-generation sequencing data. (2009). HLA-B∗ 5701 genotype is a major determinant of drug-induced liver Hum. Mutat. 38, 788–797. doi: 10.1002/humu.23230 injury due to flucloxacillin. Nat. Genet. 41, 816–819. doi: 10.1038/ng.379 Kishore, A., and Petrek, M. (2018). Next-generation sequencing based HLA typing: Das, S., Forer, L., Schönherr, S., Sidore, C., Locke, A. E., Kwong, A., et al. (2016). deciphering immunogenetic aspects of sarcoidosis. Front. Genet. 9:503. doi: Next-generation genotype imputation service and methods. Nat. Genet. 48, 10.3389/fgene.2018.00503 1284–1287. doi: 10.1038/ng.3656 Klasberg, S., Surendranath, V., Lange, V., and Schöfl, G. (2019). Bioinformatics Dina Manuel, D. N. (2019). Metamizole-Induced Neutropenia: Investigations of strategies, challenges, and opportunities for next generation sequencing- the adaptive immune system involvement. Ph.D. thesis, University of Bern, based HLA genotyping. Transfus. Med. Hemother. 46, 312–324. doi: 10.1159/ Switzerland. 000502487 Ehrenberg, P. K., Geretz, A., Sindhu, R. K., Vayntrub, T., Fernández Viña, Kotsch, K., and Blasczyk, R. (2000). The noncoding regions of HLA- M. A., Apps, R., et al. (2017). High-throughput next-generation sequencing to DRB uncover interlineage recombinations as a mechanism of HLA genotype six classical HLA loci from 96 donors in a single MiSeq run. HLA 90, diversification. J. Immunol. 165, 5664–5670. doi: 10.4049/jimmunol.165.10. 284–291. doi: 10.1111/tan.13133 5664 Erdfelder, E., FAul, F., Buchner, A., and Lang, A. G. (2009). Statistical power Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie analyses using G∗ Power 3.1: tests for correlation and regression analyses. Behav. 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923 Res. Methods. 41, 1149–1160. doi: 10.3758/BRM.41.4.1149 Lee, H., and Kingsford, C. (2018). Kourami: graph-guided assembly for novel Ewels, P., Magnusson, M., Lundin, S., and Käller, M. (2016). MultiQC: summarize human leukocyte antigen allele discovery. Genome Biol. 19:16. doi: 10.1186/ analysis results for multiple tools and samples in a single report. Bioinformatics s13059-018-1388-2 32, 3047–3048. doi: 10.1093/bioinformatics/btw354 Legge, S. E., Hamshere, M. L., Ripke, S., Pardinas, A. F., Goldstein, J. I., Rees, E., Faul, F., Erdfelder, E., Lang, A. G., and Buchner, A. (2007). G∗ Power 3: a flexible et al. (2017). Genome-wide common and rare variant analysis provides novel statistical power analysis program for the social, behavioral, and biomedical insights into clozapine-associated neutropenia. Mol. Psychiatry. 23, 162–163. sciences. Behav. Res. Methods 39, 175–191. doi: 10.3758/BF03193146 doi: 10.1038/mp.2016.97 García-Martínez, J. M., Vara, J. A. F., Lastres, P., Bernabéu, C., Betés, P. O., and Li, H., and Durbin, R. (2010). Fast and accurate long-read alignment with Burrows- Martín-Pérez, J. (2003). Effect of metamizole on promyelocytic and terminally Wheeler transform. Bioinformatics 26, 589–595. doi: 10.1093/bioinformatics/ differentiated granulocytic cells. Comparative Anal. with acetyl-salicylic acid btp698 diclofenac. Biochem. Pharmacol. 65, 209–217. doi: 10.1016/s0006-2952(02) Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). the 01511-3 sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. González-Galarza, F. F., Takeshita, L. Y. C., Santos, E. J. M., Kempson, F., Maia, doi: 10.1093/bioinformatics/btp352 M. H. T., Da Silva, A. L. S., et al. (2015). Allele frequency net 2015 update: Lima, T. H. A., Souza, A. S., Porto, I. O. P., Paz, M. A., Veiga-Castelli, L. C., Oliveira, new features for HLA epitopes, KIR and disease and HLA adverse drug reaction M. L. G., et al. (2019). HLA-A promoter, coding, and 3’UTR sequences in a associations. Nucleic Acids Res. 43, D784–D788. doi: 10.1093/nar/gku1166 Brazilian cohort, and their evolutionary aspects. HLA 93, 65–79. doi: 10.1111/ Gordon, M., and Lumley, T. (2019). forestplot: Advanced Forest Plot Using “grid” tan.13474 Graphics R package version 1.9. Available online at: https://CRAN.R-project.org/ Lucena, M. I., Molokhia, M., Shen, Y., Urban, T. J., Aithal, G. P., Andrade, R. J., package=forestplot et al. (2011). Susceptibility to amoxicillin-clavulanate-induced liver injury is Guo, Y., He, J., Zhao, S., Wu, H., Zhong, X., Sheng, Q., et al. (2014). Illumina influenced by multiple HLA class i and II alleles. Gastroenterology 141, 338–347. human exome genotyping array clustering and quality control. Nat. Protoc. 9, doi: 10.1053/j.gastro.2011.04.001 2643–2662. doi: 10.1038/nprot.2014.174 Maiers, M., Gragert, L., and Klitz, W. (2007). High-resolution HLA alleles and Hallberg, P., Eriksson, N., Ibañez, L., Bondon-Guitton, E., Kreutz, R., Carvajal, haplotypes in the United States population. Hum. Immunol. 68, 779–788. doi: A., et al. (2016). Genetic variants associated with antithyroid drug-induced 10.1016/j.humimm.2007.04.005 agranulocytosis: a genome-wide association study in a European population. Mallal, S., Nolan, D., Witt, C., Masel, G., Martin, A. M., Moore, C., et al. Lancet Diabetes Endocrinol. 4, 507–516. doi: 10.1016/S2213-8587(16)00 (2002). Association between presence of HLA-B∗ 5701, HLA-DR7, 113-3 and HLA-DQ3 and hypersensitivity to HIV-1 reverse-transcriptase He, Y., Zheng, J., Zhang, Q., Hou, P., Zhu, F., Yang, J., et al. (2017). Association inhibitor abacavir. Lancet 359, 727–732. doi: 10.1016/S0140-6736(02) of HLA-B and HLA-DRB1 polymorphisms with antithyroid drug-induced 07873-X Frontiers in Genetics | www.frontiersin.org 14 August 2020 | Volume 11 | Article 951 Cismaru et al. HLA and Pharmacogenetics of Metamizole-Agranulocytosis Mallal, S., Phillips, E., Carosi, G., Molina, J. M., Workman, C., Tomažiè, J., et al. subchronic treatment in rats. Dig. Dis. Sci. 47, 2791–2798. doi: 10.1023/A: (2008). HLA-B∗ 5701 screening for hypersensitivity to abacavir. N. Engl. J. Med. 1021077810548 2:RRN1203. doi: 10.1056/NEJMoa0706135 Shah, R. R. (2018). Metamizole (dipyrone)-induced agranulocytosis: does the risk McCarthy, S., Das, S., Kretzschmar, W., Delaneau, O., Wood, A., Teumer, A., vary according to ethnicity? J. Clin. Pharm. Ther. 44, 129–133. doi: 10.1111/jcpt. et al. (2015). A reference panel of 64,976 haplotypes for genotype imputation. 12768 bioRxiv[Preprint] doi: 10.1101/035170 Shiina, T., Hosomichi, K., Inoko, H., and Kulski, J. K. (2009). The HLA genomic McCormack, M., Alfirevic, A., Bourgeois, S., Farrell, J. J., Kasperavičiūtė, D., loci map: expression, interaction, diversity and disease. J. Hum. Genet. 54, Carrington, M., et al. (2011). HLA-A∗ 3101 and carbamazepine-induced 15–39. doi: 10.1038/jhg.2008.5 hypersensitivity reactions in Europeans. N. Engl. J. Med. 364, 1134–1143. doi: Shiina, T., Suzuki, S., Ozaki, Y., Taira, H., Kikkawa, E., Shigenari, A., et al. (2012). 10.1056/NEJMoa1013297 Super high resolution for single molecule-sequence-based typing of classical Miljkovic, M., Dragojevic-Simic, V., Rancic, N., Simic, R., Pekez-Pavlisko, T., HLA loci at the 8-digit level using next generation sequencers. Tissue Antigens Kovacevic, A., et al. (2018). Metamizole utilization and expenditure during 6- 80, 305–316. doi: 10.1111/j.1399-0039.2012.01941.x year period: Serbia vs. Croatia. Front. Public Heal. 6:213. doi: 10.3389/fpubh. Stammschulte, T., Ludwig, W. D., Mühlbauer, B., Bronder, E., and Gundert-Remy, 2018.00213 U. (2015). Metamizole (dipyrone)-associated agranulocytosis. an analysis of Nachmanson, D., Lian, S., Schmidt, E. K., Hipp, M. J., Baker, K. T., Zhang, Y., et al. German spontaneous reports 1990-2012. Eur. J. Clin. Pharmacol. 71, 1129– (2018). Targeted genome fragmentation with CRISPR/Cas9 enables fast and 1138. doi: 10.1007/s00228-015-1895-y efficient enrichment of small genomic regions and ultra-accurate sequencing Thao, M. P., Tuan, P. V. A., Linh, L. G. H., Van Hoang, L., Hen, P. H., Hoa, with low DNA input (CRISPR-DS). Genome Res. 28, 1589–1599. doi: 10.1101/ L. T., et al. (2018). Association of HLA-B – 38:02 with antithyroid drug-induced gr.235291.118 agranulocytosis in kinh vietnamese patients. Int. J. Endocrinol. 2018:7965346. Nakatani, K., Ueta, M., Khor, S. S., Hitomi, Y., Okudaira, Y., Masuya, A., et al. doi: 10.1155/2018/7965346 (2019). Identification of HLA-A∗ 02:06:01 as the primary disease susceptibility Turner, S., Armstrong, L. L., Bradford, Y., Carlsony, C. S., Crawford, D. C., HLA allele in cold medicine-related Stevens-Johnson syndrome with severe Crenshaw, A. T., et al. (2011). Quality control procedures for genome-wide ocular complications by high-resolution NGS-based HLA typing. Sci. Rep. association studies. Curr. Protoc. Hum. Genet. Chapter 1:Unit1.19. doi: 10.1002/ 9:16240. doi: 10.1038/s41598-019-52619-2 0471142905.hg0119s68 Nariai, N., Kojima, K., Saito, S., Mimori, T., Sato, Y., Kawai, Y., et al. Uetrecht, J. P. (1996). Reactive metabolites and agranulocytosis. Eur. J. Haematol. (2015). HLA-VBSeq: accurate HLA typing at full resolution from whole- Suppl. 57, 83–88. doi: 10.1111/j.1600-0609.1996.tb01651.x genome sequencing data. BMC Genomics. 16:S7. doi: 10.1186/1471-2164- Van der Klauw, M. M., Goudsmit, R., Halie, M. R., Van’t Veer, M. B., Herings, 16-S2-S7 R. M. C., Wilson, J. H. P., et al. (1999). A population-based case-cohort study of Poplin, R., Ruano-Rubio, V., DePristo, M. A., Fennell, T. J., Carneiro, M. O., drug-associated agranulocytosis. Arch. Intern. Med. 159, 369–374. doi: 10.1001/ Auwera, G. A., et al. (2017). Scaling accurate genetic variant discovery to tens archinte.159.4.369 of thousands of samples. bioRxiv[Preprint] doi: 10.1101/201178 Vlahov, V., Bacracheva, N., Tontcheva, D., Naumova, E., Mavrudieva, M., Ilieva, Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., P., et al. (1996). Genetic factors and risk of agranulocytosis from metamizol. et al. (2007). PLINK: a tool set for whole-genome association and population- Pharmacogenetics 6, 67–72. doi: 10.1097/00008571-199602000-00005 based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/51 Von Salomé, J., Gyllensten, U., and Bergström, T. F. (2007). Full-length 9795 sequence analysis of the HLA-DRB1 locus suggests a recent origin of alleles. R Development Core Team, (2019). R: A Language and Environment for Statistical Immunogenetics 59, 261–271. doi: 10.1007/s00251-007-0196-8 Computing. Vienna: R Development Core Team. Wadelius, M., Eriksson, N., Kreutz, R., Bondon-Guitton, E., Ibañez, L., Carvajal, Raymond, M., and Rousset, F. (1995). GENEPOP (Version 1.2): population A., et al. (2018). Sulfasalazine-induced agranulocytosis is associated with the genetics software for exact tests and ecumenicism. J. Hered. 86, 248–249. doi: human leukocyte antigen locus. Clin. Pharmacol. Ther. 103, 843–853. doi: 10. 10.1093/oxfordjournals.jhered.a111573 1002/cpt.805 Robinson, J., Halliwell, J. A., Hayhurst, J. D., Flicek, P., Parham, P., Warren, R. L., Choe, G., Freeman, D. J., Castellarin, M., Munro, S., Moore, R., et al. and Marsh, S. G. E. (2015). The IPD and IMGT/HLA database: allele (2012). Derivation of HLA types from shotgun sequence datasets. Genome Med. variant databases. Nucleic Acids Res. 43, D423–D431. doi: 10.1093/nar/gku 4:95. doi: 10.1186/gm396 1161 Willer, C. J., Li, Y., and Abecasis, G. R. (2010). METAL: fast and efficient meta- Rousset, F. (2008). GENEPOP’007: a complete re-implementation of the analysis of genomewide association scans. Bioinformatics 26, 2190–2191. doi: GENEPOP software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106. 10.1093/bioinformatics/btq340 doi: 10.1111/j.1471-8286.2007.01931.x Yunis, J. J., Corzo, D., Salazar, M., Lieberman, J. A., Howard, A., and Yunis, Rudin, D., Lanzilotto, A., Bachmann, F., Housecroft, C. E., Constable, E. C., Drewe, E. J. (1995). HLA associations in clozapine-induced agranulocytosis. Blood 86, J., et al. (2019a). Non-immunological toxicological mechanisms of metamizole- 1177–1183. doi: 10.1182/blood.v86.3.1177.bloodjournal8631177 associated neutropenia in HL60 cells. Biochem. Pharmacol. 163, 345–356. doi: Zheng, X., Shen, J., Cox, C., Wakefield, J. C., Ehm, M. G., Nelson, M. R., 10.1016/J.BCP.2019.01.011 et al. (2014). HIBAG—HLA genotype imputation with attribute bagging. Rudin, D., Roos, N. J., Duthaler, U., and Krähenbühl, S. (2019b). Pharmacogenomics J. 14, 192–200. doi: 10.1038/tpj.2013.18 Toxicity of metamizole on differentiating HL60 cells and human neutrophil granulocytes. Toxicology 426:152254. doi: 10.1016/j.tox.2019.15 Conflict of Interest: The authors declare that the research was conducted in the 2254 absence of any commercial or financial relationships that could be construed as a Rudin, D., Spoendlin, J., Cismaru, A. L., Liakoni, E., Bonadies, N., Amstutz, U., potential conflict of interest. et al. (2019c). Metamizole-associated neutropenia: comparison of patients with neutropenia and metamizole-tolerant patients. Eur. J. Intern. Med. 68, 36–43. Copyright © 2020 Cismaru, Grimm, Rudin, Ibañez, Liakoni, Bonadies, Kreutz, doi: 10.1016/j.ejim.2019.07.029 Hallberg, Wadelius, EuDAC Collaborators, Haschke, Largiadèr and Amstutz. This Saito, T., Ikeda, M., Mushiroda, T., Ozeki, T., Kondo, K., Shimasaki, is an open-access article distributed under the terms of the Creative Commons A., et al. (2016). Pharmacogenomic study of clozapine-induced Attribution License (CC BY). The use, distribution or reproduction in other forums agranulocytosis/granulocytopenia in a Japanese population. Biol. Psychiatry. is permitted, provided the original author(s) and the copyright owner(s) are credited 80, 636–642. doi: 10.1016/j.biopsych.2015.12.006 and that the original publication in this journal is cited, in accordance with accepted Sánchez, S., De la Lastra, C. A., Ortiz, P., Motilva, V., and Martín, M. J. (2002). academic practice. No use, distribution or reproduction is permitted which does not Gastrointestinal tolerability of metamizol, acetaminophen, and diclofenac in comply with these terms. Frontiers in Genetics | www.frontiersin.org 15 August 2020 | Volume 11 | Article 951