HLA typing from RNA-Seq sequence reads

Abstract

We present a method, seq2HLA, for obtaining an individual's human leukocyte antigen (HLA) class I and II type and expression using standard next generation sequencing RNA-Seq data. RNA-Seq reads are mapped against a reference database of HLA alleles, and HLA type, confidence score and locus-specific expression level are determined. We successfully applied seq2HLA to 50 individuals included in the HapMap project, yielding 100% specificity and 94% sensitivity at a P-value of 0.1 for two-digit HLA types. We determined HLA type and expression for previously un-typed Illumina Body Map tissues and a cohort of Korean patients with lung cancer. Because the algorithm uses standard RNA-Seq reads and requires no change to laboratory protocols, it can be used for both existing datasets and future studies, thus adding a new dimension for HLA typing and biomarker studies.

Figure 1

HLA sequence variability. (a) The helices encoded by exons 2 (α1- chain: orange) and 3 (α2-chain: red) of the HLA class I alleles bind peptides (crystal structure 3OXR from PDB). (b) HLA locus is polygenic (A, B, C) and highly polymorphic within and across the three class I genes, showing the nucleotide variability at each position of a multiple sequence alignment of all HLA class I alleles, where 0 and 2 represent minimum and maximum variability, as defined using Shannon's Entropy (2 - Information Content; see Methods). (c) Example sequences are shown for four regions in seven HLA groups.

Figure 2

Workflow of seq2HLA. (a) RNA-Seq reads are mapped against the reference HLA sequences. (b) For each locus, the allele with the greatest number of reads is determined. (c) Reads associated with the first determined group are removed and the mapping procedure is repeated. (d) The second digital haplotype is called. (e) Zygosity is determined. (f) The genotype and the associated P-values are reported.

Figure 3

Sensitivity versus specificity for different mapping and technical parameters applied to the 50 Montgomery test samples. Bowtie mapping parameters include -a (report all mappings), -m1 (report only unique mappings), -v<0|1|2> (allow zero, one or two mismatches), using the initial reference dataset consisting of exons 1 to 3, plus 75 nucleotides of exon 4. Analysis of different technical parameters comprised varying the average number of reads per sample (1e⁶, 5e⁶ and 10e⁶) while mapping against the initial reference dataset with '-a -v1', shortening the read length (from 73 nucleotides to 30 nucleotides) and using only one of the read pairs to simulate a single-end read dataset. nt, nucleotide; PE, paired-end; SE, single-end

Figure 4

Summed expression of MHC class I and class II in normal tissues in RPKM units [19]. Seq2HLA was applied to the 50 nucleotide paired-end RNA-Seq dataset from the Illumina Body Map 2.0 project. HLA, human leukocyte antigen; RPKM, reads per kilobase of exon model per million mapped reads.