Quantifying HLA allele-specific expression using arcasHLA-quant

As a consistency check, we verified that gene-specific quantification levels obtained with arcasHLA-quant by summing minor and major allele expression for each HLA-I gene (HLA-A, HLA-B, and HLA-C) were consistent with expression levels as inferred through alternate methods [36] available on the TCGA portal (Pearson’s correlation coefficients in the range 0.94–0.96, p < 10−16 across tumor subtypes, Additional file 1: Fig. S1). arcasHLA-quant also yielded consistent results with HLApers [13] in terms of HLA-I gene quantification (Additional file 1: Fig. S2) and minor allele frequency inference (Additional file 1: Fig. S3), which was calculated as the number of reads supporting the allele with fewer reads over the total number of reads per HLA-I gene locus.

In order to evaluate the accuracy of HLA-I ASE loss calling, we took the TCGA normal cases as a negative control and the TCGA tumor cases with heterozygous loss of HLA-I genes at the DNA level as a positive control. We assumed that allelic loss at DNA level in the tumor would also be detectable at the level of mRNA expression. The positive control cases were manually validated LOHHLA calls with high confidence. Among 382 positive control cases, 324 were identified to have HLA-I ASE loss, which indicates an overall sensitivity of 84.8%. The positive control samples covered a wide range of tumor purity values (from 10% to > 90%). Among the 703 negative control cases which we took as a negative control, 15 were called as HLA-I ASE loss, indicating a specificity of 97.9% (688/703). The area under the receiver operating curve (AUC) is 0.921 (Additional file 1: Fig. S4).

HLA-I allele-specific expression loss is pervasive across tumor types

We first determined that HLA-I ASE loss is pervasive across TCGA tumor types (Fig. 2a): ASE loss was detected in every tumor type analyzed, most prominently in kidney chromophobe (KICH), with a frequency of 86%, followed by another eight tumor types that exhibited frequencies greater than 40%, including cervical squamous cell carcinoma (CESC, 58%), adrenocortical carcinoma (ACC, 50%), stomach adenocarcinoma (STAD, 48%), head and neck squamous cell carcinoma (HNSC, 46%), esophageal adenocarcinoma (ESCA, 45%), diffuse large B-cell lymphoma (DLBC, 43%), thymoma (THYM, 43%) and lung adenocarcinoma (LUAD, 42%). Five tumor types showed a markedly lower incidence of HLA-I ASE loss rates, including lower-grade glioma (LGG, 16%), glioblastoma (GBM, 12%), pheochromocytoma and paraganglioma (PCPG, 11%), uterine carcinosarcoma (UCS, 7%) and testicular germ cell tumors (TGCT, 6%). In order to limit ASE loss calling errors due to extremely low purity levels, we filtered out TCGA samples with purity below 0.1 (Additional file 1: Fig. S5). Overall, HLA-I ASE loss was attributable to HLA-A in 27% of cases, to HLA-B in 22% and to HLA-C in 25% of cases, with loss at all three genes occurring at a rate of 13% (Additional file 1: Fig. S6).

Fig. 2

Pervasiveness of HLA-I allele-specific expression loss. a Proportions of HLA-I ASE loss across TCGA subtypes (orange bars) as inferred using arcasHLA-quant. Blue (orange) bars represent proportions of cases where expression loss is (not) accompanied by somatic DNA loss, as inferred by LOHHLA on WES data. b HLA-I ASE comparison between tumor and normal cases in TCGA cohorts. HLA-I ASE is captured by the minimal raw minor allele frequency among the three HLA-I genes (minimal raw MAF). Numbers in the parentheses indicate the normal and tumor case numbers respectively. Only the TCGA cohorts with more than 10 normal cases are shown. Significance labels: “ns” or nothing labeled: p > 0.05; “*”: p < 0.05; “**”: p < 0.01; “***”: p < 0.001; “****”: p < 0.0001

We then investigated whether HLA-I ASE loss was accompanied by somatic DNA lesions (e.g., chromosomal or focal deletions) at the corresponding HLA-I locus. Using LOHHLA [12], we found that only a fraction of ASE losses showed evidence of DNA haplotype loss (Fig. 2a). The maximal proportion of DNA-to-expression-only loss was found in ACC (70%), while THCA had the smallest such proportion (under 4%). Our results suggest that a major proportion of HLA-I ASE loss in cancer may occur through epigenetic or other expression regulatory mechanisms instead of somatic DNA lesions.

We further compared HLA-I expression between the TCGA tumor samples and the paired normal samples, where available. We found that TCGA tumor samples had significantly lower minor allele frequency in most cohorts (Fig. 2b). We posited that the extensive allelic imbalance observed in many tumors was primarily due to the tumor component in each of the bulk samples analyzed.

We also studied the association between HLA-I ASE loss and microsatellite instability (MSI) in TCGA cohorts [37]. Due to the low reported frequency of MSI, only three cohorts (UCEC, COAD, and STAD) had more than 20 MSI cases (Additional file 2: Table S1). We found that the MSI cases consistently had higher HLA-I ASE loss frequency compared with non-MSI cases in the three cohorts: 53% vs 30% in COAD (odds ratio, OR = 2.65, p = 1.8 \(\times\) 10−4; Fisher’s exact test); 65% vs 52% in STAD (OR = 1.69, p = 0.075; Fisher’s exact test); and 38% vs 25% in UCEC (OR = 1.75, p = 5.5 \(\times\) 10−3; Fisher’s exact test).

Disrupted HLA-I expression can also result from the accumulation of somatic mutations in HLA-I [5]. However, somatic mutations in HLA-I are relatively infrequent in TCGA, varying from below 1% incidence in BRCA and GBM, to around 5% in BLCA, LUAD, and SKCM, and up to 10% in HNSC [7, 27]. Using a comprehensive list of HLA-I mutations in TCGA [31], we identified 58 nonsense HLA-I mutation events. Nonsense mutations in HLA-I genes are expected to lead to a severe imbalance of mRNA from the mutant allele, as compared to the wild-type allele, through nonsense-mediated decay. Consistent with this assumption, we found an enrichment of HLA-I ASE loss in the cases with nonsense HLA-I mutations (OR = 2.34, p = 2.8 \(\times\) 10−3; Fisher’s exact test; Additional file 2: Table S2a). Several factors may limit such an analysis: (a) the nonsense mutations may not be clonal, but only sub-clonal; (b) there may be issues with coverage in both DNA and RNA data; and (c) nonsense mutations result in premature termination of translation and it is not a priori clear that mRNAs are degraded or eliminated at such a rate so as to be detectable in RNA-seq. Despite these obstacles in detecting HLA-I ASE loss from nonsense mutations, we were still able to observe a significant enrichment (Additional file 2: Table S2a), further reinforcing the validity of our approach for assessing HLA-I ASE loss. We performed the same analyses on HLA-I missense mutations [27]. Interestingly, HLA-I ASE loss was also significantly enriched in the cases with HLA-B missense mutations (OR = 2.46, p = 5.4 \(\times\) 10−3; Fisher’s exact test; Additional file 2: Table S2b). The enrichment maintained significance when looking solely at the HLA-B missense mutations corresponding to contact residues (OR = 3.49, p = 0.012; Fisher’s exact test; Additional file 2: Table S2c), which was defined by [27] as the mutations in HLA positions that are in actual physical contact with the peptide. We also found that HLA-I ASE loss had significant enrichment in the cases with B2M missense mutations (OR = 2.35, p = 0.018; Fisher’s exact test; Additional file 2: Table S2d).

We assessed whether the presence of HLA-I ASE loss resulted in any survival time impact across tumor types. Specifically, a multivariate Cox regression stratified by tumor type was conducted as a pan-cancer analysis to investigate the clinical significance of HLA-I ASE loss. A total of 27 tumor sample features were taken as covariates, including age at diagnosis, tumor purity and ploidy estimates, and several immune-related and microenvironmental features (Additional file 1: Fig. S7). While we noticed that age at diagnosis, tumor ploidy, and macrophages were significantly associated with shorter overall survival, HLA-I ASE loss did not show a significant pan-cancer effect in this model (Additional file 1: Fig. S7). We also conducted univariate cox-regression models for every TCGA cohort but found no significant cohort-wide association after correcting for multiple hypothesis testing (Additional file 1: Fig. S8). However, HLA-I ASE loss showed a trend towards worse survival in KIRP (n = 264 cases; HR = 1.38, nominal p = 0.021; Additional file 1: Fig. S8a). The same trend was observed when the predicted neoantigen count was added to the regression model (HR = 1.43, p = 7.11 \(\times\) 10−3; Additional file 1: Fig. S8b). Additionally, HLA-I ASE loss showed a trend towards poorer prognosis in PDAC when the predicted neoantigen count was considered (n = 130 cases; HR = 1.24, nominal p = 0.058; Additional file 1: Fig. S8b). The difference in prognosis between cases with HLA-I ASE loss and those without HLA-I ASE loss was sharper among the patients with overall survival shorter than 24 months (Additional file 1: Fig. S9). Finally, among patients with heterozygous genotypes at all three HLA-I loci, there was a marked difference in prognosis as well (LR = 6.28; p = 0.01; log-rank test; Additional file 1: Fig. S10). Next, we evaluated these association trends with survival, observed at cohort level, by further dissecting our analysis into molecular subtypes and tumor stage.

HLA-I allele-specific expression loss is associated with decreased overall survival in the basal-like PDAC subtype

We interrogated the potential associations between HLA-I ASE loss in pancreatic cancer and the well-characterized classical/basal-like transcriptional subtypes [22, 38,39,40]. We found that ASE loss was present in both PDAC subtypes with some enrichment in the basal-like tumors (OR = 1.59, p = 0.21; Fisher’s exact test) which was not significant. Moreover, detection of HLA-I ASE loss in the basal-like subtype, but not in the classical one, was associated with worse survival (LR = 6.88, p = 0.01; log-rank test; Fig. 3a, b), suggesting the existence of a basal-like subcategory of PDAC characterized by HLA-I ASE loss and poorer prognosis. This result was basically consistent with alternate definitions of transcriptional subtypes in PDAC (Additional file 1: Fig. S11-12). We also noticed that HLA-I ASE detected at AJCC stage 2B had a significant association with shorter survival (LR = 5.19, p = 0.02; log-rank test; Additional file 1: Fig. S13). Additionally, HLA-I ASE loss was detected in earlier stages too, although without a noticeable effect (Additional file 1: Fig. S13).

Fig. 3

HLA-I allele-specific expression loss in basal-like and classical transcriptional subtypes of pancreatic adenocarcinoma. Survival curves in a TCGA-PDAC basal-like subtype; b TCGA-PDAC classical subtype; c CUMC-Epithelial PDAC basal-like subtype; and d CUMC-Epithelial PDAC classical subtype. Cases with ultra-low purity (< 0.1) were filtered out. Log-rank scores and p-values are indicated

Next, we validated widespread HLA-I ASE loss in an independent cohort of 96 laser-capture micro-dissected pancreatic ductal adenocarcinoma samples [21] where RNA-seq was performed separately on cleanly delineated epithelial and stroma compartments (CUMC cohort: CUMC-E for epithelial samples and CUMC-S for the stroma; see Cohort descriptions in Methods). Indeed, HLA-I ASE loss was strongly associated with the tumor epithelial compartment (OR = 3.95, p = 9.7 \(\times\) 10−6; Fisher’s exact test), which further supports our hypothesis that HLA-I ASE loss occurs in the cancer cells. Consistent with our previous TCGA analysis, HLA-I ASE loss was linked with shorter survival when detected in CUMC-E (LR = 3.97, p = 0.05; log-rank test; Additional file 1: Fig. S14). At the transcriptional subtype level in CUMC-E, HLA-I ASE loss was significantly associated with shorter survival in the basal-like tumors (LR = 8.03, p = 0.005; log-rank test; Fig. 3c), while no significant trend was observed in the classical subtype (LR = 0.00, p = 0.95; log-rank test; Fig. 3d). Combining the power of both the TCGA-PDAC and the CUMC-E basal-like cohorts strengthened the association with poor survival for HLA-I ASE loss cases (p = 1.02 \(\times\) 10−3, Fisher’s combined probability test). Altogether, our findings indicate that HLA-I ASE loss is a prognostic marker of shorter overall survival in the basal-like subtype of pancreatic ductal adenocarcinoma.

We performed a similar analysis in TCGA-KIRP stratified by subtypes (P-e.1a, P-e.1b, P-e.2, and P.CIMP-e) and stages (I, II, III, and IV), without identifying any new subtype associations (Additional file 1: Fig. S15).

HLA-I allele-specific expression loss is associated with poor outcomes in anti-PD-1 immunotherapy-treated metastatic melanomas

Finally, we hypothesized that HLA-I ASE loss may be a factor in the efficacy of immune checkpoint blockade immunotherapies. We revisited a previously published metastatic melanoma cohort [23] with pre- (n = 46) and on- (n = 29) Nivolumab therapy samples (see the “Cohort descriptions” section) and inferred HLA-I ASE loss as described before. Excluding samples with ultra-low tumor purity (below 0.1), we found ASE loss in both pre- and on-therapy samples with frequencies around 37% (Fig. 4). Furthermore, ASE loss was associated with worse overall survival regardless of whether it was assessed before or during therapy. The group with on-therapy ASE loss showed a slightly greater effect on prognosis (LR = 2.85, p = 0.09; log-rank test; Fig. 4).

Fig. 4

HLA-I allele-specific expression loss and poor outcomes to anti-PD-1 immunotherapy treatment for metastatic melanoma. Extensive HLA-I ASE loss was found in melanoma cohort pre- and on-treatment with Nivolumab. Heterozygous cases are those with heterozygous genotypes for all the three HLA-I genes

It has been previously reported that HLA class I homozygosity can reduce overall survival with immune checkpoint blockade [29] (Additional file 1: Fig. S16). As such, we also analyzed the impact of HLA-I ASE loss separately for individuals heterozygous at all three HLA-I genes. For these individuals, HLA-I ASE loss was associated with significantly worse prognoses (LR = 6.03, p = 0.01; log-rank test; Fig. 4) when expression loss occurred on-therapy (1 month after the start of therapy [23]). To a large extent, the association with decreased survival for heterozygous individuals was observed even before therapy, although these results are not as conclusive (LR = 2.76, p = 0.10; log-rank test; Fig. 4). The survival impact for the fully heterozygous cohort was also observed when we took neoantigen predictions into account (LR = 6.91, p = 0.01; log-rank test; Additional file 1: Fig. S17). Results with the full cohort (including cases with ultra-low tumor purity) showed the same trend towards worse prognosis, particularly for heterozygous individuals with on-treatment ASE loss (LR = 4.06, p = 0.04; log-rank test; Additional file 1: Fig. S18). Interestingly, among the heterozygous individuals with on-treatment samples and RECIST v1.1 [41]-defined response (n = 17), there were only 3 responders (complete or partial), none of which exhibited HLA-I ASE loss (OR = inf., p = 0.21; Fisher’s exact test). Among pre-treatment samples, HLA-I ASE loss resulted in slightly lower odds of responding to subsequent treatment (OR = 0.67, n. s.; Fisher’s exact test). In addition, survival associations were not explained by factors such as sample purity (Additional file 1: Fig. S19). In conclusion, our results highlight a potential significant clinical impact of HLA-I ASE loss on the efficacy of anti-PD-1 immunotherapy in metastatic melanoma.