A machine learning algorithm with subclonal sensitivity reveals widespread pan-cancer human leukocyt

Human leukocyte antigen loss of heterozygosity (HLA LOH) allows cancer cells to escape immune recognition by deleting HLA alleles, causing the suppressed presentation of tumor neoantigens. Despite its importance in immunotherapy response, few methods exist to detect HLA LOH, and their accuracy is not well understood. Here, we develop DASH (Deletion of Allele-Specific HLAs), a machine learning-based algorithm to detect HLA LOH from paired tumor-normal sequencing data. With cell line mixtures, we demonstrate increased sensitivity compared to previously published tools. Moreover, our patient-specific digital PCR validation approach provides a sensitive, robust orthogonal approach that could be used for clinical validation. Using DASH on 610 patients across 15 tumor types, we find that 18% of patients have HLA LOH. Moreover, we show inflated HLA LOH rates compared to genome-wide LOH and correlations between CD274 (encodes PD-L1) expression and microsatellite instability status, suggesting the HLA LOH is a key immune resistance strategy.

The success of immune checkpoint blockade (ICB) therapies confirms that the immune system is capable of recognizing and eliminating cells that present tumor neoantigens via major histocompatibility complexes (MHC)¹. However, a large percentage of patients do not respond to these therapies, driving the need for methods that can predict response². Elucidating the cause of resistance has proven more challenging than initially anticipated because of the complex tumor intrinsic and extrinsic mechanisms underlying ICB evasion^3,4.

Antigen presentation can influence the evolution of tumors^5,6,7. Each individual has up to six different classical class I human leukocyte antigen (HLA) alleles capable of presenting a unique set of antigens to the immune system. Tumor mutations are enriched in HLA and B2M genes, suggesting that the functionality of these genes is critical to mounting an immune response^8,9,10. Germline sequence diversity of HLA alleles also impacts tumor evolution; moreover, the impact of this evolution appears to be even more pronounced in the presence of ICB^11,12. Thus, somatic LOH in the HLA region has the potential to cause extreme reduction in HLA sequence diversity and is increasingly being recognized as a pathway of ICB escape¹³ and an indicator of response^14,15.

Detecting HLA LOH accurately from sequencing data is of interest given the growing ubiquity of tumor molecular profiling. Many copy number detection algorithms estimate tumor purity and ploidy to establish the expected relationships between b-allele frequency (sequencing depth ratio between two alleles— B/(A+B)), genomic sequencing depth and copy number deletions^16,17. However, general copy number detection algorithms used for detecting HLA LOH from sequencing data have proven unreliable because the polymorphic nature of the genes causes poor alignment to the reference genome. Moreover, the complexity of the sequence variation obscures the specific HLA allele that has been deleted, information crucial for neoantigen therapy design. Though genomic interrogation of the HLA region has been a challenge across applications, both alignment to allele-specific references and application of genome graphs have improved HLA typing and HLA somatic mutation identification^8,18,19. Although allele-specific alignment has been used for several HLA LOH detection algorithms, these algorithms are still based on tumor-only sequencing¹⁴ or on standard copy number variant (CNV) algorithms that fail to account for HLA-specific challenges such as differences in exome probe capture between alleles²⁰. Furthermore, CNV detection tools are notoriously poor for samples with low tumor purity and have trouble detecting subclonal deletions, raising concerns regarding the sensitivity and accuracy of these tools¹⁴. Thus, despite growing interest in the field, many studies still supplement deletions in regions surrounding the HLA region as a proxy for HLA LOH instead of using an HLA LOH-specific algorithm and forfeit information about the identity of the deleted allele^12,21,22.

Validating performance of HLA LOH detection algorithms has been an additional hurdle in the field. Two main approaches have been taken. First, concordance has been assessed between HLA LOH calls and copy number calls made by a standard CNV algorithm in regions flanking each HLA gene¹⁷. Second, primers have been designed to capture the regions surrounding the HLA genes and PCR has been used to test for copy number loss in patients predicted to have all HLA loci impacted. However, neither of these approaches validate the specific allele lost nor address the accuracy of calls for low tumor purity samples or samples with subclonal HLA LOH²⁰, hindering their utility as orthogonal validation approaches for securing clinical approval.

In this work, we design a machine learning algorithm, Deletion of Allele-Specific HLAs (DASH), to more accurately profile HLA allele-specific LOH by capturing the unique features of the HLA region and demonstrate the limit of detection for subclonal events and low tumor purity using sequencing data mixtures from cell lines. Moreover, we develop an allele-specific digital PCR assay to orthogonally validate DASH and probe the mechanistic impact of HLA LOH on peptide presentation with quantitative immunopeptidomics. Using a cohort of 610 patients spanning 15 tumor types, a range of tumor purity and all clinical grades, we apply DASH and find widespread HLA LOH, with the majority of tumor types having at least 20% of patients impacted. Finally, we find evidence for the evolutionary selection of HLA LOH, suggesting its importance in tumor development and immune escape.

DASH employs machine learning with specialized features to detect HLA LOH detection

Traditional copy number detection tools are unable to accurately detect allele-specific deletions in the HLA region; thus, we developed a machine learning approach to overcome the unique challenges this region presents. To create a dataset for learning genomic features important to accurately calling allele-specific LOH in HLA genes, we profiled 279 paired tumor and healthy normal samples (either adjacent tissue or PBMC) from patients. To type the HLA alleles, we used whole-exome sequencing on the ImmunoID NeXT Platform^TM, which has additional capture probes covering specific HLA alleles and has been validated for accurate HLA typing (Fig. 1a, Supplementary Tables 1 and 2). After determining the HLA types for each patient, we mapped all HLA-specific reads to each patient-specific HLA reference with stringent alignment parameters²³. From the alignment data, we calculated allele-specific sequencing depth and developed four features based on the genomic positions containing allele-specific germline variations (Fig. 1b). First, we wanted to capture a measurement of b-allele frequency, but we found that the sequencing depth in the data from adjacent normal tissue was often different between alleles. To address this issue, we divided the tumor b-allele frequency by the adjacent normal b-allele frequency, creating an adjusted b-allele frequency (Supplementary Fig. 1A). Second, we calculated an allele-specific sequencing depth ratio by dividing the tumor depth of an allele by the adjacent normal depth of the same allele and normalizing by the depth across the rest of the exome (Supplementary Fig. 1B). Third, we found that alleles with consistently lower sequencing depth than the alternate allele across the entire gene were likely to be deleted, whereas sporadically lower depth may be due to stochastic variation (p < 2.2e−14, paired T test with the null hypothesis that samples with and without HLA LOH would have the same consistency of depth across the gene, Supplementary Fig. 1C). Thus, we added a feature, consistency of sequencing depth, to capture this observation. Finally, in order to distinguish allelic imbalance driven by LOH from allelic imbalance driven by a large amplification of the alternate allele, we added a feature, total sequencing depth ratio, that captures the combined depth of the two alleles. We found that patients with HLA LOH have significantly lower total sequencing depth ratios (p = 0.0004, paired T test with the null hypothesis that patients with and without HLA LOH will have the same distribution of total depth ratios, Supplementary Fig. 1D).

a Tumor and adjacent normal samples were collected for 279 patients. DNA was extracted and exome sequencing was performed. HLA typing was derived from the normal sample and HLA somatic mutations were obtained from the tumor sample. Reads mapping to an HLA reference database were assembled and mapped onto the patient-specific HLA alleles. Horizontal gray lines denote positions where the homologous alleles differ from each other. b Features used to train the DASH algorithm include the adjusted b-allele frequency, sequencing depth ratio, total sequencing depth ratio, consistency of sequencing depth, tumor purity, tumor ploidy and deletion of flanking regions. Dashed gray lines represent the expected values for a sample without copy number change. All features were used to train an XGBoost prediction model. c A bar plot showing the sensitivity and specificity of HLA LOH detection by DASH and LOHHLA for samples in the 10-fold cross validation data set. Both algorithms were tested on the HLA-enhanced ImmunoID NeXT Platform. Only samples with at least 20% tumor purity were included. Source data are provided with this paper.

Full size image

To supplement the allele-specific features, we also included two patient-specific features, tumor purity and tumor ploidy¹⁶, that are identical across HLA-A, HLA-B, and HLA-C of a particular patient (Supplementary Fig. 1E, F). Finally, since 73% of copy number alterations causing HLA LOH are deletions of greater than one megabase (compared to 63% of genome-wide copy number alterations being greater than one megabase, Supplementary Fig. 1G), we reasoned that the genomic regions flanking the genes of interest can provide useful information to supplement the within-gene data. Thus, we utilized the whole-exome nature of our dataset to generate a final feature, deletion of flanking regions, that measures deletions in the 10 kb region surrounding each gene.

DASH outperforms competing algorithms in cross-validation evaluation

From the 279 patients, we collected these six features for 720 heterozygous HLA loci, excluding homozygous alleles because LOH cannot occur. To add training labels, we visualized all 720 alleles and manually curated deleted alleles and identified that 19.6% of the heterozygous loci had HLA LOH (Supplementary Fig. 2A, B, see the “Methods” section). Though generic LOH calling tools are unable to identify the specific allele that has been lost, we first wanted to assess the ability of a generic LOH calling tool to detect a deletion in HLA genes as a baseline on ImmunoID NeXT. Using Sequenza¹⁶, we found 92.9% specificity and 95.0% sensitivity (F1-Score = 0.848; Supplementary Fig. 3A). Then, using 10-fold cross validation and the features outlined above, we trained a gradient boosted regression (XGBoost) algorithm, which outperformed several other algorithms, to predict deleted alleles (Supplementary Fig. 3B). DASH reaches 98.7% specificity and 92.9% sensitivity while LOHHLA²⁰, another published HLA LOH algorithm, only achieves 94.3% specificity and 78.8% sensitivity (F1-Scores = 0.939 and 0.777, respectively; DASH AUROC = 0.939, DASH AUPRC = 0.940; Supplementary Fig. 3C–E). Though Sequenza is unable to detect the lost allele, we observe that Sequenza outperforms LOHHLA, but not DASH, in its ability to detect a deletion.

Of note, the majority of the incorrect calls have low tumor purity, highlighting the difficulty of accurately predicting HLA LOH at low tumor purity levels (Supplementary Fig. 3F). When samples with tumor purity below 20% are removed from the analysis, DASH continues to outperform LOHHLA at 99.7% specificity and 100% sensitivity compared to 94.3% specificity and 91.8% sensitivity (F1-Scores = 0.995 and 0.857, respectively; DASH AUROC = 0.990, DASH AUPRC = 0.992; Fig. 1c). The algorithm revealed that all seven of our features were independently contributing to the model, with deletion of flanking regions and adjusted BAF impacting the outcome most significantly (Supplementary Fig. 4A). We also found that none of DASH’s features alone could achieve this level of performance (Supplementary Fig. 4B).

All three of these algorithms were tested on the HLA-enhanced ImmunoID NeXT Platform, suggesting that performance may decrease with other exome platforms. One of the differences between the ImmunoID NeXT Platform and other exomes for this application is the boosted sequencing of the HLA loci (see the “Methods” section). ImmunoID NeXT provides an exonic region sequencing depth of 1586x prior to stringent post-processing via DASH and a full gene-sequencing depth of 587x after post-processing, while typical cancer exomes are performed at lower sequencing depths. For instance, example cohorts from dbGaP provide sequencing depths of 373x and 161x before and after post-processing, respectively (Supplementary Fig. 5A, B)^24,25. Thus, we evaluated the performance of DASH with decreasing sequencing depths across the HLA genes. We found that the F1-score drops slightly as sequencing depth decreases, dropping 0.06 from the ImmunoID NeXT depth to 100x depth (Supplementary Fig. 5C). Furthermore, the performance is highly dependent on tumor purity, with the F1-score dropping more dramatically at lower sequencing depths for samples with lower tumor purity (Supplementary Fig. 5D).

In silico cell line mixtures demonstrate the subclonal sensitivity of DASH

Though more study is required to understand the general timing of HLA LOH, it has been observed in some instances as occurring late in tumor progression as a resistance mechanism^13,20. Furthermore, tumor types that tend to be most responsive to ICB (e.g. lung, skin) also tend to produce lower purity samples²⁶. Finally, assessment of performance with manually curated labels has some limitations. For these three reasons, we thought that it was critical to assess the limit of detection of DASH with gold standard samples in both subclonal and low tumor purity settings. We sequenced 22 tumor-normal paired cell lines with ImmunoID NeXT and identified four with HLA LOH in at least one locus (see the “Methods” section, Supplementary Table 3, Supplementary Figs. 6–8). For three of these cell lines (CRL-5911, CRL-5922, and CRL-2314), we sub-sampled the tumor sequencing data and mixed it with complementary normal sequencing data to achieve simulated purity levels in replicates of ten (see the “Methods” section, Fig. 2a, Supplementary Data 1).

a A schematic showing the tumor and normal cell line mixing approach for simulating low-purity sample pairs. b and c Heatmaps showing the b specificity and c sensitivity of DASH to capture HLA LOH in simulated samples of differing purity and clonality. Dark blue denotes high sensitivity or specificity, light blue denotes low sensitivity or specificity and gray denotes no data. Source data are provided with this paper.

Full size image

With all three cell lines analyzed together, >98% sensitivity was observed at all tumor purity levels higher than 27% while >98% specificity was observed across all tumor purity levels (Supplementary Fig. 9A, B). While LOHHLA reached the same level of specificity as DASH, LOHHLA shows weaker sensitivity, with >97% above 35% tumor purity (Supplementary Fig. 9C, D). Notably, performance is much higher for both algorithms in two of the cell lines (CRL-5911 and CRL-2314) than the third (CRL-5922), with 100% sensitivity for tumor purities above 8% (Supplementary Fig. 10A–D, I–L and Q–T). Next, we mixed the HLA-mapping reads with the whole-exome reads across a range of different ratios to simulate the potential spectrum of tumor purities and subclonalities (see the “Methods” section). Both DASH and LOHHLA retained >99% specificity across all subclonality and tumor purity levels (Fig. 2b, Supplementary Fig. 9E). As expected, both their sensitivities decreased with lower purity and clonality; however, DASH still retained >97% sensitivity until the mixture of reads was at least 25% derived from the HLA LOH event (as observed in the tumor purity and HLA LOH clonality combination: 50% purity and 50% clonality). In contrast, LOHHLA was unable to recognize the vast majority of events that were <80% clonal, highlighting the accuracy and low detection limit of DASH to capture subclonal HLA LOH events (Fig. 2c, Supplementary Fig. 9F, Supplementary Fig. 10E–H, M–P, U–X).

Genomic validation of DASH in tumor samples using allele-specific digital PCR with patient-specific primer designs

After establishing DASH’s ability to detect HLA LOH in challenging samples, we sought to orthogonally validate our prediction ability with digital PCR (dPCR). Due to the highly polymorphic nature of HLA alleles, PCR primers and probes had to be independently designed for each pair of alleles in each patient (Fig. 3a, Supplementary Data 2). In addition to having to target the alleles unique to each patient, the probes were also designed to avoid targets on other patient-specific alleles. As a proof of concept and to determine the limit of detection of the dPCR assay, we designed primers for the predicted retained and lost HLA-C alleles for one of the cell lines previously interrogated with in silico mixtures (CRL-5911). Then, we mixed the tumor and blood (normal) cell line DNA to emulate decreasing levels of tumor purity in triplicate and performed dPCR with both allele-specific designs. Importantly, we also multiplexed the allele primers with RNase P (RPP25) primers as a diploid control. To ensure the specificity of the primers, we first confirmed that both the predicted lost and predicted retained allele copies (normalized by half of the diploid RNase P copies) resulted in a single copy in the normal sample with a tumor purity of zero (Fig. 3b, Supplementary Fig. 11A). Then, we compared the copy number of each allele in the tumor sample to the normal sample with a one-sided t-test and the null hypothesis that the tumor sample would have greater than or equal copy number than the normal sample. As anticipated, we found zero copies of the lost allele in the tumor sample, confirming the allele-specific LOH event, and dPCR sensitivity for 10% tumor purity and above, confirming the sensitivity and reproducibility of allele-specific dPCR as an orthogonal method (Fig. 3b, Supplementary Fig. 11A).

a Allele-specific genomic validation was performed using paired tumor and adjacent normal fresh frozen samples. DNA was extracted from each sample. Allele-specific primers were designed specifically for each patient. Digital PCR was performed on each sample to orthogonally determine the allele-specific copy number. b Bar plot showing the allele-specific copy number of the predicted lost allele, relative to RNase P, as measured by dPCR for cell line mixtures of 11 varying tumor purities, with n = 3 technical replicates examined over three independent experiments per tumor purity (tumor purity 2 replicate 3 only had 1 independent experiment; tumor purity 0 and 100 only have 1 technical replicate). The dashed line denotes the expected value for no change in copy number. Asterisks denote a statistically significant difference (p < 0.05) from the copy number in the normal sample with a one-sided Student Ttest. c Bar plots denoting the HLA allele dPCR copy number relative to the multiplexed RNase P for 8 patient tumor samples with n = 3 technical replicates. The alleles predicted by DASH to be retained are shown on the top plot while the alleles predicted to be deleted are shown on the bottom plot. The dashed gray lines denote the expected copy number of one if there are no copy number alterations. Asterisks denote samples with p-values < 0.05 as determined by a one-sided Student T test. d Bar plots showing the HLA allele dPCR copy number relative to the multiplexed RNAse P for 13 patient tumor samples where both alleles were predicted to be retained with n = 3 technical replicates. Asterisks denote samples with p-values < 0.05 as determined by a one-sided Student T test. 95% confidence intervals are shown in gray. Source data and exact p values are provided with this paper in the Source Data file.

Full size image

Since the detection of HLA LOH in the cell line with dPCR was successful, we profiled a set of 21 patient tumors: 4 confirmatory samples from the training data and 17 independent samples. Of these 21 samples, 8 were predicted to have HLA LOH by both DASH and LOHHLA, 1 was predicted to have HLA LOH by LOHHLA and not DASH, and 12 were predicted not to have HLA LOH by either algorithm. For one locus of each patient, we designed primers for the predicted retained and lost alleles (Supplementary Data 2) and performed dPCR on DNA from normal (adjacent) and tumor samples (Supplementary Fig. 11B–E, Supplementary Data 3). First, we ensured the specificity of the allele-specific primers (Supplementary Fig. 12A, B). We found that the RNase P-adjusted copy number for 30 of the 44 primers were highly specific and in close proximity to one copy; however, we excluded several alleles due to low specificity. Next, we evaluated the RNase P-adjusted copy number in the tumor samples. For the samples with predicted HLA LOH, we found a significant reduction of copy number in all eight of the predicted lost alleles by DASH and no reduction in copy number in all of the predicted retained alleles with specific primers in the tumor samples (Fig. 3c). For the samples predicted to not have HLA LOH by both algorithms, none of the alleles showed a significant reduction in copy number by digital PCR (Fig. 3d). Finally, for the sample predicted to have HLA LOH by LOHHLA but not by DASH (Patient K), both alleles appear to have a reduction in the tumor. However, genomic copy number evidence suggests an amplification leading to four copies of the RNase P control gene, suggesting that these dPCR results are likely false positives. Excluding samples from the training dataset (A, B, I, K), which includes the ambiguous results from patient K, all of the retained and lost alleles were identified correctly by both DASH and LOHHLA (100% sensitivity and 100% specificity). The patient-specific dPCR presented here represents a major technological advancement in allele-specific genomic HLA LOH validation.

Quantitative immunopeptidomics is limited as a validation approach

Mechanistically, HLA LOH is expected to effectively reduce the neoantigen load by eliminating surface presentation of neoantigens that would bind to specific HLA alleles. This mechanistic link has been demonstrated with organoids²², but it has not been shown in complex patient tumor samples. Thus, we sought to measure the quantitative changes in peptide presentation between adjacent-normal samples predicted by DASH to have or not have HLA LOH (Fig. 4a). We performed quantitative immunopeptidomics using TMT labeling on three patients with paired (normal and tumor) samples without any HLA LOH and on three patients with paired samples that were predicted to have HLA LOH (Supplementary Data 4). Across the six samples, we found robust peptide yields (1039–8319, average 4777.29 unique peptides), the expected length distribution (63.0–69.5%, average 65.7% of peptides having nine amino acids), and a high percentage of observed peptides predicted to bind to at least one of the patient-specific alleles (Supplementary Fig. 13A, B). Then, we mapped the peptides to their predicted respective alleles using an MHC-binding algorithm and looked for differences in peptide intensity for the peptides bound to lost, retained and homozygous alleles (Supplementary Figs. 14–16)²⁷. With all six samples taken together, we found that peptides predicted to bind to lost alleles had increased peptide intensity in tumor samples compared to normal samples for HLA-A, HLA-B and HLA-C alleles (Fig. 4b). However, the differences are very small and inconsistent across samples (Supplementary Fig. 13C). In fact, one of the samples that did not show the expected trend was previously validated as having HLA LOH with digital PCR (patient C), suggesting that other factors also influence peptide presentation. While a lack of consistently lower peptide intensities for lost alleles in tumor samples could be due to low tumor purities (Supplementary Fig. 13D), variable differences in allele-specific expression (Supplementary Fig. 13E) or the influence of immune infiltration, we wanted to ensure that the method was robust to variabilities in immunoprecipitation and the limit of detection using mass spectrometry. Thus, we performed the same quantitative immunopeptidomics experiment with the CRL-5911 paired cell line used in the previous experiments. Unfortunately, the tumor cell line had very low expression of MHC and 10-fold higher peptide intensity values in the normal cell line than the tumor cell line. Due to these limitations, further work will be needed to functionally demonstrate the impact of HLA LOH on peptide presentation in patient tumor samples.

a Functional immunopeptidomic validation was performed using paired tumor and adjacent normal fresh frozen samples. HLA-B2M complexes were purified from each sample and peptides were gently eluted. Peptides from each sample were labeled with TMT tags and measured using quantitative mass spectrometry. b Box plots showing the log2 fold change of peptide intensity between lost, kept, and homozygous alleles across all six patients. HLA-A is represented by 1490 peptides in the ‘lost’ category, 6613 peptides in the ‘kept’ category, and 2814 peptides in the ‘homozygous’ category. HLA-B is represented by 3823 peptides in the ‘lost’ category and 9910 peptides in the ‘kept’ category. HLA-C is represented by 703 peptides in the ‘lost’ category, 3421 peptides in the ‘kept’ category and 620 peptides in the ‘homozygous’ category. The center of the box denotes the median value, the box denotes the quartiles and the whiskers denote the remainder of the distribution apart from outliers. Statistical significance was assessed using a two-sided student T-test with asterisks denoting p-values < 0.05. Source data are provided with this paper.

Full size image

HLA LOH is pervasive across tumor types and exhibits signs of immune pressure on tumors

To assess the pervasiveness of HLA LOH as a potential immune escape mechanism, we applied DASH to 610 tumors across 15 tumor types (see the “Methods” section). The tumors from the patients were predominately stage II and III, with a median tumor purity of 31% and a median MATH²⁸ score (mutant-allele tumor heterogeneity) of 32.8 (Supplementary Fig. 17A–C). The fraction of patients with at least one incidence of HLA LOH ranged from 40% of patients in head and neck squamous cell carcinoma (HNSCC) to 4% of patients in liver cancer (Fig. 5a, Supplementary Table 4). We found a similar incidence of HLA LOH in non-small cell lung cancer adenocarcinoma (NSCLC-A) as compared to a previous study (24% and 29%, respectively), but found a significantly lower incidence in non-small cell lung cancer squamous cell carcinoma (NSCLC-SCC) (28% and 61%, respectively)²⁰. Of note, the overall frequency we report is similar to another study with a very large cohort size¹⁴. Moreover, we show that patients with HLA LOH often lose all three genes and more rarely lose only one gene or two genes (76.4%, 21.1%, and 5.5% of patients, respectively; Fig. 5b) and that HLA-A and HLA-B are preferentially lost compared to HLA-C (Supplementary Fig. 17D). While some specific alleles appear to be lost much more frequently than other alleles, this trend is likely due to low starting frequencies of these alleles (Supplementary Fig. 17E).

a Bar plots denoting the number of patients and the frequency of HLA LOH in each tumor type cohort. Only cohorts with at least 10 patients are shown (see source data for precise n of each cohort). 95% confidence intervals are shown in dark gray. b A bar plot showing the number of patients with 1, 2, or 3 genes impacted by HLA LOH. Only patients that are fully heterozygous across HLA-A, HLA-B, and HLA-C are shown. c Boxplots showing the distribution of the fraction of each genome impacted by LOH. Each tumor type is divided into patients with HLA LOH and without HLA LOH. Only tumor types with at least 10 patients impacted by HLA LOH are shown (with n = 88 for NSCLC-A; 80 for Colorectal; 31 for NSCLC-SCC; 55 for Kidney; 38 for Bladder; 30 for Pancreatic; 20 for HNSCC). The center of the box denotes the median value, the box denotes the quartiles and the whiskers denote the remainder of the distribution apart from outliers. Statistical analyses are performed with two-sided Mann–Whitney U tests and are Bonferroni corrected (p = 7.2e−20). d A scatter plot showing the relationship between the average fraction of the genome impacted by LOH and the frequency of HLA LOH in each tumor type. The gray dashed line denotes x = y. e and f Box plots denoting the average difference for patients without HLA LOH (green) and patients with HLA LOH (blue) for e neoantigen burden (n = 468 patients) and f CD274 (PD-L1) expression (n = 605 patients). The center of the box denotes the median value, the box denotes the quartiles and the whiskers denote the remainder of the distribution apart from outliers. Some outliers are excluded from the plot in order to focus on the majority of the distribution. Statistical significance was assessed using a two-sided Mann–Whitney U test. g Bar plot denoting the average difference for patients without HLA LOH (green) and patients with HLA LOH (blue) for percentage of microsatellite sites with instability (n = 576 patients). Source data are provided with this paper. NSCLC-A non-small cell lung cancer adenocarcinoma, NSCLC-SCC non-small cell lung cancer squamous cell carcinoma, HNSCC head and neck squamous cell carcinoma.

Full size image

Though high frequencies of HLA LOH in specific tumor types are of interest due to impairment of the antigen presentation pathway, the high frequencies alone do not necessitate an evolutionary advantage of the LOH event. For example, we found that patients with HLA LOH have significantly higher estimated rates of LOH across their genome, suggesting that some LOH in the HLA region may happen by chance (pan-cancer p = 7.2e−20 with the null hypothesis that patients with and without HLA LOH would have the same rate of LOH across their genome, Fig. 5c). To investigate if HLA LOH frequencies across cancer types would occur by chance, we compared the average estimated rate of LOH across the genome with the frequency of HLA LOH in a given tumor type cohort. If LOH was randomly occurring in the HLA region, we would expect to observe the rate and frequency to be similar. However, we see that nearly all tumor types individually and collectively have a higher frequency of HLA LOH than genome-wide LOH (Fig. 5d, Supplementary Fig. 18A). While this difference is small for some tumor types, we observed that colorectal cancer and HNSCC had substantial enrichment of HLA LOH, suggesting that HLA LOH may provide a greater evolutionary advantage in these tumor types than others. Alternatively, HLA may be more prone to deletion than the rest of the genome.

Since HLA LOH impacts the ability of a tumor cell to present antigen on the cell surface for recognition by the immune system, we reasoned that tumors with a greater ability to display neoantigens would be under higher selective pressure to incur HLA loss. Considering patients with high HLA evolutionary diversity (HED) are able to present a larger immunopeptidome and respond better to checkpoint inhibitors¹², we explored if they were more susceptible to HLA loss. However, we did not see this trend in our data (Supplementary Fig. 18B). Furthermore, high mutation rates and neoantigen burdens may also present pressure for cells to lose HLA. Indeed, we found correlations between mutation burden and HLA LOH pan-cancer (p = 0.006 with the null hypothesis that patients with and without HLA LOH would have the same distribution of mutation burdens, Supplementary Fig. 18C) and between neoantigen burden and HLA LOH pan-cancer (p = 0.02 with the null hypothesis patients with and without HLA LOH would have the same distribution of neoantigen burdens; Fig. 5e, Supplementary Fig. 18D). Interestingly, we did confirm the “goldilocks effect” for mutation burdens with intermediate mutation burden tumors having the highest rate of HLA LOH¹⁴, suggesting very high mutation rates may pose negative selection against HLA LOH (Supplementary Fig. 18E). Furthermore, we also corroborate correlations between HLA LOH and CD274 expression (encodes PD-L1) and microsatellite instability (MSI) status (p = 0.006 and 0.006, respectively, with the null hypotheses that patients with and without HLA LOH would have the same distribution of CD274 expression and MSI status, Fig. 5f, g, Supplementary Fig. 18F, G). The relationship between MSI status and HLA LOH suggests that some proportion of the HLA LOH events may be driven by genomic instability instead of immune pressure. We further hypothesized that HLA LOH may exert selective pressure on the mutational process as tumors develop. Indeed, we do observe that more neoantigens are predicted to bind to lost HLA alleles than their homologous counterparts (Wilcoxon rank sum test, p = 0.01 with the null hypothesis that the same distribution of neoantigens would be predicted to lost and retrained HLA alleles), suggesting that HLA LOH is a key mechanism for the selective exposure of antigen to the immune system.

Accurate detection of HLA LOH is critical for its use as a biomarker for cancer immunotherapy. Although allele-specific alignment improves current heuristic methods²⁰, variability in exome capture across alleles and the relatively short sequences of the HLA genes introduce additional hurdles for identifying changes in copy number. Moreover, the field lacks methods to accurately and comprehensively determine the limit of detection, sensitivity and specificity of HLA LOH detection algorithms.

To address these limitations, we have implemented a machine learning approach for HLA LOH detection, incorporating specialized features that result in improved performance. The features include (1) a modified b-allele frequency equation to account for differences in probe capture, (2) consistency of sequencing depth to capture variance across the allele, (3) total sequencing depth ratio to detect hyper amplifications and (4) information from the regions surrounding the HLA genes because the vast majority of HLA LOH events are large deletions. Though DASH was trained using sequencing data from the ImmunoID NeXT Platform, our downsampling experiment demonstrates that DASH can be applied to datasets with lower sequencing depths with a limited decrease in performance.

To address the lack of well-developed validation approaches for HLA LOH, we developed two methods to comprehensively evaluate algorithmic performance. Using in silico cell line mixtures, we utilized DASH’s high sensitivity and specificity across a range of tumor purities and subclonality levels and demonstrated superior ability to detect subclonal HLA LOH events than other available tools. To assess performance in actual tumor samples, we developed a PCR-based, allele-specific approach with sensitivity to detect subclonal HLA LOH. Of note, the scalability, sensitivity and allele-specific nature of this approach lends credibility to its potential clinical application. Moreover, we tested the applicability of quantitative immunopeptidomics to observe the functional impact of HLA LOH. However, more work must be done to test the sensitivity and robustness of this method as well as measure the impact of other biological factors, including infiltrating immune cells, variable expression or translation and peptide loading conditions, that may disrupt the signal.

Our HLA LOH prevalence data shows that a large percentage of patients are impacted by HLA LOH in several tumor types. Though NSCLC was previously reported to have a high incidence of HLA LOH²⁰, we also identified a large fraction of HLA LOH in cervical cancer (38%) and HNSCC (40%). NSCLC and HNSCC are also known for having a high mutational burden; however, we observe only 14% of patients with HLA LOH in melanoma, which also has a high mutational burden²⁹. Cervical cancer is strongly associated with human papillomavirus (HPV), which may play a role in the high frequency of HLA LOH. Of note, we observed a lower frequency of HLA LOH in NSCLC-SCC than one previously published study²⁰ and a similar frequency to another study¹⁴. Differences in tumor stage, clonality and treatment status may contribute to these differences. Larger cohorts and comparisons of tools within the same dataset are needed to understand all underlying covariates. Our analyses also revealed that patients tend to lose more than one HLA-allele at a time, potentially having stronger implications on tumor evolution.

Our study has several limitations. Larger cohorts of patients will allow for more scalable validation and increase confidence in the accuracy of the manual annotations used for training. Further, we have explored the detection of the functional impact of HLA LOH with quantitative immunopeptidomics. However, more work must be done to evaluate the robustness of this approach and impact of other biological factors. Finally, detecting HLA LOH alone does not capture all allelic-imbalance driven immune escape mechanisms. Methods to detect significant amplification of a single allele, potentially allowing the drowning out peptides from the homologous allele, and allele-specific expression are also needed.

In conclusion, DASH is a validated HLA LOH algorithm. We demonstrated DASH’s ability to sensitively detect subclonal events in samples with low tumor purity, enabling comprehensive profiling of widespread HLA LOH across tumor types.