Inference of high resolution HLA types using genome-wide RNA or DNA sequencing reads

Abstract

Background: Accurate HLA typing at amino acid level (four-digit resolution) is critical in hematopoietic and organ transplantations, pathogenesis studies of autoimmune and infectious diseases, as well as the development of immunoncology therapies. With the rapid adoption of genome-wide sequencing in biomedical research, HLA typing based on transcriptome and whole exome/genome sequencing data becomes increasingly attractive due to its high throughput and convenience. However, unlike targeted amplicon sequencing, genome-wide sequencing often employs a reduced read length and coverage that impose great challenges in resolving the highly homologous HLA alleles. Though several algorithms exist and have been applied to four-digit typing, some deliver low to moderate accuracies, some output ambiguous predictions. Moreover, few methods suit diverse read lengths and depths, and both RNA and DNA sequencing inputs. New algorithms are therefore needed to leverage the accuracy and flexibility of HLA typing at high resolution using genome-wide sequencing data.

Results: We have developed a new algorithm named PHLAT to discover the most probable pair of HLA alleles at four-digit resolution or higher, via a unique integration of a candidate allele selection and a likelihood scoring. Over a comprehensive set of benchmarking data (a total of 768 HLA alleles) from both RNA and DNA sequencing and with a broad range of read lengths and coverage, PHLAT consistently achieves a high accuracy at four-digit (92%-95%) and two-digit resolutions (96%-99%), outcompeting most of the existing methods. It also supports targeted amplicon sequencing data from Illumina Miseq.

Conclusions: PHLAT significantly leverages the accuracy and flexibility of high resolution HLA typing based on genome-wide sequencing data. It may benefit both basic and applied research in immunology and related fields as well as numerous clinical applications.

Figure 1

PHLAT algorithm workflow. The algorithm consists of read mapping via Bowtie 2 to a reference sequence comprising the human genome and a plurality of genomic sequences of HLA alleles (I), selection of candidate alleles based on the number of mapped reads (II-IV), and log-likelihood scoring (V) over every pair of selected candidate alleles (e.g. a pair of a and b alleles). The pair of alleles with the best likelihood score is reported as the inferred HLA type at a given locus.

Figure 2

Analysis of frequently mistyped alleles. (A) The histograms illustrate the type (x-axis) and the number (y-axis) of the misidentified alleles at the HLA-DQA1 (left panel) and HLA-DQB1 (right panel) loci, summarized over the HapMap RNAseq, the 1000 Genome WXS and the HapMap WXS datasets. (B) Visualization of the mapped reads in one representative sample (subject NA12156, Additional file 1: Table S1) where the HLA-DQA1*03:01 allele is mistyped as the HLA-DQA1*03:03 allele. The mapped reads are shown around the single SNP position (chr6: 32609965, highlighted in between two vertical dashed lines) that distinguishes the two alleles. The hg19 reference sequence of the HLA-DQA1 gene is shown at the bottom of the panel. The nucleotide bases A, C, G, T are colored in green, red, blue grey and blue, respectively. The bases in the reads, if different from the reference sequence at the aligned positions, are visualized in the same color code. The pileup counts of the A, C, G, T bases at the highlighted SNP are 141, 117, 0 and 0, respectively. (C) The alignment of a 135-nucleotide segment from the HLA-DQA1*03:03 allele, noted as the query, with the HLA-DQA2 reference sequence in human genome hg19. The query sequence is simplified as a horizontal bar with only the mismatches indicated. The existing dbSNP record at the mismatch is labeled with a red vertical marker and the associated identification numbers (e.g. rs62619945) followed by a parenthesis indicating the major and the alternative base sequences. The alignment of the SNP that differ the DQA1*03:01 and DQA1*03:03 alleles is boxed.

Figure 3

Impact of read length, coverage and sequencing protocols on HLA typing accuracy. The plot summarizes the HLA typing accuracy of PHLAT using samples from the HapMap RNAseq (top panel), the 1000 Genome WXS (middle panel) and the HapMap WXS (bottom panel) datasets. Prediction accuracies are calculated by considering the sequencing data as either paired-end (close symbols and solid lines) or single-end (open symbols and dashed lines). The symbols represent the mean accuracy at four-digit resolution of the samples that are binned by their fold coverage at the HLA loci, with the error bars indicating the variance. The post-mapping fold coverage is calculated regarding to the CDS regions of the major class I and II HLA loci, excluding the reads suboptimal or not aligned to the candidate alleles. The smooth lines by spline interpolation illustrate the trend of the symbols.