The HLA (Human Leukocyte Antigens) genes are well-documented targets of balancing selection, and variation at these loci is associated with many disease phenotypes. Variation in expression levels also influences disease susceptibility and resistance, but little information exists about the regulation and population-level patterns of expression. This results from the difficulty in mapping short reads originated from these highly polymorphic loci, and in accounting for the existence of several paralogues. We developed a computational pipeline to accurately estimate expression for HLA genes based on RNA-seq, improving both locus-level and allele-level estimates. First, reads are aligned to all known HLA sequences in order to infer HLA genotypes, then quantification of expression is carried out using a personalized index. We use simulations to show that expression estimates obtained in this way are not biased due to divergence from the reference genome. We applied our pipeline to the GEUVADIS dataset, and compared the quantifications to those obtained with reference transcriptome. Although the personalized pipeline recovers more reads, we found that using the reference transcriptome produces estimates similar to the personalized pipeline (r ≥ 0. 87) with the exception of HLA-DQA1. We describe the impact of the HLA-personalized approach on downstream analyses for nine classical HLA loci (HLA-A, HLA-C, HLA-B, HLA-DRA, HLA-DRB1, HLA-DQA1, HLA-DQB1, HLA-DPA1, HLA-DPB1). Although the influence of the HLA-personalized approach is modest for eQTL mapping, the p-values and the causality of the eQTLs obtained are better than when the reference transcriptome is used. We investigate how the eQTLs we identified explain variation in expression among lineages of HLA alleles. Finally, we discuss possible causes underlying differences between expression estimates obtained using RNA-seq, antibody-based approaches and qPCR.
The HLA region is the most polymorphic in the genome, and also shows the greatest number of disease associations, which has made it very well characterized at the genomic, population and functional levels [1,2]. Decades of research have also shown that the HLA genes are targets of natural selection, likely a consequence of their role in responding to pathogens [2,3]. This combination of evolutionary and biomedical interest has resulted in an extensive catalogue of HLA variation in human populations, with the frequency of HLA alleles defined for various populations [4–7]. The most intensely studied genes in the region are the classical HLA loci. These include Class I genes (HLA-A, HLA-C, and HLA-B), whose proteins present endogenous antigens to CD8+ T cells, and also regulate innate immune responses by interacting with killer cell immunoglobulin-like receptors (KIR) expressed on natural killer (NK) cells, and Class II genes (HLA-DRA, HLA-DRB1, HLA-DQA1, HLA-DQB1, HLA-DPA1, and HLA-DPB1) whose proteins present exogenous antigens to CD4+ T cells [1,8]. Variation within or near HLA genes has been convincingly linked to resistance and susceptibility to both autoimmune and infectious diseases [2,9, 10]. Although in many cases the mechanistic basis of the associations remains poorly understood, there has been an effort to identify whether associations can be linked to features such as variation at the level of specific amino-acids, HLA alleles, HLA haplotypes, non-coding variants near HLA genes, or HLA expression levels (reviewed in [2,10,11]). In some instances, the association of one feature results from the fact that it tags another feature. Several studies have shown that HLA alleles may have a protective effect because they mark overall gene expression [12–15]. Other studies showed that a combination of two features may be important. For example, mortality of transplant recipients was associated with increased expression of HLA-C alleles which harbor specific amino-acids [16]. An understanding of how HLA expression varies among individuals, and the identification of genetic variants involved in the regulation of expression, will play a central role in understanding the contribution of HLA genes to normal and disease phenotypes. However, until recently, little information existed about the regulatory variation and population-level expression patterns of HLA genes, a result of the difficulty in quantifying expression for genes which show an unusually high polymorphism and are members of a multi-gene family [9,17]. Efforts have been made to develop antibody-based methods to quantify HLA protein on the cell surface [12–14,18–20], or hybridization-based approaches to quantify mRNA, such as qPCR [13,21,22] and microarray methods [23]. However, the design of PCR primers, microarray probes, or antibodies that span the diversity of possible variants represents a technically challenging and labor-intensive undertaking. In addition, qPCR technologies are not appropriate for comparison of expression levels among different loci, an important concern when seeking to understand how expression of HLA genes responds to environmental challenges. The results obtained to date using qPCR [13,21,22], antibody-based approaches [12–14,18–20] and customized arrays [23] have contributed to the understanding of HLA expression and its underlying genetic regulation. However, they do not take advantage of the large amount of RNA-seq data generated by studies of whole transcriptomes in large samples [24–26], often involving populations from various regions of the world, which represents an attractive resource to investigate HLA expression. Although such whole-transcriptome RNA-seq studies do provide expression estimates for HLA genes, they bring new challenges. RNA-seq pipelines may provide biased expression estimates for two reasons: 1) many short reads originating from genes with extreme polymorphism fail to map to the reference genome, due to high degree of variation (which results in a large number of mismatches between the reference genome and that of most individuals being analyzed), and 2) the presence of paralogues makes it difficult to map a read uniquely to a specific gene, leading to the exclusion of many reads. This raises concerns about the reliability of RNA-seq approaches to quantify HLA expression, given that these loci represent both the extreme of polymorphism in the human genome and are part of a multi-gene family [27–29]. A strategy to overcome these challenges is the mapping of reads to an HLA-personalized reference (an index containing sequences of the individual being quantified), rather than to a single reference genome. For example, seq2HLA is a tool developed by Boegel et al. [30] to provide in-silico HLA types and expression estimates, and later applied to demonstrate that different tumor types are associated with different HLA expression levels [31], and also to provide a large catalog of HLA expression in 56 human tissues and cell types [32]. AltHapAlignR [33] is another software which infers the HLA references which are the closest to the individual’s HLA haplotypes, and maps reads to them. The authors reanalyzed a large RNA-seq dataset [24], and provided comparisons with conventional read mapping, showing an improvement in accuracy with the HLA-tailored pipeline. In this article we present HLApers, a personalized pipeline which we have developed to reliably quantify HLA expression from RNA-seq data. We compare HLApers to conventional pipelines, and discuss for the first time the impact of accurate estimation of HLA expression on downstream analyses such as eQTL (expression Quantitative Trait Loci) mapping and allele-specific expression. We show that it is possible to adapt different computationally efficient methods to work under the personalized reference framework, providing reliable quantification of HLA expression from RNA-seq data. We find that implementations with either a conventional read mapper [34] or a pseudoaligner [35] show similar expression estimates. We use simulations to assess accuracy, showing that the HLA-personalized pipeline is more accurate than conventional mapping, and apply the tool to reanalyze RNA-seq data from the GEUVADIS Consortium [24], which made available whole-transcriptome RNA-seq data for Lymphoblastoid Cell Lines (LCLs) from 462 European and African individuals. We evaluate the impact of more accurate expression estimates obtained with HLApers on downstream analyses by carrying out a detailed survey for allele-specific expression and eQTL mapping at the classical HLA loci. Surprisingly, we find that conventional RNA-seq pipelines provide gene-level expression estimates and identify eQTLs which are highly correlated with those obtained under the HLA-personalized approach. However, we identify gains of using a pipeline tailored for HLA expression: higher accuracy, expression estimates at the HLA allele-level, and eQTLs with higher probabilities of being causal.
We developed the HLApers pipeline (for HLA expression with personalized genotypes) to measure HLA expression from whole-transcriptome RNA-seq data. The pipeline can use either (1) a suffix array-based read mapper (STAR [34]) followed by quantification with Salmon [36] (henceforth called STAR-Salmon), or (2) a pseudoaligner with built-in quantification protocol (kallisto [35]). The key feature of our implementation is the use of an index supplemented with a set of sequences covering the breadth of known HLA sequences (see Materials and methods: Index supplemented with the HLA diversity). We implemented a two-step quantification approach, where (1) we align reads to reference sequences corresponding to all known HLA alleles, and identify those which maximize the read counts at each locus to infer the genotype which is present (in-silico genotyping), and (2) we use this inferred HLA genotype to create a personalized index which we use to quantify expression (Fig 1). A problem for the analysis of HLA is that of multimaps: reads originating from a particular allele which map to other alleles of the same locus, or to other loci. To deal with multimaps different approaches have been used. One possibility is to discard them (as in [33]), and another is to split them evenly among the compatible references (as in [30]). In HLApers, we use maximum likelihood estimates of expression obtained by an expectation-maximization (EM) algorithm (implemented within Salmon [36] and kallisto [35]), which infers the quantities of each HLA allele that maximize the probability of observing the set of sequenced reads. Because the in-silico typing is an important step for accurate expression estimates, we assessed the concordance between our RNA-seq based HLA typing and the HLA allele calls experimentally determined for 5 HLA loci using Sanger sequencing [6]. The concordance was higher than 97% for all of the HLA genes compared (S1 Table). This is consistent with previous results showing that RNA-seq provides reliable HLA alleles calls [30,37–40]. Our HLApers pipeline integrates widely used quantification tools such as Salmon and kallisto, and provides both allele and locus level estimates of HLA expression. The pipeline is available at https: //github. com/genevol-usp/HLApers. We next investigated how using the HLA-personalized index affects the quantifications of HLA expression. To this end we simulated an RNA-seq experiment using the Polyester package [41], with a dataset of 50 individuals, with the read length and counts matching those of the observed data (see Materials and methods: Simulation). We analyzed the simulated dataset with three different methodologies: (1) the two-step approach in HLApers, where we first inferred the personalized HLA genotype and then aligned reads to it; (2) alignment to the reference transcriptome (Gencode release 25; primary assembly); and (3) alignment to the reference genome (GRCh38). In approaches (1) and (2) the alignment is followed by expression quantification using Maximum Likelihood (ML), which provides a statistical framework for dealing with multimap reads, whereas in approach (3) quantification is performed using only uniquely mapped reads. For each HLA locus and methodology, we assessed the proportion of simulated reads which successfully aligned. Because previous studies for genome sequencing identified a correlation between mapping success and the number of mismatches between the HLA allele an individual carries and the reference genome [28], we analyzed how alignment success behaves as a function of the number of mismatches between each HLA allele and the reference genome (Fig 2). The use of an HLA-personalized index results in the largest proportion of successfully aligned reads, no matter how different the allele carried by the individual is from the allele in the reference genome. This is expected, since the personalized HLA component guarantees that a sequence close or identical to that originating the read will be present. When the alignment was performed using the reference transcriptome, there was a marked reduction in the proportion of successfully aligned reads for HLA-DRB1, HLA-DQA1, HLA-DQB1, driven by decreased alignment success for alleles with a greater proportion of mismatches with respect to the reference genome. When using uniquely mapped reads there was a massive read loss for HLA-A, HLA-B and HLA-DPB1, regardless of the divergence to the reference genome, as well as a lower proportion of successfully aligned reads across other loci. This shows that both discarding multipmaps, as well as not including a personalized index, have a negative impact on mapping success. Finally, for the least polymorphic HLA loci, mapping should not be sensitive to the specific reference used. This is precisely what we find, with all pipelines performing similarly for HLA-DRA and HLA-DPA1. Having demonstrated that including an individual’s HLA alleles in the index improves the success of read alignment in the simulated data (Fig 2), we set out to address two questions with real data by applying HLApers to the GEUVADIS dataset [24]. First, we examined how expression varies among HLA loci, when the personalized index is used (Fig 3). Secondly, we compared expression estimates with and without the use of the personalized index, so as to evaluate the impact of its usage on real data (Fig 4). By summing the estimates for the 2 alleles at each HLA locus, we obtain gene-level expression estimates (Fig 3). We observe that HLA-B is the highest expressed gene overall. Among the Class I genes, HLA-B is followed by HLA-A with similar levels, and by HLA-C which has about 50% of the expression levels of HLA-B. For Class II genes, HLA-DRA is the most highly expressed. Although we observe a general concordance with the original GEUVADIS quantifications [24], there are some notable differences: in the original quantifications, HLA-B is twice as expressed as HLA-A, and HLA-DPA1 is more expressed than HLA-DRB1 (S1 Fig). We found that the correlation between results using the reference transcriptome or the personalized index was greater than 0. 87 for every locus except for HLA-DQA1 (Fig 4). The loci with the lowest correlations between indices (HLA-DRB1, HLA-DQA1 and HLA-DQB1), are also those with the greatest read loss when divergence from the reference allele is high in the simulation (Fig 2). We then investigated if the quantification tool used in HLApers (STAR-Salmon or kallisto) influences expression estimates. Correlations between STAR-Salmon and kallisto approaches were on average r > 0. 99 for read counts, dropping to r = 0. 8 for TPM estimates for Class I genes, likely due to different bias correction models (S2 Fig). These results show that the key features influencing alignment success are the use of a personalized index and the statistical treatment of multimaps (as opposed to discarding them), with the specific alignment tool being less influential. Finally, we evaluated the reproducibility of the HLA expression estimates obtained with HLApers. To this end, we analyzed replicates for 97 European individuals in GEUVADIS. We observed an average correlation of expression estimates between replicates of 0. 92 over the 9 classical HLA loci. This shows that HLA quantifications from RNA-seq with HLApers are highly reproducible (S3 Fig). The key role played by HLA loci in the immune response and their strong and abundant associations with infectious and autoimmune diseases have motivated studies to uncover their regulatory architecture [15,20,21,23,42–45]. Here we use accurate expression estimates obtained with HLApers, together with genotype data from the 1000 Genomes Project [46], to identify SNPs which are associated with variation in expression levels (eQTLs). Because multiple SNPs can affect expression, it is interesting to identify independent contributions made by distinct SNPs. This has previously been done by using a best eQTL (i. e. , the one with the most extreme p-value) as a covariate in subsequent searches for an additional variant. Here we use a related approach implemented in QTLtools, a collection of tools to perform eQTL analysis [47]. The module QTLtools cis allows for the identification of groups (or “ranks”) of SNPs associated with independent signals. For each rank, we identified the site with the most extreme association (Fig 5). Transcriptome studies can quantify expression for various biological features: individual SNPs, exons, isoforms, genes. In the case of HLA loci, a natural unit of interest is the HLA allele. HLApers provides expression estimates for individual alleles, since it is the allelic sequences which are included in the index. The immunogenetics literature has shown that many HLA alleles are associated with specific phenotypes of evolutionary and medical importance (reviewed in [2,10,11]). Gauging information about the expression levels of alleles can therefore provide an additional layer of information. In order to explore allele-level estimates, we first grouped alleles in “lineages”, which comprise groups of alleles which are evolutionarily and functionally related, since the large number of alleles would make sample sizes per allele too sparse. Although the expression of individual allelic lineages is highly variable among individuals, there is an overall significant difference in expression among lineages (Welch’s ANOVA p-values ranging from 3. 7 × 10−8 for HLA-DPB1 to 6 × 10−51 for HLA-DQA1). Enhanced coordination of gene expression has been proposed as an advantage for gene clustering, as seen at the HLA region [1,2]. We found a high correlation of expression both within the group of Class I and Class II genes, and lower levels between Class I and Class II genes, which are more than 1Mb apart (Fig 8A) (but see [22] for a result of no co-expression among Class I genes). A possible cause for co-expression of genes which are physically close to each other is that regulatory activity is structured in domains (CRDs, for Cis Regulatory Domains) [52]. Such domains comprise contiguous regions along a chromosome, and their existence predicts co-expression along haplotypes for genes associated with the same CRD. In fact, previous studies did find evidence that expression of some genes in the HLA region was a feature associated with haplotype membership [23,62]. We used our inferences of HLA haplotype structure to investigate if there is an haplotypic effect on coordination of expression among nearby HLA genes. Specifically, we tested the hypothesis that co-expression is stronger between alleles located within the same haplotype than between those on different haplotypes. We did not find a consistently higher correlation of expression for alleles on the same haplotype (correlation within haplotypes being higher in only 7 out of 18 locus pairs surveyed within Class I or Class II, Fig 8A). This result suggests that correlation of expression among HLA loci in LCLs is a result of factors acting at the gene level, and is not driven predominantly by properties of the haplotype. For example, we observe correlation of expression between Class II genes and CIITA, the Class II Major Histocompatibility Complex Transactivator. This correlation is not driven by proximity (CIITA is on chromosome 16), but rather by a trans regulatory mechanism (Fig 8B).
The contribution of HLA variation to normal and disease phenotypes goes beyond peptide specificity, and includes other factors which influence the strength of immune responses, such as HLA expression levels. However, we are still only starting to understand the regulation of expression of these genes, a consequence of their extreme polymorphism which imposes challenges to the methodologies available to measure mRNA and protein levels. There has been an increasing effort to use RNA-seq data for in-silico HLA typing ([63] and other methods reviewed in [64]) and to estimate expression levels for HLA genes [30,33]. In this paper, we present HLApers, an HLA-personalized pipeline that quantifies HLA expression based on RNA-seq. We used these expression estimates to identify eQTLs for HLA genes, to estimate the degree of allelic imbalance in expression, and to investigate how the alleles at eQTLs explain variation in expression among HLA alleles. When new technologies/methods are developed, it is a good practice to evaluate if the results are in agreement with the current knowledge, which in the case of HLA was mostly derived from qPCR or antibody-based studies. However, in our survey of this literature, we found it difficult to compare our newly obtained expression estimates based on RNA-seq with those from previous studies. There are several reasons underlying this difficulty, which we now discuss. First, expression estimates based on qPCR and antibody-based approaches are in general not comparable among different HLA loci, since different primers/antibodies are often developed for each locus. Thus, the information for expression differences which we observe among HLA loci are usually not available from studies using these techniques. A few studies developed antibodies with comparable affinities for HLA-A/B/C, allowing comparison of expression between loci. Using such a method, Apps et al. [19] found an at least 12-fold reduction in HLA-C expression with respect to HLA-A and HLA-B as measured by flow cytometry, whereas the difference we documented based on RNA-seq was only twofold. However, the Apps et al. [19] study was restricted to a specific haplotype, and because our results show substantial within locus variation in expression levels, it is difficult to compare studies. Second, many additional layers of biological factors make a direct comparison across studies challenging. For example, the abundance of total protein, of proteins on cell surface or of mRNA in the cell, represent different molecular phenotypes associated with gene expression, given the existence of various post-transcriptional and post-translational regulatory processes [18,20,42,43,65,66]. As a consequence, there are good biological reasons to expect differences when expression is estimated at each of these levels. In addition, it is well known that gene expression varies among tissues and cell types [26,32], posing an additional challenge in comparisons between our findings and those previously reported in the literature, which are often based on analyses of different cell types or tissues. Third, regarding mRNA quantification, qPCR and RNA-seq each have their own sources of bias, and the finding that there are differences between studies using these approaches should not, on its own, provide definitive evidence that one or the other is inaccurate. Future studies will need to evaluate the degree to which these two technologies can be compared for quantification of HLA expression. Nevertheless, given the availability of qPCR based quantifications for HLA expression, as well as the interest in obtaining allele-level estimates of expression, we compared the expression for individual HLA lineages obtained by RNA-seq with those from qPCR from previous studies. We found that RNA-seq and qPCR show high concordance for the relative ordering of lineage-level expression for HLA-B and HLA-C (66% and 100%, respectively), whereas for HLA-A concordance was low, at 28%. Therefore results suggest that for HLA-B and HLA-C, when we examine lineages that are sufficiently different in their expression relative to one another (see S1 Text), the two methods are in high agreement, but further research is required to identify the sources of discordance between RNA-seq and qPCR for HLA-A. It is important to recall that these analyses rely on lineage-level expression estimates obtained by qPCR which are indirect, since they estimate locus-level expression (whereas RNA-seq data analysis with HLApers distinguishes among alleles), adding another source of differences among methods. Our perspective is that differences across studies call for further investigation, and highlight the potential biological variation that needs to be accounted for when comparing cell types, infection status, methodological approach, and molecule being quantified. In this study we present HLApers, and use simulations and empirical data analyses to show that it produces estimates for HLA expression from RNA-seq with high precision and accuracy. HLApers incorporates a Maximum-Likelihood estimator to deal with instances of reads mapping to multiple alleles or loci, following a strategy that has been widely used in RNA-seq studies [26,67]. Our pipeline can use different alignment strategies (e. g. STAR-Salmon [34,36] or kallisto [35]), but we show that the accuracy of expression depends on the sequences contained in the index, and is less dependent on the program used. The impact of using an HLA-personalized index on expression estimates varies markedly among loci. We found that for HLA-DQA1 there are large differences between gene-level expression estimates obtained using the reference transcriptome and those obtained using the HLA-personalized index. However, this difference is quite low for Class I genes, and intermediate for Class II loci other than HLA-DQA1 (Figs 2 and 4). However, even for the Class I genes, using the personalized index does result in changes in expression estimates with respect to the reference transcriptome. Therefore, we asked whether these changes have an impact on downstream analyses, such as eQTL mapping. When we compare the eQTLs identified using either HLApers or the reference transcriptome, 6 out of 8 loci had non-overlapping eQTLs. However, in most cases the same biological signal was being captured (S2 Table), and the p-values (Fig 5) and the causality (S5 Fig) of the eQTLs obtained with HLApers are only modestly better than when the reference transcriptome is used. We show that most of the eQTLs we identified (14/17) share a signal with previously reported regulatory SNPs (either experimentally validated SNPs and eQTLs (S3 Table) or CRD-QTLs (S5 Table) ). This indicates that, despite of the improvement in expression estimates that personalized pipelines can generate, larger sample sizes are necessary in order to identify eQTLs with greater probability of being causal, and to identify novel eQTLs with smaller effects. We note that the use of transformed cell lines (such as LCLs) may lead to expression profiles specific to this environment, and as a consequence some of the regulatory architecture underlying the expression of the HLA genes may not be shared with other tissues, cell types or treatments. The HLA-personalized approach provides expression estimates at the HLA allele level, which is not a product of standard RNA-seq pipelines. We integrate allele-level information with the eQTLs mapped for the genes, showing that the HLA allele is a relevant layer of information to understand the regulation of gene expression, because in some instances the regulatory architecture is linked to specific HLA alleles. This joint mapping of regulatory variants and assessment of expression of HLA alleles can illuminate the understanding of the HLA regulation, and contribute to disentangle specific contributions to disease phenotypes.
All data were obtained from third party sources and no additional ethical approval was required. In order to create the index, we downloaded 16,187 nucleotide sequences for 22 HLA loci (HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, HLA-G, HLA-DMA, HLA-DMB, HLA-DOA, HLA-DOB, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQB1, HLA-DRA, HLA-DRB1, HLA-DRB2, HLA-DRB3, HLA-DRB4, HLA-DRB5, HLA-DRB7, HLA-DRB8) from the International Immunogenetics/HLA database (IMGT) (release 3. 31. 0 available at https: //github. com/ANHIG/IMGTHLA). For many alleles, sequence data is not available for the entire coding region (e. g. , only for exons 2 and 3 for class I, and exon 2 for class II genes, which are called ARS exons). Because the lack of sequences for much of the coding region would cause the exclusion of many reads, including those mapping to the boundaries of the available exons, for each allele with partial sequence we used the available sequence to find the closest allele which has the complete sequence, and attributed the sequence from this allele. This is expected to introduce little bias in either the genotyping step (because ARS exons are the most polymorphic and sufficient to distinguish specific alleles) or the expression estimation (because the non-ARS sequence attributed is likely very similar to the real one). For the final index file, we replaced the HLA transcripts in the reference transcriptome (Gencode v25, primary assembly) with the HLA diversity described above. STAR’s module genomeGenerate, Salmon’s index, and kallisto’s index compile an index from these sequences. In order to select the alleles to be used in the HLA-personalized index, we observed that a simple procedure of selecting the 2 alleles with the largest number of estimated read counts after applying a zygosity threshold is sufficient to produce calls with accuracy of ≥ 95%. However, in order to avoid false homozygotes and false heterozygotes, we implemented additional steps. First, we selected the top 5 alleles and applied an intra-lineage threshold of 0. 25, meaning that only alleles which had at least 25% of the total expression in their lineage were considered for further steps. For each individual, we compiled an index containing only these (up to) 5 alleles and estimated their expression. We then determined if the individual was heterozygote at the lineage level by applying a threshold of 0. 15 on lineage expression levels. The lead allele from each lineage was selected to compose the genotype. A zygosity threshold of 0. 15 was applied to decide whether the genotype was heterozygous at the allele-level. For each locus, the reads mapped to the lead allele were removed, and another step of alignment and quantification was performed in order to determine if the second allele was real, or just noise due to extensive similarity to the lead allele. If the second allele had at least 1% of the locus read counts, it was kept, otherwise the genotype was considered to be homozygous for the lead allele. The thresholds described above were chosen because they maximized the concordance with the Sanger sequencing typings [6], while also minimizing the rate of false homozygotes and heterozygotes. We implemented two versions of the HLApers pipeline: (1) one using STAR (v2. 5. 3a) [34] to map reads followed by Salmon (v0. 8. 2) [36] to quantify the expression, and (2) using kallisto (v0. 43. 1) [35], which performs pseudoalignment and quantification. The quantification pipeline is structured in a two-stage process, first identifying the most expressed allele (s) at each HLA locus in order to infer the genotype which is present, and next quantifying expression for these inferred genotypes as well as for the rest of the transcriptome (Fig 1). Reads were aligned directly to the transcriptome. STAR alignments were passed to Salmon for quantification (module quant under alignment mode), whereas kallisto directly produces quantifications with the quant module. In both the HLA typing step (for which the index contains all HLA sequences in the IMGT database) and the quantification step (for which an HLA-personalized index is used), short reads can map to more than one locus, or more commonly to multiple alleles of the same locus (multimaps). The quantification methods we are using deal with multimaps by inferring maximum likelihood estimates optimized by an expectation-maximization algorithm to probabilistically assign reads to each reference in the index, and also include models to account for sequencing bias. For the mapping with STAR, we tuned parameters in order to avoid discarding multimaps and to accommodate mismatches. For quantification, we used all bias correction options available (–seqBias and –gcBias in Salmon, and –bias in kallisto). Polyester is an R package designed to simulate RNA-seq datasets [41]. We used the function simulate_experiment_countmat to simulate transcriptome data for 50 randomly chosen GEUVADIS individuals. Simulations were based on the read lengths and counts in the original data, with library sizes of 30 million reads, sampling without bias from a normal distribution of read start sites (with average fragment length of 250bp and sd of 25bp, and error rate of 0. 005). The code for the generation of the simulated datasets is available at https: //github. com/genevol-usp/hlaexpression/tree/master/simulation/data. We processed the simulated reads with STAR-Salmon to perform the quantifications using different indices: To investigate the relationship between quantifications and sequence divergence with respect to the reference, we used the R function adist to calculate the proportion of mismatches between the HLA alleles carried by the individuals and the alleles in the reference genome. We quantified HLA expression based on RNA-seq data for 358 European individuals, in samples of LCLs (Lymphoblastoid Cell Lines), made available by the GEUVADIS Consortium [24] (we excluded samples from the original dataset which are not in the 1000 Genomes phase 3). We performed expression quantification using the HLApers and reference transcriptome pipelines. We evaluated the reproducibility of the HLA quantifications by analyzing a subset of 97 individuals for which a replicate was available (https: //github. com/genevol-usp/hlaexpression/blob/master/geuvadis_reanalysis/replicates/data/write_sample_info. R). For the eQTL analysis, we used only autosomal genes which are expressed in a large proportion of samples, exploring the thresholds of TPM > 0 in at least 25%, 50%, or 75% of samples. In order to correct the expression data for technical effects, we sequentially removed the effect of the first 0 to 100 PCs and ran an eQTL analysis for each condition (S6 Fig). The configuration of thresholds and number of PCs which maximized the eQTL discovery (at FDR = 5%) was considered. This resulted in the use of genes expressed in ≥ 50% of samples (19,613 genes), and 60 PCs. The PCA analysis and data correction were performed with QTLtools v1. 1 [47], using the modules pca and correct respectively. For the genetic variant data, we used the 1000 Genomes Phase 3 biallelic variants, lifted to GRCh38 coordinates, after filtering for MAF ≥ 0. 05 in the individuals included in this study (6,837,505 variants in total). In order to control for population structure in the eQTL analysis, we ran a PCA on the variant genotype data and assessed the PCs which captured the structure. We used QTLtools pca requiring that variants should be at least 6kb apart. After visual inspection of the plots in S7 Fig, PCs 1–3 were used as covariates in the eQTL analysis. We used QTLtools cis to conduct the cis-eQTL analysis using the following model: PCA - corrected and standard normal expression ∼ SNPs + covariates (PCs for population stratification) The permutation pass was performed with 1000 permutations and a cis-window of 1Mb. P-values were computed by beta approximation and significance was determined by running the script runFDR_cis. R provided by QTLtools with FDR of 5%. Multiple eQTLs with independent effects on a particular gene were mapped with a conditional analysis based on step-wise linear regression (see Supplementary method 8 in [47]). The method automatically learns the number of independent signals per gene and provides sets of candidate eQTLs per signal. In order to investigate the putative function of the eQTLs we mapped, we investigated whether these eQTLs were present in ENCODE [53] regulatory elements annotated for LCLs. We used three types of functional annotations: open chromatin regions given by DNAse footprinting, transcription factor binding sites (TFBS) assayed by ChIP-seq, and histone modifications. We performed an RTC [48] analysis as described in [47] to investigate whether our eQTLs tagged the same causal variant as a GWAS variant or previously reported eQTL. We downloaded the GWAS catalog data (v1. 0. 1) from https: //www. ebi. ac. uk/gwas/api/search/downloads/alternative and selected associations with p-value < 10−8. We obtained the coordinates of recombination hotspots from http: //jungle. unige. ch/QTLtools_examples/hotspots_b37_hg19. bed. We applied the RTC module implemented in QTLtools (QTLtools rtc), selecting the HLA region only, using a D’ threshold of 0. 5, and turning on the conditional flag (–conditional) to test all independent eQTLs for a gene. Following the recommendation of Delaneau et al. (see Supplementary Note 7 in [47]), we considered that two SNPs tagged the same functional signal if the RTC score was >0. 9. To investigate whether there is a haplotypic coordination of expression at HLA, we used phased HLA genotype data to verify if there was more correlation of expression between alleles on the same haplotype than on different haplotypes (Fig 8). To also assess the phasing between HLA alleles and the eQTLs (Fig 6), we included the eQTLs mapped for each HLA gene in the phasing procedure, accounting for the fact that their phasing was already known from 1000 Genomes Project. In order to be conservative in the phase estimation, we used only the haplotype calls which were concordant between two approaches for phasing. First, we used PHASE [68] to determine the haplotype of each allele in the genotype, providing HLA allele designations and phased eQTL genotypes as input. Second, we checked the compatibility of the individual HLA genotypes with the phased SNP haplotypes from 1000 Genomes. Given all possible haplotypes for each individual at HLA-G∼A∼E∼C∼B∼DRA∼DRB1∼DQA1∼DQB1∼DPA1∼DPB1, we checked which combination of 2 haplotypes minimized the number of differences to the 1000 Genomes data to infer the haplotypes present in the individual. This resulted in 417 haplotypes completely concordant between the two approaches. The HLApers pipeline is available at https: //github. com/genevol-usp/HLApers. The entire analysis, including simulations, index compilation, quantification of expression, eQTL mapping, etc is available at https: //github. com/genevol-usp/hlaexpression.
The level at which a gene is expressed can have important influence on the phenotype of an organism, including its predisposition to develop diseases. One way to estimate gene expression is by quantifying the abundance of RNA. RNA-seq has become the method of choice to provide such estimates at the genomewide scale. However, the application of RNA-seq to HLA genes —key players in the immune adaptive response— has remained a rarely explored approach. This is due to the problem of mapping bias, which causes deficient read alignment at genes which are very polymorphic and different from the reference genome. This has motivated approaches that replace the single reference genome with personalized sequences, comprised of the individual’s specific HLA genotype. Here we explore the use of computational frameworks to obtain reliable expression levels for HLA genes from RNA-seq datasets. We present a pipeline in which the quantification of HLA expression is carried out using methods which account for HLA diversity, avoiding the biases of standard approaches. We then evaluate the impact of this form of quantifying HLA expression on downstream analyses. The pipeline also allows us to integrate information on eQTLs with expression levels at the HLA allele-level, which can help disentangle different contributions to disease phenotypes and help understand the regulatory architecture at the HLA region.
genome-wide association studies
variant genotypes
alleles
genetic mapping
genome analysis
research and analysis methods
sequence analysis
sequence alignment
bioinformatics
gene expression
genetic loci
haplotypes
heredity
database and informatics methods
genetics
transcriptome analysis
biology and life sciences
genomics
computational biology
human genetics
Expression estimation and eQTL mapping for HLA genes with a personalized pipeline