TruNeo: an integrated pipeline improves personalized true tumor neoantigen identification

Background Neoantigen-based personal vaccines and adoptive T cell immunotherapy have shown high efficacy as a cancer treatment in clinical trials. Algorithms for the accurate prediction of neoantigens have played a pivotal role in such studies. Some existing bioinformatics methods, such as MHCflurry and NetMHCpan, identify neoantigens mainly through the prediction of peptide-MHC binding affinity. However, the predictive accuracy of immunogenicity of these methods has been shown to be low. Thus, a ranking algorithm to select highly immunogenic neoantigens of patients is needed urgently in research and clinical practice. Results We develop TruNeo, an integrated computational pipeline to identify and select highly immunogenic neoantigens based on multiple biological processes. The performance of TruNeo and other algorithms were compared based on data from published literature as well as raw data from a lung cancer patient. Recall rate of immunogenic ones among the top 10-ranked neoantigens were compared based on the published combined data set. Recall rate of TruNeo was 52.63%, which was 2.5 times higher than that predicted by MHCflurry (21.05%), and 2 times higher than NetMHCpan 4 (26.32%). Furthermore, the positive rate of top 10-ranked neoantigens for the lung cancer patient were compared, showing a 50% positive rate identified by TruNeo, which was 2.5 times higher than that predicted by MHCflurry (20%). Conclusions TruNeo, which considers multiple biological processes rather than peptide-MHC binding affinity prediction only, provides prioritization of candidate neoantigens with high immunogenicity for neoantigen-targeting personalized immunotherapies.

Next-generation sequencing data has been widely applied to predict neoantigens, and a number of bioinformatic tools have already been developed. For example, NetMHCpan [9] and MHCflurry [10] can predict the ability of mutant peptides to bind to class I or class II HLAs. These tools achieve good performance as predictors of binding affinity but poor performance as predictors of actual HLA presentation, let alone immunogenicity, or T cell response against neoantigens. Previous studies reported that less than 5% of neoantigens identified using these methods can be successfully found on the surfaces of tumor cells. The in-silico neoantigen prediction models based mainly on MHC binding affinity are limited by low predictive performance for actual immunogenicity, likely because these models considered one of the multiple steps in the neoantigen presentation processes while ignoring other steps, hence potentially underestimating the complexity of forming of true immunogenic neoantigens. The typical pipeline for the identification of immune-relevant neoantigens consists of six main steps [11].Therefore, it is necessary to add more biological factors to these prediction algorithms. We have developed an integrated pipeline called TruNeo to predict neoantigens by considering the following biological factors: peptide-MHC class I binding affinity, proteasomal C terminal cleavage, transporter associated with antigen processing (TAP) transport efficiency, expression abundance, tumor heterogeneity, clonality and HLA LOH (loss of heterozygosity).
Immunotherapies are partially limited in the number of specific antigens that can be targeted. In clinical trials [2,5,12], the top-ranked 10-20 mutations predicted by bioinformatic tools are incorporated into neoantigen-based personal vaccines. Thus patients can benefit if more immunogenic neoantigens of greater immunogenicity are included in the top-ranked predicted neoantigens. The aim of our study was to improve the positive immunogenicity rate of top-ranked predicted neoantigens, so we compared the predictive performance of TruNeo and other algorithms according to this evaluation criterion. We compared several predictive methods by counting the number of immunogenic neoantigens in the top-ranked 5, 10 or 20 mutations of 13 patients from published data. We also compared the top 10 personalized neoantigens in a lung cancer patient as predicted by TruNeo and MHCflurry. We show that the integrated TruNeo pipeline model improves personalized true tumor neoantigen identification.

TruNeo pipeline
The TruNeo pipeline required two data types as input; raw DNA sequencing FASTQ files from the paired tumor and normal samples, and RNA-seq data from a tumor sample. In the first step annotated somatic mutation information, HLA genotype and gene expression information were prepared. Then, the candidate neoantigens were predicted based on peptide-MHC binding affinity. In the third step, candidate neoantigens were scored by integrating information from multiple neoantigen presentation processes. Lastly, high confidence neoantigens were filtered and output.

Generation of annotated mutation information, HLA genotype and gene expression information
Raw sequencing data of DNA from paired tumor and normal samples were aligned to a reference genome for generating bam files, and then somatic single nucleotide variants (SNVs) and InDels were identified and annotated. Normal bam was used for HLA genotyping. RNA FASTQ file from tumor samples were used for fusion identification and gene expression quantification.
The GATK Somatic InDel Detector (v3.3, GATK IndelRealigner) was used to identify somatic insertions and deletions (InDels) with the default parameters. InDels with high confidence were filtered using the following steps: (1) local realignment was performed with combined normal and tumor BAM files for each predicted somatic InDel; (2) depth at mutation position ≥ 10×, VAF of InDels > 10% in tumor and < 2% in normal.
Finally, all SNVs and InDels were annotated using an in-house annotation software based on snpeff [16].

Genotyping and loss of heterozygosity in HLA class I genes
The HLA genotype was identified with the combined use of Polysover (v1.0) [17] and BWA-HLA (v1.3). If Polysover identified the same genotype in both tumor and normal samples, Polysover's result was taken as the HLA genotype. If not, the result of BWA-HLA was checked. If BWA-HLA identified the same genotype across tumor and normal samples, the BWA-HLA result was used. If not, both the result of Polysover and BWA-HLA were considered. If Polysover and BWA-HLA identified the same genotype in the normal sample, the normal HLA genotype result was used. If Polysover and BWA-HLA were not in agreement, Polysover's result in the normal sample was taken as the HLA genotype and marked as low confidence.
The tumor and matched normal sequences were mapped to the HLA reference, and an HLA-LOH event was reported if paired t test of HLA allelic imbalance was significant (p value < 0.0002). The allelic imbalance is tested using the ratio of log2 (Tumor unique reads/Normal unique reads) [18].

Fusion identification and gene expression quantification from RNA sequencing data
The raw RNA-seq data were processed with STAR v2.5.3a [19], and the gene expression level was estimated as the transcripts per million (TPM) via RSEM v1.3.0 [20]. In addition, RNA-based gene fusion was detected by STAR fusion [21], which provided another source of neoantigens.

Prediction of candidate neoantigens based on peptide-MHC binding affinity
21-mer polypeptides centred on mutated residues were scanned to identify candidate peptides binding to class I HLAs, such as peptide sequences surrounding mutated amino acids resulting from missense mutations and frame-shift or non-frame-shift InDels. The binding affinity of 8-11-mer peptides for class I HLAs was predicted using the NetM-HCPan 3.0 [22] binding algorithm. Epitopes were filtered if the following conditions were met: (1) mutations were not expressed according to RNA-seq data (mutations with mutant allele reads ≥ 1 in RNA sequencing data [23] were confirmed as expressed); (2) the sequence was homologous to self; (3) the half-maximum inhibitory concentration (IC50) according to NetMHCPan 3.0 was larger than 500 nM.

Scoring of candidate neoantigens by integrating information from multiple neoantigen presentation processes
We first combined biological processes including MHC binding, proteasomal cleavage efficiency and TAP transport efficiency, and then integrated variant allele frequency, expression abundance and type of neoantigen. Thus, the final score for each neoantigen was calculated as follows: The proteasomal cleavage efficiency was predicted by netChop [24], and the TAP transport efficiency was predicted by netCTLpan [25].
Expression score was identified by TPM, and normalised based on ranking status.
The deep learning-based model was trained on a large mass spectrometry HLA peptide data set from various human tumors using a neural network structure as follows: (1) Training dataset: 8-11mer peptide paired with HLA genotype from published data on 74 patients by mass spectrometry [26] were selected as positive neoantigen. Some random 8-11mer peptides from reference proteome (Uniprot protein database) paired with HLA genotype were used as negative neoantigen; (2) Data process: Peptides were vectorized using a one-hot encoding scheme; Embedding HLA genotypes to vectors; (3) Model architecture: A 256, 74 neurons fully connected neural networks, using relu and sigmoid as activation function. A 74 long embedding vector of HLA type to control the output from fully connected layer from each HLA type of a patient; (4) Training: Split 10% data as validation set. Use binary-crossentropy as loss function to optimize model until the loss function value of validation set stop decreasing.

Validation cohort from published studies
As raw sequencing data was not available, single-nucleotide variants of 13 patients from published articles were collected. Enzyme-linked immune-spot (Elispot) assays were used to mark SNVs as immunogenic. Finally, 1599 assayed single-nucleotide variants from 13 patients were collected from published studies, 19 of which were immunogenic. Neoantigen prediction using the TruNeo pipeline, NetMHCpan [20], MHCflurry, PSSMHCpan, and DeepHLA started from input of annotated mutation list, HLA genotype, and gene-level TPM. Prediction using EDGE was collected from a published article [33].
Recall rate (true positives/19) was used as evaluation criterion of predictive performance.

Whole exome DNA and RNA sequencing of a lung cancer patient
A 68-year-old patient (patient 01) with squamous cell carcinoma of the lung (SCLC) was enrolled in the study at Shanghai Tenth People's Hospital in 2018. This study was approved by the Shanghai Tenth People's Hospital ethics committee and the patient provided written informed consent. A tumor biopsy and peripheral blood samples were collected for whole exome sequencing and transcriptome sequencing to identify mutations and potential neoantigens. PBMCs were also used to conduct the Elispot assay.
DNA and RNA from fresh tumor that were isolated pre-treatment and DNA from paired blood samples were extracted, purified, and hybridized using the Agilent Sure-Select Target Enrichment System kit (Qiagen, USA) according to the manufacturer's instructions, and paired-end multiplex sequencing of samples was performed on the Illumina Novaseq 6000 sequencing platform. The average sequencing depth was 258× in tumor tissue and 126× in paired peripheral blood.

Evaluating the performance of neoantigen prediction methods by Elispot assays of PBMCs from a cancer patient
To assess whether the multidimensional pipeline performed better than the single-factor methods, we validated the immune responses to neoantigen candidates identified with each method. The peripheral blood cells were obtained from patient 01, and somatic variants identified as described above. We chose two methods to predict neoantigens, the multidimensional TruNeo and open-sourced deep-learning-based MHCflurry methods [10]. Afterwards, the immunogenicity of the 10 top-ranked candidate peptides identified by each software were validated by the Elispot assay [27] DCs were cultured as previously described [28]. CD8+T cells sorted from PBMCs were stimulated in 24-well cell culture plates with autologous DCs pulsed with individual neoantigen peptides (10 μg/mL) and IL-7 (10 ng/mL; PeproTech). On day 3, IL-2 (5 μg/mL; PeproTech) was added. Half of the medium was changed, and the addition of cytokines was performed every 3 days, as described previously [2]. After 10 days, the IFN-γ response of the prestimulated T-cells was tested against neoantigens by Elispot assays with a Human IFN-γ Elispot kit (MabTech). The Elispot plates were washed five times with PBS. Prestimulated CD8+T cells and DCs pulsed with neoantigen peptides were added to individual wells of the plates and incubated at 37 °C for 18-24 h in the presence of 5% CO 2 . The plates were washed five times with PBS and incubated for 2 h with 100 μL/well anti-human IFN-γ (7-B6-ALP) at 37 °C. Then, the plates were washed five times with PBS and incubated with BCIP/NBT-plus substrate at room temperature. The resulting spots were counted using a computer-assisted Elispot image analyser (Biosys Bioreader 4000), and custom software was designed to detect spots using predetermined criteria based on size, shape, and colorimetric density. The measurement of the spot-size distribution is a built-in function of the software. According to the established guidelines [27], a positive response was defined when the mean of the antigen-stimulated replicates was greater than or equal to ten spots per well, and the mean of the antigen-stimulated replicates was greater than two times the mean of the replicates of the negative control wells.
True positive rate (true positives/(true positives + false positives)) was used as the evaluation criterion.

T cell receptor (TCR) sequencing and neoantigen-specific TCR clone analysis
The methods of stimulation and expansion of neoantigen-specific T cells were the same as the methods used for the preparation of prestimulated CD8+T cells. CD8+T cells stimulated by DCs without pulsed peptides or CD8+T cells were used as negative controls. Cells were harvested on day 10 and washed twice with PBS. Cultured T-cell pellets were flash-frozen in liquid nitrogen and stored at − 196 °C.
RNA was extracted from flash-frozen peptide-stimulated T cells using TRIzol reagent (Invitrogen. USA). The CDR3 region of the TCRβ chain was amplified by using the iRepertoire multiplex primer set (iRepertoire, Inc), and sequencing was performed using the Illumina 4000 system (Illumina Inc.). Bioinformatic analysis of productive clones was performed to identify antigen-specific expansion using the following criteria: (1) significant expansion (Fisher's exact test with a Benjamini-Hochberg FDR, p < 0.05) compared to that of T cells cultured without peptide, (2) no significant expansion of the relevant clone in any other peptide-stimulated culture, and (3) an odds ratio > 1 (default value).

Overview of the TruNeo pipeline from the processing of next-generation sequencing data to the identification of candidate neoantigens
In this study, we built an integrated pipeline to predict neoantigens called TruNeo (Fig. 1). The main purpose of TruNeo is to score candidate neoantigen and prioritize 10-20 top-ranked neoantigens for personalized neoantigen-based immunotherapy. The pipeline began with FASTQ data from paired tumor and normal DNA, along with tumor RNA expression data. The first step was to prepare annotated mutation, fusion list, the alleles of human leukocyte antigen (HLA), as well as RNA expression quantification. The second step was to predict candidate neoantigens according to MHC binding affinity (IC50 threshold < 500 nm). All possible 8-11-mer amino acid fragments were derived from SNVs, InDels and fusions, MHC binding affinity was predicted using NetMHCpan. The third step was to integrate a variety of biological factors contributing to immunogenicity for ranking. These factors involved proteasomal cleavage, TAP-mediated peptide transportation, the anchoring residue of HLA, homologous sequences, expression abundance, variant allele frequency and HLA LOH. Candidate peptides were further filtered for high confidence neoantigens based on high affinity of HLA binding, high expression abundance and high VAF. The top-ranked epitopes were considered to be more immunogenic as candidate neoantigens (see "Methods" section).

The performance of TruNeo and other algorithms on published data
We collected 1599 non-redundant SNVs from 13 patients through published articles [29][30][31], among which 19 neoantigens were identified as immunogenic true positives that caused a T-cell response using IFN-γ ELISPOT assays. Predictive performance was compared among of TruNeo and other 5 algorithms, including NetMHCpan 4.0, MHCflurry, PSSMHCpan [32], DeepHLA [33] and EDGE [26]. NetMHCpan, MHCflurry, and PSSMHCpan are machine learning models trained on peptide-HLA binding dataset, which identify neoantigen based only on peptide-HLA binding affinity. DeepHLA and EDGE are deep learning models trained on HLA specific tumor mass spectrometry data, which can predict the presented neoantigens or immunogenic neoantigens. To compare the 6 algorithms, the recall rate was calculated in the top ranked 20, 10 and 5 neoantigens provided by each algorithm. Our results (Fig. 2, Additional file 1: Tables S1) showed that TruNeo could rank the immunogenic true positive neoantigens better than HLA binding-based prediction methods (p = 0.098, one-sided paired Wilcoxon rank sum test). For example, TruNeo was able to rank 13 immunogenic neoantigens in the top 20,

Performance of TruNeo and MHCflurry in a SCLC patient
To further estimate the performance of TruNeo, we predicted the candidate neoantigens in a real case: a patient with advanced squamous cell carcinoma (patient 01). Fresh tumor tissue and blood were collected pre-treatment followed by next generation sequencing. TruNeo were used to call somatic mutations, fusions, HLA genotyping and expression quantification. 451 somatic mutations were identified, including 313 nonsilent mutations (Additional file 2: Table S2) (Fig. 3). The non-silent mutations consisted of 297 missense mutations, 2 in-frame mutations, and 14 frameshift mutations. The class-I HLAs of patient 01 were typed as HLA-A*11:01, HLA-A*02:10, HLA-B*40:01, HLA-B*40:01, HLA-C*08:01, and HLA-C*07:02, which were double-checked by assaying tumor and blood samples. These results were further used in neoantigen prediction and ranking through TruNeo and MHCflurry. 254 short candidate peptides (Additional file 2: Table S2) were identified as candidate neoantigens by TruNeo and 395 by MHCflurry. Then, candidate neoantigens were ranked by TruNeo and MHCflurry separately.  The top 10 ranked neoantigens selected by TruNeo and MHCflurry separately are shown in Table 1. There was one neoantigen identified by both methods. In total, 19 unique neoantigens peptides were synthetized by GenScript Corporation (Nanjing, China), and then validated using Elispot assays. Five immunogenic neoantigens (#1, #3, #4, #6, and #8) of the top 10 predicted neoantigens of patient 01 identified by TruNeo were showed immunogenic activity, and 2 (#5 and #6) of the top 10 predicted neoantigens identified by MHCflurry were showed immunogenic activity (Fig. 4). The number of true positive neoantigens among the top 10 predicted by TruNeo was 2.5 times higher than that predicted by MHCflurry (50% vs 20%).
If a neoantigen is immunogenic, there should be specific corresponding TCR clones among PBMCs. To verify the immunogenicity of the neoantigen, we chose the top 1 ranked neoantigen validated by Elispot assays and found that it had two specific 5 immunogenic neoanƟgens of top 10 neoepitopes ranking by TruNeo 2 immunogenic neoanƟgens of top 10 neoepitopes ranking by MHCflurry   Fig. 4 Immunogenic neoantigens of patient 01 predicted by TruNeo and MHCflurry validated by Elispot assays. We have predicted the top10 MHC class I neoantigens for patient01 by TruNeon and MHCflurry. According the rank of the neoantigens were labeled as #1-#10, the immunogenicity of the neoantigen peptides predicted by two models have been validated by the Elispot assay. #1, #3, #4, #6, #8 of the top 10 neoantigen predicted by TruNeo were Elispot positive; #5, #6 of the top 10 neoantigen predicted by MHCflurry were Elispot positive corresponding TCR clones (Table 2), which further supports the immunogenicity of this neoantigen from the perspective of the interaction between the peptide and the TCR.

Discussion
Recent studies have shown that neoantigen peptides predicted by current bioinformatic tools such as NetMHCpan or MHCflurry were found on the surface of cells was lower than 5% [34,35], likely because the training data capture information about only one of multiple steps in the HLA class I processing pathway [25]. The pathway from DNA mutations to neoantigens is a complex biological process and the typical pipeline of neoantigen prediction consists of 6 main steps. These steps include identification of somatic mutation, transcription into mutated mRNA, proteasomal processing of the mutated protein, transport of TAP-mediated peptide into the ER lumen, binding of this peptide to the MHC I protein on the endoplasmic reticulum. Finally, the neoepitope MHC complex can be transported to the cell membrane for recognition by TCRs. While current neoantigen identification algorithms rely primarily on the prediction of peptide-HLA binding affinity, it is not sufficient for neoantigen prediction. Secondly current algorithms are mainly machine learning models trained on large in vitro peptide-HLA binding dataset. They have excellent performance as predictors of peptide-HLA binding affinity, but poor performance as predictors of actual neoantigen presentation. Another issue to be considered is tumor cell heterogeneity. Neoantigens may not be expressed in all tumor cells, so the tumor cell fraction and expression abundance of neoantigen should be considered. expression abundance of neoantigen have been proved have relation to the positive rate of neoantigen validation, while current algorithms don't take into account. Other features, such as mutation due to anchor residues [36], self-proteome homologs, and diversity of HLA molecules have also been shown to be associated with immunogenicity [35]. Thus, the use of a single predictor is less accurate when prioritizing potential neoantigens. We have demonstrated that TruNeo, a pipeline which considers multiple biological factors, can predict and rank high-quality actionable neoantigens from whole-exome and transcriptome data. TruNeo can predict neoantigens derived not only from point mutations but also from insertions, deletions and fusion genes. Important biological steps, including proteasomal processing and TAP-mediated peptide transport, HLAbinding affinity, presence of homologous sequences, clonality, and gene expression status are considered during annotation and ranking to select the top neoantigens that are most suitable for vaccine development or adoptive cell therapy.
Consideration of features likely contributes to the improved predictive performance of TruNeo compared to existing methods, among which gene expression might be the most critical factor. Thus, the performance of existing methods could be significantly improved by including a gene expression threshold. By raising the minimal TPM (Transcripts Per Million) threshold from 0 to 2, the proportion of CD8-recognized mutations in the top 10 neoantigens was increased from 21 to 42% by MHCflurry, 26-32% by netMHCpan 4.0, 37-47% by DeepHLA, and 37-43% by PSSMHCpan. This result was also demonstrated in previous studies [26,33], suggesting that expression greatly contributes to high-confidence neoantigen identification.
For class I HLA antigens, we analysed and compared the published data set of the Elispot validation cohort and found that 52.6% of the confirmed positive neoantigens were ranked in the top 10 by TruNeo, which outperformed MHCflurry, NetMHCpan, PSSMHCpan and DeepHLA. For a single case, we found that the positive rate of the top 10 ranked by TruNeo was 50%, compared to 20% using MHCflurry. These results show that the top-ranked neoantigens identified by TruNeo have an increased true positive rate compared with those ranked using standard HLA binding affinity prediction methods.
We have not only developed a neoantigen prediction pipeline but also an experimental platform for neoantigen validation (Elispot assay and TCR-seq) (Fig. 5). TCRseq can provide the TCR CDR3 sequences of immunogenic neoantigens validated by Elispot assays, which are useful in dynamically monitoring the immune function status of patients treated with immunotherapy.
There are also some limitations to our study. As described previously [7], the immunogenicity of 68% of neoantigens was validated post vaccination, and some naive T cells had transformed into neoantigen-specific T cells after vaccination. However, in our study, the PBMCs used for the Elispot assay were collected prior to vaccination, which might influence the true positive rate during validation. Another limitation is that we did not verify HLA class II neoantigens. Two reported clinical studies [2,7] found that CD4+T cells induced by MHC II comprise the main response by T cells. However, antigenic peptide binding to MHC class II is affected by the long length, poor conservation, and multiple motifs of these antigens, so it is more difficult to predict than binding of MHC class I to antigenic peptides. Work is ongoing to optimize the TruNeo pipeline for prediction of neoantigens presented by HLA class II. A final limitation is that we use Fig. 5 Overview of neoantigen prediction and identification pipeline. To identify the predicted neoantigens, T cells are sorted from PMBC and stimulated by DC pulsed with peptides. Next, CDR3 clones are analysed through TCR Seq, along with applying Elispot assays to confirm immunogenicity a single real patient in comparing TruNeo's performance to other tools using raw data. Future work should focus on recruiting additional patients to repeat validation experiments for comparing the predictive performance of current tools. This can provide stronger evidence of the relative performance of each method. However, the validation of neoantigen immunogenicity is costly. Many neoantigen peptides must be selected and synthetized for algorithm comparison. Large amounts of peripheral blood is needed for PBMC isolation, or fresh tumor samples are needed to culture TIL for immunogenic validation. Seeking an appropriate cancer patient and blood donor is still a costly and time-consuming work. We hope to improve this aspect with high-throughput, low original input validation platforms in the future.
TruNeo was also compared with other deep learning-based algorithms including EDGE and DeepHLA. We found that EDGE ranked best among the 6 methods. MSIn-trinsicEC, another deep learning method, outperformed standard methods by twofold as described [37] which was similar to TruNeo. As described in MSIntrinsicEC, only 16 single-HLA-expressing cell lines were collected for training, which means that MSIn-trinsicEC works quite well only with 16 HLA alleles. On the contrary, TruNeo was not limited in terms of HLA alleles. A deep learning method with mass spectrometry data, instead of in vitro HLA-peptide binding affinity data, can deliver neoantigen probabilities without the tediously biological features assessment processes. Deep learning models can be applied to improve TruNeo in the aspect of filtering and ranking when the neoantigen experimental validation data set become large enough. Moreover, neoantigen validated data from spectrometry, Elispot and multiplexed tetramer binding assays are helpful in improving the accuracy of the algorithm.
Cancer immunotherapies which target neoantigens are of growing interest and are in the early stages of human trials, but methods to identify neoantigens require invasive or difficult-to-obtain clinical specimens, require the screening of hundreds or thousands of synthetic peptides or tandem minigenes, or may only be relevant to specific HLA alleles. We have created a neoantigen identification and ranking pipeline that considers multiple factors. Our model increases the immunogenicity rate of the top 10 predicted neoantigens to 50%. We hope that in the future, with the accumulation of positive neoantigen databases that can be used as training sets, developing a prediction method for neoantigens will help to optimize the composition of personalized cancer vaccines and mass spectrometry T cell therapy with high precision and will speed up vaccine and ACT design to meet growing clinical needs.