The Riaz and Hugo bulk RNA datasets used in this study are available under BioProject accession numbers PRJNA356761 (Riaz et al., 2017) and PRJNA312948 (Hugo et al., 2016), respectively. Raw RNA-Seq reads were downloaded and processed using an in-house RNA-Seq mapping and analysis pipeline to calculate gene-level TPM values.

Pre-calculated gene-level and transcript-level TPM values for the TCGA Pan-cancer cohort for 33 cancer types were downloaded from the UCSC Xena data pages (Goldman et al., 2020).

Pre-calculated gene-level and transcript-level TPM values for 256 healthy tissues were downloaded from the Human Protein Atlas (HPA) (Uhlen et al., 2010).

Pre-calculated gene-level and transcript-level TPM values for 54 healthy tissue subtypes were downloaded from The Genotype-Tissue Expression (GTEx) project data portal (Carithers and Moore, 2015). Median TPM values were calculated for each of the 31 main tissue types.

Pre-calculated gene-level and transcript-level TPM values for 1,019 cell lines were downloaded from the Cancer Cell Line Encyclopedia (CCLE) (Ghandi et al., 2019).

All datasets were downloaded in July 2022.

The Trolle dataset consisted of 15,524 non-redundant HLA class I ligands eluted from mono-allelic HeLa cells transfected with five different HLA class I alleles (Trolle et al., 2016). This dataset was downloaded from the IEDB (Vita et al., 2019) under the accession number 1000685 (http://www.iedb.org/subID/1000685).

The Abelin dataset contained 22,310 non-redundant eluted ligands from mono-allelic B721.221 cells transfected with 16 different HLA class I alleles. The dataset was retrieved from the supplementary materials of the original publication (Abelin et al., 2017). Abelin et al. also provided matched RNA-Seq data for four replicates under BioProject accession number PRJNA360601. Raw RNA-Seq reads were downloaded and processed using an in-house RNA-Seq mapping and analysis pipeline to calculate gene-level TPM values. Median TPM values of the four replicates were used.

The HLA Ligand Atlas consisted of tissue-specific HLA ligands from 23 healthy tissue types (Marcu et al., 2021). The dataset contained 223,246 non-redundant peptides and 675,346 peptide tissue pairs. We downloaded the data from the HLA Ligand Atlas data pages (downloaded in September 2022).

The Shinkawa dataset contained 2,352 non-redundant HLA class I eluted ligands from a HCT15/β2 cell line. The dataset was retrieved from the supplementary materials of the original publication (Shinkawa et al., 2021).

The Pyke dataset contained 34,090 ligands eluted from mono-allelic K562 cell lines transfected with 25 different HLA alleles. The Pyke Cancer dataset contained 31,660 ligands eluted from 12 tissue samples of colorectal and lung cancer patients. Both datasets were retrieved from the Supplementary Material S1 of the original publication (Pyke et al., 2021).

The Sarkizova dataset contained 140,918 eluted ligands from B721.221 cells transfected with 79 different HLA class I alleles. The dataset was retrieved from the supplementary materials of the original publication (Sarkizova et al., 2020).

The backend of pepX is a PostgretSQL database that is populated with all possible 8-15mers from the human proteome linked to TPM data from several publicly available databases.

A partial entity-relation diagram of the core tables in the pepX database is shown in Figure 1. The proteome tables of the database were populated by retrieving all possible 8-15mers from the GRCh38 Ensembl proteome (release-106). All peptides are linked to their associated protein sequences, as well as gene, transcript, and protein identifiers. The ‘gene’ table includes the unique Ensembl gene (ENSG) identifiers, their gene symbols, and the number of proteins encoded by the gene. The ‘gene2tx2protein’ table contains one row per protein/transcript and maps the Ensembl gene identifier to their corresponding protein (ENSP) and transcript (ENST) IDs. This table also includes the full protein sequence. The ‘peptide2protein’ table maps each 8-15mer to their corresponding Ensembl protein ids, keeping note of the zero-indexed start position of the kmer within the full protein sequence. This schema allows pepX to quickly lookup which genes and transcripts are linked to a given peptide.

We included expression data from six public databases in pepX, namely HPA, GTEx, TCGA, CCLE, and RNA-Seq data of a B721.221 cell line (Abelin et al., 2017) as well as a HeLa cell line (Cantarella et al., 2019). Each has several subtypes, which are cancer types in the case of TCGA (e.g. BRCA, COAD, PAAD), tissue-types in the case of HPA and GTEx (e.g. Skin, Stomach, Thyroid Gland), and cell-lines in the case of CCLE (e.g. HELA_CERVIX, HCC56_LARGE_INTESTINE). These subtypes and their external dataset source are listed in the ‘expression_dataset’ table, each with Boolean values to indicate if gene- and transcript-level data are available (Figure 1). All entries in the ‘expression_dataset’ table contain a unique ‘dataset_id’ that is used to map TPM values in the ‘gene_TPM’ and ‘transcript_TPM’ tables back to the associated dataset. The ‘gene_TPM’ table maps TPM’s from a given subtype to an Ensembl gene id, while the ‘transcript_TPM’ table maps TPMs to the Ensembl transcript id. This structure allows pepX to quickly grab all TPMs linked to a given gene or transcript for a given study.

Several views have been created that efficiently perform the joins necessary to provide peptide-level and gene/transcript-level TPM output, given a list of peptides and dataset ID as input. In addition to the raw TPMs, scaled TPMs are calculated for each of the source proteins as: s c a l e d t p m = T P M ∗ #   o f   p r o t e i n s   f o r   t h e   g e n e   i n   w h i c h   t h e   p e p t i d e   i s   f o u n d # o f   t o t a l   p r o t e i n s   f o r   t h e   g e n e

For each peptide, the total TPM, median TPM, and maximum TPM are also provided.

R/Bioconductor was used for all statistical analyses. The following significance levels were used in all figures: ns: p > .05, *: p≤.05, **: p≤.01, ***: p≤.001, ****: p≤.0001. All statistical tests are paired Wilcoxon tests, unless otherwise indicated.

Reads mapping to tRNA, rRNA, adapter sequences, and spike-in controls were filtered with Bowtie 2 (v2.1.0). Remaining reads were mapped go the GRCh38 reference genome with Gencode v27 annotations using STAR (v2.6.1). Low complexity reads (DUST >4) were removed from the BAMs with PRINSEQ Lite (v0.20.3) before counting reads with FeatureCount (v1.6.5).

The TCGA program sequenced thousands of tumor and matched normal samples spanning 33 cancer types. For each cancer type, the number of patients analyzed varies, from 1,211 patients with breast cancer (TCGA-BRCA) to 45 patients with Cholangiocarcinoma (TCGA-CHOL, Supplementary Figure S1). This means when looking for the expression of a gene of interest in a specific cancer type, there are several expression values from different patients which need to be aggregated and summarized to derive one expression value per gene.

We sought to investigate what the best method is to aggregate TPM values from different patients and to determine how well the RNA-Seq data from TCGA correlates with patient-matched RNA-Seq. We obtained patient-matched RNA-Seq data from two melanoma studies published by Hugo et al. (Hugo et al., 2016) and Riaz et al. (Riaz et al., 2017) and compared the TPM values to the TCGA skin cancer samples (TCGA-SKCM). We aggregated the data over the 470 TCGA-SKCM patients by calculating the mean, median, and geometric mean TPM for each gene. Given the statistical background of how these metrics are calculated, it is expected that the values can significantly vary. The TPM values for PD-1 (PDCD1), for example, range between 0 and 60 in the TCGA-SKCM cohort (Figure 2A). The mean is 4.5, the median is 1.6, and the geometric mean is 1.3. For each patient in the Hugo and Riaz datasets, we considered all genes and calculated how well the TPMs correlate to the TCGA-SKCM mean, median, and geometric mean. We found that the values of the three metrics significantly vary from each other (Kruskal-Wallis test, p < .0001) and that the median TCGA-SKCM correlates best with patient-specific TPM values in both datasets (Wilcoxon test, p < .001, Figure 2B).

To determine the overall correlation between the patient-specific TPM values and the TCGA-SKCM median TPM values, we combined the Hugo and Riaz datasets and found a significant correlation (Spearman correlation coefficient r = .82). To get a better overview of the correlation, we separated the TPM values for each patient into ranges for both the patient-specific and the TCGA median TPM values to generate a 2-dimensional matrix. We then analyzed each TPM range combination and calculated the fraction of genes expressed within the corresponding TPM ranges (Figure 2C). We found that 67% of genes are expressed at low levels in both the Hugo and Riaz cohort and in the TCGA-SKCM cohort, with TPM values of <1. The majority of genes are expressed at similar levels in both cohorts, as demonstrated by an enrichment of genes on the diagonal in Figure 2C. We observed similar results when we analyzed a smaller set of in-house patients with six different cancer types (Supplementary Figure S2).

Taken together, these findings show that TPM values from a public database like TCGA are suitable for estimating gene expression in a patient sample if patient-specific RNA-Seq is not available.

We developed pepX, a tool for estimating a peptide’s expression level based upon the source antigen(s) in which it is contained. pepX takes a list of peptides and a public dataset identifier as input and returns the expression level of each protein the peptide was found in. pepX also provides aggregated expression levels for peptides that can be retrieved from multiple transcripts and proteins. The expression levels can be retrieved from a number of public databases, including The Cancer Genome Atlas (TCGA), CCLE (The Cancer Cell Line Encyclopedia) (Ghandi et al., 2019), HPA (The Human Protein Atlas) (Uhlen et al., 2010), and GTEx (The Genotype-Tissue Expression Project) (Carithers and Moore, 2015). These datasets provide expression values on the gene level as well as the transcript level. We used pepX to investigate different ways of estimating peptide abundance and the differences in using gene-level and transcript-level TPM values.

As illustrated in Figure 3A, it is possible that the exact same peptide can be found in different proteins encoded by different genes (e.g., Peptide A in Figure 3A can be retrieved from two proteins of Gene A and from one protein of Gene B). To analyze the extent of this, we considered the set of unique peptides in the HLA Ligand Atlas (n = 223,246) (Marcu et al., 2021) and investigated the number of possible source proteins for each peptide. We found that 88% of peptides can be retrieved from protein sequences corresponding to exactly one Ensemble gene id. We investigated the remaining 12% of peptides that could be retrieved from different gene ids and found that, for the majority of cases (96%), the corresponding genes belonged to the same gene family. As gene families are formed by duplication of a single original gene, genes that are categorized into families usually share nucleotide and protein sequences. It is thus not surprising that a peptide can occur in multiple proteins that are encoded by genes that are part of a gene family. It is, however, not clear how the abundance of such peptides should be measured, as there are several options: 1) using the median TPM of all genes, 2) using the maximum TPM among all genes, or 3) summing up the TPM values of all genes.

We used peptides from the HLA Ligand Atlas for validation and analysis of pepX performance. The HLA Ligand Atlas contains tissue-specific HLA ligands from 23 healthy tissue types (Marcu et al., 2021). The dataset consisted of 223,246 non-redundant peptides and 675,346 peptide tissue pairs. We wanted to investigate how the abundance level of these peptides could be best estimated using gene-level TPM values. We generated a set of decoy peptides by randomly selecting length-matched peptides from the human proteome and assessed the performance of different metrics in distinguishing true ligands from the set of random decoy peptides. The Human Protein Atlas was used as the expression dataset by matching the tissue types represented in the HLA Ligand Atlas ( Supplementary Table S1 ).

We evaluated all three options for aggregating expression values across different source antigens: 1) using the median TPM of all genes, 2) using the maximum TPM among all genes, or 3) summing up the TPM values of all genes. We used pepX to identify from which proteins each peptide could be derived from and retrieved expression levels of those source proteins. We evaluated the performance of each TPM aggregation method in distinguishing peptides from the HLA Ligand Atlas from the set of decoy peptides. The Area under the Receiver Operating Characteristics (ROC) Curve (AUC) was used to measure performance. We found that the values of the three metrics significantly vary from each other (Kruskal-Wallis test, p = .01034) and that with a mean AUC of .805, summing up TPM values (total TPM) performs best, followed by using the maximum TPM among all genes (mean AUC = .803) and using the median TPM of all genes (mean AUC = .785, Wilcoxon test, p < .0001, Figure 3B).

Another detail we wanted to investigate was that each gene can be transcribed into different transcripts and thus translated into different proteins. The different transcripts correspond to different splice variants that are found in different tissue types, developmental stages, etc. However, the expression values of the different transcripts of a gene are collapsed into a single expression value when generating gene-level expression data. We hypothesized that, when using such gene-level TPM data, it might be important to consider in how many of a gene’s transcripts the peptide occurs (e.g., Peptide B in Figure 3A occurs in three different proteins from Gene A, while Peptide C occurs only in one). We developed a ‘scaled TPM’, which considers the number of proteins in which the peptide is found and the total number of proteins for the gene (detailed in the Methods section). However, with a mean AUC of .802, this scaled TPM did not improve the performance in predicting peptides from the HLA Ligand Atlas (Figure 3C).

To validate these findings, we gathered additional ligand elution datasets and matched them to their corresponding expression datasets (Figure 3D). Again, the total gene TPM significantly outperformed the median gene TPM (Wilcoxon Test, p < .05, Figure 3E). The total gene TPM was also slightly higher than the maximum TPM in these datasets; the difference was, however, not significant. Also, the values of the three metrics were distributed in a similar way (Kruskal-Wallis test, p = .4637). Scaling the TPM values considering the number of proteins in which the peptide is found did again not improve performance (Wilcoxon Test, p > .05 Figure 3F).

It is also possible that a peptide occurs multiple times in a single protein sequence, e.g., in the case of repeating amino acid sequences (Luo and Nijveen, 2014). In our calculations above, such peptides were only counted once and we did not see an increase in performance when we considered duplicate peptides in a protein (Supplementary Figure S3). This is likely due to the fact that 99.6% of the peptides in the human proteome do not occur multiple times in the same protein.

Taken together, we used pepX successfully to retrieve abundance levels of peptides’ source proteins and showed that in cases where peptides can be retrieved from multiple proteins, summing up the TPM values of the encoding genes performs best in distinguishing ligands from decoy peptides.

We next investigated how well peptide abundance can be estimated using transcript-level instead of gene-level TPM values. We again used pepX to retrieve transcript-level TPM values from HPA for the peptides from the HLA Ligand Atlas and the set of random decoy peptides. We evaluated the performance of using 1) the median TPM of all transcripts, 2) using the maximum TPM among all transcripts, and 3) summing up TPM values of all transcripts. Similar to what we have observed when using gene-level TPM, summing up the TPM values of all transcripts (total TPM) significantly outperformed the median and the maximum TPM in the HLA Ligand Atlas dataset (Kruskal-Wallis test, p < .0001, Wilcoxon Test, p≤.0001, Figure 4A). In the validation datasets, the total TPM also significantly outperformed the median TPM (Kruskal-Wallis test, p < .05, Wilcoxon Test, p < .05, Figure 4B). With a mean AUC of .816 the total TPM performed slightly better than the maximum TPM with a mean AUC .814; the difference was, however, not significant (Wilcoxon Test, p > .05).

Comparing the total gene TPM and the total transcript TPM for the peptides from HLA Ligand Atlas and the set of random decoy peptides showed that transcript-level TPM values perform significantly better than gene-level TPM values (Wilcoxon Test, p≤.001, Figure 4C). On the validation datasets, with a mean AUC of .816, transcript-level total TPM performed better than the gene-level total TPM with an AUC of .811, however not significantly (Wilcoxon Test, p > .05, Figure 4D).

The Genotype-Tissue Expression Project (GTEx) is another database that provides tissue-specific gene expression data from healthy tissue samples. We wanted to compare the performance of using GTEx and HPA transcript-level expression data for estimating the abundance of the peptides from the HLA Ligand Atlas. We focused on the 18 tissue types that we could clearly match between the three datasets ( Supplementary Table S1). For this subset of peptides, using TPM values from HPA significantly outperformed using TPM values from GTEx (AUC of .812 vs. .805, Wilcoxon Test p < .01, S Supplementary Figure S4).

All expression databases provide transcript-level TPM values calculated with RSEM (Li and Dewey, 2011). The TCGA also provided TPM values calculated with Kallisto (Bray et al., 2016). We compared the performance of the two metrics using the Pyke Cancer dataset. The two metrics performed very similarly: the TPM calculated using RSEM had an AUC of .787 and the one calculated using Kallisto had an AUC of .786.

Taken together, we have shown here that, if available, transcript-level TPM data should be used to estimate peptide abundance, regardless if RSEM or Kallisto was used to calculate the TPM values. In the case of expression data of healthy tissue, HPA seems to be slightly more accurate for estimating peptide abundance of ligands eluted from healthy tissue samples.