Long-read sequencing reveals widespread novel splicing and neojunction-derived neoantigens in nasopharyngeal carcinoma

LR-seq unveils a vast repertoire of novel transcripts and extensive splicing diversity in the NPC transcriptome

To comprehensively profile the AS landscape in NPC, we performed long-read RNA-seq on a cohort of 18 samples, including both nonmalignant and malignant samples. This cohort included two NPC cell lines, 10 patient-derived xenografts (PDXs), two NPC tumor tissues, and four immortalized nasopharyngeal epithelial (NP) cell lines. All samples were analyzed using SMRT circular consensus sequencing (CCS) on the PacBio Sequel II platform. On average, about 783,000 CCS reads were obtained per library, yielding about 350,000 clustered isoforms after downstream processing (Table 1; Supplemental Table S1). For each sample, the numbers of genes and full-length transcripts identified in LR-seq after collapsing were reported (Supplemental Fig. S1B). The observed differences across samples were expected and likely reflected variation in RNA quality and biological heterogeneity. After merging transcripts across samples and collapsing redundant isoforms, we identified 51,115 high-quality full-length transcripts with a median transcript length of ∼3.0 kb (Fig. 1A,C), constructing a comprehensive LR-seq-derived transcriptome for NPC.

Figure 1.

Profiling of full-length transcripts in nasopharyngeal carcinoma (NPC) by long-read sequencing (LR-seq). (A) Pie chart illustrating the distribution of full-length transcripts across different subcategories. A total of 51,115 full-length transcripts were identified across 18 LR-seq samples, with the numbers and proportions of different transcript subcategories highlighted in different colors. (B) Bar plots showing LR-seq isoforms detected in individual NPC tumor samples and immortalized NP cell lines, colored by the subcategories defined in A and grouped by sample origin. (C) Histograms comparing length distributions of long-read transcripts between different subcategories. (D) Bar plot showing the number of transcripts identified per gene by LR-seq. (E) Box plot comparing the proportions of long-read transcripts supported by short-read Illumina RNA-seq among different subcategories in matched samples.

Open in new tabDownload PowerPointLink to figure

Transcripts were categorized using SQANTI3 based on their structural alignment to the GENCODE v38 reference transcriptome. This classification revealed 24.10% (12,319) as full splice match (FSM), perfectly matching a known reference transcript; 26.10% (13,341) as incomplete splice match (ISM), aligning to a known transcript but lacking complete 5′ or 3′ ends; 22.65% (11,578) as novel in catalog (NIC), containing known splice sites in novel combinations; and 21.73% (11,109) as novel not in catalog (NNC), featuring at least one previously unannotated splice site. The remaining 2768 transcripts were classified as antisense, genic, fusion, or intergenic transcripts (Fig. 1A). Novel transcripts (NIC/NNC) constituted a substantial portion of the transcriptome, representing between 17% and 36% of all transcripts per sample (Fig. 1B), underscoring the substantial unannotated transcriptional complexity in NPC. Collectively, these transcripts were mapped to 11,807 annotated genes and 2259 novel genomic loci. Notably, 56.66% (6690) of the genes in the LR-seq-derived transcriptome expressed two or more transcripts, indicating widespread AS (Fig. 1D).

To assess the reliability of the LR-seq-derived transcriptome, we aligned short-read RNA-seq data from matched samples to our LR-seq-derived transcriptome. This validation revealed a clear hierarchy of support (transcripts per million [TPM] ≥1): FSM transcripts showed the highest validation rate (∼68%–80%); NIC showed intermediate validation rates (∼40%–60%); NNC and ISM transcripts showed lower support (∼30%–40%); and antisense, genic, fusion, or intergenic transcripts exhibited the lowest validation rates (∼10%–25%) (Fig. 1E). We next evaluated the distribution of isoforms across samples. We quantified the number of samples supporting each isoform using two criteria: long-read presence/absence calls (binary “0/1” per sample) and short-read RNA-seq expression (counting samples with TPM above the multiple thresholds). FSM transcripts exhibited low heterogeneity, with many isoforms detected across multiple samples. Conversely, novel transcripts exhibited high sample-to-sample variability, with the majority detected in only a few samples, highlighting their greater heterogeneity relative to well-annotated FSM isoforms (Supplemental Fig. S2A,B).

Comprehensive characterization of novel transcripts with quality evaluation reveals their structural diversity, functional relevance, and potential roles in NPC pathogenesis

Next, we performed quality assessments of the identified transcripts, focusing on coding potential, susceptibility to nonsense-mediated mRNA decay (NMD), and splice site characteristics. We observed that the vast majority of FSM and novel transcripts were predicted to be protein coding (94%–97%), whereas antisense, genic, fusion, or intergenic transcripts largely lacked coding potential. NMD prediction revealed that NNC, antisense, genic, fusion, or intergenic transcripts were more likely to contain premature termination codons (PTCs) and noncanonical splice sites (Fig. 2A), consistent with increased susceptibility to RNA surveillance and degradation (Supek et al. 2021). Investigating the mechanisms triggering NMD indicated that the use of novel splice sites and intron retention were the primary causes, often generating PTCs that mark transcripts for degradation (Supplemental Fig. S2C). For example, two novel ANO9 transcripts (TCONS_00009816 and TCONS_00009828) were predicted to undergo NMD owing to a novel splice site and intron retention, respectively, both introducing PTCs (Supplemental Fig. S2D). Furthermore, we observed a positive correlation between the number of exons per gene and the number of novel transcripts produced by that gene, suggesting that genes with higher exon complexity have a greater potential for generating novel splice variants (Supplemental Fig. S2E).

Figure 2.

Characterization of full-length transcripts identified in LR-seq. (A) Bar plot comparing the coding potential, predicted NMD status, and splice site canonicality between transcript subcategories. (B) Box plots comparing isoform length, CDS length, and junction coverage across the FSM, NIC, and NNC transcript categories. (C) Dot plot showing pathways significantly enriched among genes with novel transcripts detected by LR-seq (FDR < 0.05), based on the MSigDB C2, C5, and Hallmark collections. Enrichment remained significant after controlling for gene expression using expression-matched permutation analysis. (D) Dot plot illustrating the numbers of novel LR-seq transcripts compared with annotated GENCODE isoforms detected in selected oncogenes (left) and tumor suppressors (right). (E) Pie chart showing the proportion of novel transcripts sharing the same translation start site as annotated protein isoforms. (NS) not significant, (*) P < 0.05, (**) P < 0.01, (***) P < 0.001, (****) P < 0.0001.

Open in new tabDownload PowerPointLink to figure

Comparison with reference transcripts suggested that NNC isoforms often arose from unannotated transcriptional start sites (TSSs) and transcriptional termination sites (TTSs) located >100 bp from known transcript boundaries (Supplemental Fig. S2F). In addition to RNA-seq, we integrated multiple orthogonal data sets to validate novel isoforms identified by LR-seq, including cap analysis of gene expression (CAGE) and 3′-seq peaks from the poly(A) site database. The 5′ ends of many novel isoforms overlapped with FANTOM5 CAGE-defined TSSs, and their 3′ ends were validated by bona fide TTSs mapped by 3′-seq poly(A) databases (Supplemental Fig. S2G). Comparative analysis revealed distinct structural features distinguishing FSM transcripts from novel transcripts. NNC isoforms were typically shorter in total length compared with FSM and NIC transcripts, yet novel transcripts often possessed longer coding sequences (CDSs) compared with FSM isoforms. Furthermore, splice junctions in FSM transcripts showed higher read coverage, whereas novel isoforms, particularly those with novel splice sites, had lower read support (Fig. 2B).

To evaluate the biological relevance of novel transcripts, we performed pathway enrichment analysis. Genes associated with novel transcripts were significantly enriched in pathways related to RNA processing, nuclear export, and intracellular transport, based on a matched background gene set (Fig. 2C). Notably, several established oncogenes and tumor suppressors, including FGFR1, which can be activated by LMP1 to promote glycolysis and invasion in NPC (Lo et al. 2021), DDX5, implicated in NAT10-driven immune suppression (Xie et al. 2025), and MST1R, which is linked to early-onset NPC and innate immunity (Dai et al. 2016; Cazes et al. 2022), exhibited substantial numbers of novel isoforms, suggesting a crucial role in NPC pathogenesis (Fig. 2D). Most novel transcripts (67.8%) shared the same translation start sites as known protein isoforms, indicating that the associated splicing events mainly affected internal exons. In contrast, 32.2% potentially used alternative translation start sites, possibly associated with alternative promoter usage or first exon selection, which may alter the 5′ UTR or the exon harboring the canonical start codon (Fig. 2E).

In addition to novel host transcripts identified in NPC samples, we also detected previously unannotated EBV transcripts, including transcripts with an alternative first exon in LMP2 and exon skipping of exons 2 and 3 in RPMS1 (Supplemental Fig. S2H). To assess the differences between sequencing approaches, we compared the transcriptome derived from LR-seq with the transcriptome assembled from short-read data using StringTie. The StringTie-assembled transcriptome comprised 239,952 transcripts, far exceeding the number identified by LR-seq. However, only 14,198 transcripts were commonly identified by both methods, most of which were FSMs. This discrepancy underscored the substantial variation in transcript identification between short-read and long-read sequencing technologies, highlighting the limitations of short-read-based assembly in reconstructing full-length isoforms (Supplemental Fig. S2I).

Moreover, a key advantage of LR-seq lies in its ability to accurately identify RTTs by capturing full-length transcript structures. This enabled precise delineation of fusion transcripts spanning adjacent genes. These RTTs and their encoded products have emerged as potential diagnostic biomarkers and therapeutic targets in various cancers (Alpert et al. 2020; Hu et al. 2022). In our study, we identified and experimentally validated four such RTTs, including SFT2D2-TBX19, FOXRED2-TXN2, ELAC1-SMAD4, and CLDN15-FIS1, using RT-PCR with specific primers targeting the fusion junctions (Supplemental Fig. S3A,B).

Next, we annotated seven types of alternative splicing events (ASEs) using SUPPA2, including the skipping exons (SEs), alternative 3′ splice sites (A3s), alternative 5′ splice sites (A5s), alternative first exon (AF), alternative last exon (AL), intron retention (RI), and mutually exclusive exons (MX) categories (Fig. 3A), in the GENCODE v38 transcriptome and its merged transcriptome incorporating novel transcripts identified by LR-seq. Notably, incorporating novel isoforms increased the number of detected ASEs, particularly AF and SE, whereas the counts of AL and MX remained largely unchanged (Fig. 3B). Differential ASE analysis between NPC tumor samples and immortalized NP cell lines identified 102 differentially spliced events (P < 0.05, |ΔPSI| > 0.1) (Supplemental Table S2), with AF (43, 42.16%) and SE (25, 24.51%) being the most prevalent (Fig. 3C). A prime example was TLN1, which exhibited a tumor-specific switch involving inclusion of a cancer-associated cassette exon (17b). This isoform, which enhances vinculin binding and cell motility (Gallego-Paez et al. 2023), was significantly upregulated in NPC samples (Fig. 3C,D). We validated this isoform switch by RT-qPCR, confirming that the novel isoform (TCONS_00049374) accounted for a significantly larger fraction of total TLN1 expression in tumors than the canonical isoforms (ENST00000314888.10) (Fig. 3E).

Figure 3.

Novel transcripts contribute to isoform-level diversity and differential alternative splicing events (ASEs). (A) Schematic illustration of the seven types of ASEs: alternative 5′ splice site (A5), alternative 3′ splice site (A3), alternative first exon (AF), alternative last exon (AL), mutually exclusive exons (MX), retained intron (RI), and skipping exon (SE). (B) Bar plot comparing ASEs before and after incorporating novel transcripts into GENCODE v38 transcriptome. (C) Dot plot showing differential ASEs between NPC tumor samples and immortalized NP cell lines, identified by SUPPA2 with P < 0.05 and |ΔPSI| > 0.1. (D) Schematic illustrations of an example of a differential ASE involving the inclusion of a cassette exon (17b) in TLN1 and the resulting transcripts. (E) Bar plot comparing differential isoform usage of TLN1 (left), with and without exon 17b, between NPC samples and immortalized NP cell lines, validated by RT-qPCR (right). Each RT-qPCR experiment was repeated three times with independent RT preparations, with the standard deviation of each triplicate indicated by the error bars.

Open in new tabDownload PowerPointLink to figure

Long-read transcriptomics identifies candidate transmembrane therapeutic targets and tumor-specific transcripts

To identify transcripts enriched in NPC tumors, we quantified expression by mapping matched short-read RNA-seq data to our LR-seq-derived transcriptome. Analysis of transcript expression density curves and distribution patterns revealed that most FSMs exhibited relatively higher expression levels compared with novel transcripts (Fig. 4A,B). Using differential transcript expression (DTE) analysis between NPC tumor samples and immortalized NP cell lines, we identified 1812 significantly dysregulated transcripts (|log₂FC| > 2, BH-adjusted P < 0.05). This included 460 upregulated and 325 downregulated FSM transcripts, alongside 241 upregulated and 117 downregulated NIC transcripts and 211 upregulated and 86 downregulated NNC transcripts (Supplemental Table S3). Highly expressed alternatively spliced isoforms in NPC included known isoforms of CSF2RB, LUM, and COL1A2, as well as novel isoforms of PRRX1, LGALS9, and CD74 (Fig. 4C). Notably, overexpression of CSF2RB and PRRX1 contributes to promoting oncogenic activity and maintaining epithelial–mesenchymal transition (EMT) in human cancers (Du et al. 2021; Charlet et al. 2022). High expression of LGALS9 (Klibi et al. 2009; Kam et al. 2025) and CD74 (Chow et al. 2022) has been reported to play important roles in promoting immune evasion and remodeling the tumor microenvironment in NPC. Although the well-known transcripts were used in previous functional studies, the identification of alternative and novel transcripts, especially the novel LGALS9 and CD74 transcripts, in NPC indicated the potential for new oncogenic activities of these transcripts in NPC tumorigenesis. Because of the potential functional difference between novel and known transcripts, further studies on the NPC-specific transcripts are needed to clarify their key oncogenic activities in tumorigenesis.

Figure 4.

Identification and validation of tumor-specific transcripts. (A) Density curves comparing the distribution of transcript expression levels across different transcript subcategories. (B) Bar plot showing the numbers of FSM, NIC, and NNC transcripts across differential median expression thresholds. (C) Volcano plot illustrating differential transcript expression analysis between NPC tumor samples and immortalized NP cell lines (BH-adjusted P < 0.05 and |log₂FC| > 2). (D) Ring charts showing the classification of transcripts according to their expression specificity. (E) Heatmap showing the expression of a subset of predicted tumor-specific novel transcripts identified in LR-seq. (F) Bar plots depicting RT-qPCR validation results of transcripts in E. RT-qPCR assays were carried out in three independent RT replicates, with error bars representing the standard deviation.

Open in new tabDownload PowerPointLink to figure

We next prioritized transcripts predicted to encode plasma membrane transmembrane proteins because these proteins are stably expressed at the cell surface, present extracellularly accessible domains, and therefore represent attractive targets for therapeutic antibodies (Hu et al. 2021; Bardia et al. 2024; Neri et al. 2024; Zaidi et al. 2024; Boixareu et al. 2025). Based on predictions of transmembrane topology and plasma membrane localization, 1333 known transcripts from 859 genes were predicted to encode cell surface transmembrane proteins. Additionally, 2159 novel transcripts from 595 genes were predicted to encode cell surface proteins containing transmembrane domains (Supplemental Fig. S4A; Supplemental Table S4). Notably, we identified 294 novel transcripts from 95 genes that were predicted to gain cell surface localization with transmembrane helices, even though these genes were absent from the TCSA reference surfaceome set (Hu et al. 2021). In contrast, we identified 1378 novel transcripts from 360 genes in the curated TCSA reference surfaceome whose encoded protein products were predicted to lose cell surface localization relative to the annotated protein products of the same genes (Supplemental Table S5). DTE analysis revealed significant upregulation of several cell membrane–associated transcripts in NPC tumor samples (e.g., IGSF9, CXADR, TMPRSS2, INSR) (Supplemental Fig. S4B; Supplemental Table S3). Survival analysis of an NPC public data set (obtained from the NCBI Gene Expression Omnibus [GEO; https://www.ncbi.nlm.nih.gov/geo/] under accession number GSE102349) indicated that high expression of these cell membrane–associated transcripts correlated with poorer progression-free survival (PFS), highlighting their potential clinical relevance (Supplemental Fig. S4C).

Among the predicted plasma membrane transmembrane isoforms, predominant expression of the NTRK2-T1 protein isoform was reported in our recent publication (Wang et al. 2025). The expression of the NTRK2-T1 transcript was also confirmed by short-read RNA sequencing, and the corresponding protein expression was confirmed by western blotting in NPC. By IHC, we demonstrated the membrane staining of NTRK2-T1 protein isoform in the NPC xenografts and primary tumors (Supplemental Fig. S5A–C). To further validate the newly identified transmembrane CSF2RB protein isoform TCONS_00033496, we examined its protein expression and subcellular localization in vitro. Western blot analysis confirmed expression of FLAG-tagged TCONS_00033496 in both HEK293T and NP69 cells after transfection. Immunofluorescence staining showed that the FLAG signal was predominantly localized at the cell periphery and colocalized with the plasma membrane marker WGA in both cell lines, whereas no specific anti-FLAG signal was detected in the corresponding negative control cells. These findings indicated that TCONS_00033496 encoded a novel CSF2RB protein isoform that was properly translated and targeted to the plasma membrane (Supplemental Fig. S5D).

We next assessed the tumor specificity of all AS variants. Although the majority of transcripts (34,416, 67.33%) were expressed in both NPC tumors and immortalized NP cell lines and 6487 (12.69%) were lowly expressed in all samples, a significant fraction (9349 transcripts, 18.29%) was exclusively detected in NPC tumor samples (Supplemental Table S6). In contrast, only 1.69% (863 transcripts) were specific to immortalized NP cell lines (Fig. 4D). Tumor-specific novel transcripts constituted a substantial proportion of this group, suggesting they may play a functional role in tumor biology and represent a reservoir of novel biomarkers.

To validate these findings, we selected a subset of highly expressed, tumor-specific novel isoforms that showed TPM ≥ 10 in at least one tumor sample and harbored isoform-specific novel junctions with short-read splice junction coverage ≥10. A heatmap of their expression patterns in our RNA-seq cohort confirmed their tumor-specific expression, although substantial intertumor heterogeneity was observed (Fig. 4E). This tumor specificity was further validated by RT-qPCR (Fig. 4F), and the isoform structures were depicted in Supplemental Figure S6. We further validated these results by integrating and analyzing three public NPC RNA-seq data sets (GEO: GSE118719, GSE134886, GSE68799), comprising 52 NPC tumors and 11 nontumor controls after batch-effect correction (Supplemental Fig. S7A). The expression patterns of the novel transcripts in these independent cohorts were consistent with our in-house data, confirming their specific and elevated expression in NPC tumors (Supplemental Fig. S7B).

A computational pipeline identifies potential tumor-specific neojunction-derived neoantigens for NPC immunotherapy development

To expand the repertoire of immunotherapeutic targets, we investigated neojunctions derived from tumor-specific novel transcripts. These neojunctions represented a source of noncanonical tumor-specific antigens, significantly broadening the landscape of potential targets for cancer therapy. Among all detected splice junctions, the majority (94.10%) were known, whereas 38,193 (5.90%) were novel (Fig. 5A). Most transcripts identified with novel splice junctions contained only a single novel junction, with progressively fewer transcripts carrying multiple novel junctions. Notably, most novel junctions were supported by canonical splice sites (85.3%), whereas only a minority involved noncanonical splice sites (14.7%), indicating that novel splicing largely conformed to established splice site rules (Fig. 5B). Moreover, most novel junctions were derived from novel transcripts, with few originating from intergenic, genic, fusion, or antisense regions (Fig. 5C).

Figure 5.

Identification of potential neojunction-derived neoantigens from tumor-specific novel transcripts using TS-SNAD. (A) Pie chart showing the distribution of splice junction categories detected by LR-seq, with known splice junctions (blue) and novel splice junctions (yellow). (B) Bar plot showing the number of transcripts with different numbers of novel splice junctions. Inset pie chart depicting the proportion of canonical and noncanonical splice sites. (C) Bar plot showing the distribution of transcripts containing novel splice junctions, categorized as NNC, NIC, intergenic, genic, fusion, or antisense. (D) Schematic overview of the TS-SNAD computational workflow for identifying neojunction-derived neoantigens from tumor-specific novel transcripts. (E) Ring plot showing that a total of 3326 neojunction-derived neoantigens were identified, of which 533 were predicted to bind strongly to at least one NPC HLA class I allele. (F) Bar plot showing the distribution of 533 strong-binding neoantigens by the number of tumor-specific novel transcripts from which they were derived. (G) Bar plot showing the number of distinct HLA class I alleles bound by the 533 neoantigens.

Open in new tabDownload PowerPointLink to figure

To systematically identify tumor-specific neoantigens derived from these neojunctions, we developed a customized computational workflow, TS-SNAD. This pipeline identified novel splice junctions from tumor-specific novel transcripts by comparison with reference genome annotations, as well as excluding any splice junctions present in normal tissues from GTEx short-read RNA-seq and long-read RNA-seq data from 88 samples from GTEx tissues (Glinos et al. 2022), thereby reducing the likelihood of selecting peptides with potential expression in vital normal organs and minimizing off-tumor toxicity risk. Tumor specificity was further supported by significantly higher read coverage of these novel splice junctions in NPC tumor samples compared with four immortalized NP cell lines in our cohort. The pipeline subsequently generated neojunction-spanning 9-mer peptides as candidate neoantigens. The binding affinity of these candidate neoantigens to HLA-I molecules was predicted using NetMHCpan, with an IC₅₀ < 500 nM and a percentile rank <0.5 set as selection criteria (Fig. 5D). To assess the binding potential of neojunction-derived neoantigens, we compiled HLA-I genotypes from 70 NPC patients from our in-house cohort with whole-genome sequencing data (Bruce et al. 2022) and predicted HLA-I types from 113 Chinese patients with RNA-seq data from the GSE102349 data set. Despite the high polymorphism of HLA-I, both cohorts showed a consistent presence of high-frequency alleles, including HLA-A*02:07, HLA-A*11:01, HLA-B*46:01, HLA-B*40:01, HLA-C*01:02, and HLA-C*07:02 (Supplemental Table S7), likely reflecting the population-specific genetic background.

Using TS-SNAD, we identified 944 tumor-specific novel transcripts harboring novel splice junctions with read coverage enriched in tumor samples. From these, we generated a total of 3326 neojunction-derived 9-mer neoantigens, of which 242 neojunction-derived 9-mer neoantigens matched peptide sequences encoded by other coding regions in the UniProtKB/Swiss-Prot human proteome. Of the 3326 neojunction-derived 9-mer candidates, 533 candidate neoantigens were predicted to bind patient-specific HLA-I alleles with high affinity, suggesting a high likelihood of being presented (Fig. 5E). Notably, most of these neoantigens were derived from unique tumor-specific novel transcripts, whereas 70 were produced by multiple isoforms (Fig. 5F). Analyzing the diversity of HLA-I interactions, we found that most neoantigens could bind multiple HLA-I alleles (Fig. 5G), suggesting their potential to be presented across a broader patient population with diverse immune contexts, thereby enhancing their therapeutic relevance.

Neojunction-derived neoantigens effectively elicit immunogenic responses

To evaluate the potential of neojunction-derived neoantigens as targets for broadly applicable immunotherapies in NPC, we selected 34 neojunction-derived peptides with high predicted binding affinity for HLA-B*40:01 and assessed their immunogenicity. The corresponding HLA-B*40:01 allele was carried by ∼29.20% of patients with NPC. Using peripheral blood mononuclear cells (PBMCs) from three HLA-B*40:01-positive healthy donors, we measured interferon-γ (IFNG) responses following stimulation with peptide-MHC (pMHC) complexes. Six of the 34 predicted neojunction-derived neoantigens induced robust and reproducible IFNG responses in all three donors (Fig. 6A; Supplemental Fig. S8A; Supplemental Table S8). The absence of response in DMSO controls and the robust activation induced by a CD3-specific antibody confirmed the specificity and reliability of the assay. The most immunogenic peptide, REPFVQAKL, was derived from a novel TYK2 transcript (TCONS_00025200) generated through skipping of exons 7–9. We observed considerable interdonor variability in the magnitude of IFNG responses, with donor 3 showing consistently stronger responses across all peptides, likely reflecting differences in individual immunological history.

Figure 6.

Neojunction-derived neoantigens effectively elicit immunogenic responses. (A) ELISpot assay results demonstrating that six out of the 34 splicing-derived neoantigens elicited immunogenic responses. PBMCs from three independent HLA-B*40:01 donors were stimulated with peptides (20 µg/mL), and IFNG secretion was measured by ELISpot (n = 3). DMSO served as the negative control, and anti-CD3 mAb as the positive control. Data are presented as mean ± SD. (B) Schematic illustrations of the exon–intron structures of the novel transcripts from which the six immunogenic neoantigens were derived.

Open in new tabDownload PowerPointLink to figure

Of note, each immunogenic neoantigen originated from a distinct tumor-specific novel transcript in our cohort, with the majority of underlying novel splice junctions resulting from exon skipping or cassette exon inclusion (Fig. 6B; Supplemental Fig. S8B(a)). The novel transcripts giving rise to these immunogenic neoantigens exhibited multiple features supporting their translational potential, including canonical splice sites, methionine-initiated ORFs, and no predicted susceptibility to NMD. Analysis of public NPC data sets confirmed that these novel isoforms exhibited high tumor-specific expression (Supplemental Fig. S8B(b)), minimizing potential risks of off-target reactivity against control tissues. Neojunction-derived neoantigens largely recapitulated the canonical HLA-B*40:01 binding motif, with a dominant glutamate at P2 and hydrophobic residues at the C terminus (Supplemental Fig. S8C).

As each neoantigen was found to be present in only 10%–20% of NPC patients, these results support the development of personalized multiepitope vaccines tailored to individual splicing profiles to maximize immune targeting and overcome immune evasion. Furthermore, we confirmed the presence of neojunction-derived unique peptides in NPC samples via mass spectrometry (Fig. 7A). It is noteworthy that junction-spanning peptides are typically underrepresented in mass spectrometry data sets owing to the cleavage specificity of trypsin (Kahles et al. 2018), underscoring the significance of these findings. Additionally, to assess the protein-coding potential of these novel RNA isoforms that generated validated immunogenic neoantigens, we overexpressed their ORFs fused to a 3×FLAG tag in both HEK293T and NP69 cells. Western blot analysis showed that all selected novel ORFs generated detectable proteins in vitro (Fig. 7B). Collectively, our results demonstrated that neojunction-derived neoantigens constituted a rich source of immunogenic targets, highlighting their strong potential for NPC immunotherapy development.

Figure 7.

Proteomic and western blot validation of ORFs encoding immunogenic neojunction-derived neoantigens in NPC. (A) Mass spectrometry provided protein-level evidence supporting the existence of neojunction-derived neoantigens. (B) Western blot analysis showing the expression of FLAG-tagged ORFs encoding the validated immunogenic neojunction-derived neoantigens in HEK293T and NP69 cells after transfection.

Open in new tabDownload PowerPointLink to figure

Long-read library preparation and sequencing

A total of 18 samples were collected and subjected to LR-seq, including four immortalized nasopharyngeal cell lines (NP361, NP460, NP550, and NP69) as controls, and 14 NPC samples comprising two cell lines (NPC43 and NPC53), 10 PDXs (X2117, X17, C15, X111, X113, X666, X76, X47, X23, and X32), and two tumor tissues (NPC-T9 and NPC-T64). This sample set enabled transcript discovery across diverse biological sources. Following RNA extraction, full-length cDNA was synthesized from poly(A)-tailed RNA using the SMARTer PCR cDNA synthesis kit (Clontech/Takara 634926). The cDNA was PCR-amplified to yield 1 to 2 µg of product and purified using AMPure XP magnetic beads (Beckman Coulter A63881). Size selection was performed using the Sage Science BluePippin system to remove short cDNA fragments and reduce preferential loading and sequencing of shorter molecules. SMRTbell adapters were ligated to cDNA ends, and libraries were purified using magnetic beads with the SMRTbell template prep kit (PacBio), followed by sequencing on PacBio Sequel II platform.

RNA sequencing

RNA was extracted using the RNeasy mini prep kit (Qiagen 74106), quantified on a Qubit 2.0 fluorometer (Thermo Fisher Scientific), and assessed for quality using a 2100 Bioanalyzer (Agilent) with the RNA 6000 Pico kit. cDNA libraries were prepared using the Illumina TruSeq stranded mRNA library prep kit. The final libraries were sequenced on an Illumina NovaSeq 6000 platform to generate 150 bp paired-end reads.

Quantitative PCR with reverse transcription

Reverse transcription was performed with a high-capacity RT kit (Applied Biosystems, Thermo Fisher Scientific 4368813). The resulting cDNA was diluted 1:3 using nuclease-free water, and 2 µL was used per quantitative PCR (qPCR) reaction. qPCR with reverse transcription was performed using the POWER SYBR green master mix (Applied Biosystems, Thermo Fisher Scientific 4367659). All samples were analyzed in technical triplicates using the QuantStudio 5 (Thermo Fisher Scientific), and all gene expression data were normalized to human GAPDH. Primer sequences were listed in Supplemental Table S9.

LR-seq data analysis

Raw PacBio subreads were processed using the Iso-Seq3 pipeline (version 3.8.0; https://github.com/PacificBiosciences/IsoSeq/) to generate high-quality, full-length transcript sequences for downstream analysis. Briefly, raw subreads were processed to generate CCS reads with a minimum predicted accuracy of 0.99. Full-length reads were generated by removing the 5′ and 3′ cDNA primers and performing demultiplexing using the “lima” module with the default parameters. Artificial concatemers and poly(A) tails were removed using “isoseq3 refine” to produce full-length, nonchimeric (FLNC) reads. These FLNC reads were clustered into high-quality transcripts using “isoseq3 cluster.” The filtered transcripts were aligned to the human reference genome (GRCh38/hg38) using pbmm2 (version 1.10.0; https://github.com/PacificBiosciences/pbmm2). Redundant isoforms were collapsed within each sample using the collapse_isoforms_by_sam.py script from the cDNA_Cupcake ToFU toolkit (https://github.com/Magdoll/cDNA_Cupcake). Finally, transcripts from all samples were merged into a unified, nonredundant transcriptome using gffcompare (version 0.11.2) (Pertea and Pertea 2020) for subsequent analysis (see Supplemental Fig. S1A).

Isoform annotation and quality control

Full-length isoforms were annotated and subjected to quality control using SQANTI3 (version 5.2) (Pardo-Palacios et al. 2024) with GENCODE v38 as the reference annotation. SQANTI3 was used to annotate isoforms and generate quality-control metrics, including CAGE peaks (Noguchi et al. 2017), poly(A) sites (Moon et al. 2025), and NMD predictions, to help filter potential artifacts and ensure high-quality transcript annotations. Isoform validity was further supported by short-read RNA-seq data. The coverage of isoform splice junctions was calculated using the junction file (commonly referred to as SJ.out.tab) generated from RNA-seq data. ORFs were predicted with ORFfinder (https://github.com/Chokyotager/ORFFinder), retaining only those encoding 100 or more amino acids.

De novo transcriptome assembly

Sorted BAM files (hg38-aligned) were assembled de novo using StringTie (version 2.1.1) (Pertea et al. 2015), with the following parameters: -m 200 -f 0.05 -c 1.5 -p 10. Individual assemblies from the NPC samples were merged with StringTie ‐‐merge with the default parameters. Both the de novo and the LR-seq-derived transcriptomes were compared with GENCODE v38 using gffcompare -r.

Short-read RNA-seq analysis

The quality control of short-read RNA-seq data was performed using FastQC (version 0.11.9). The low-quality reads and adapters were filtered and trimmed by using cutadapt (version 4.4) (Kechin et al. 2017). Clean reads were mapped to the human reference genome GRCh38 using STAR (version 2.7.10b) (Dobin et al. 2013). Specifically, a two-pass strategy was used to extract splice junctions from each sample. The TPM abundances of all isoforms were calculated using RSEM (version 1.3.1) (Li and Dewey 2011) using the LR-seq-derived transcriptome as reference. Public RNA-seq data from different cohorts were integrated, followed by batch-effect removal using ComBat from the R sva package.

ASE analysis

We used SUPPA2 (version 2.3) (Trincado et al. 2018) to calculate seven types of ASEs including A3, A5, AF, AL, SE, RI, and MX. To incorporate novel transcripts, a unified GTF file (merged from GENCODE v38 and long-read-derived novel isoforms) was used with the generateEvents command (-f ioe). Differential splicing analysis was performed with SUPPA2 diffSplice, with significance thresholds set at |ΔPSI| > 0.1 and P < 0.05. To improve accuracy, splice sites were held to stringent boundaries, whereas transcription start/end sites (AF/AL events) allowed flexible boundaries (±48 nt).

Pathway enrichment analysis

Pathways significantly enriched (false-discovery rate [FDR] < 0.05) among genes with novel isoforms detected by LR-seq were identified using the MSigDB C2, C5, and Hallmark collections (Liberzon et al. 2015). Enrichment analyses were performed using the R package clusterProfiler (version 3.16.1) (Wu et al. 2021). A matched background gene set was used for all enrichment tests and was defined as the intersection of (1) genes detectable in our long-read data and (2) genes that could be mapped to the selected MSigDB collections. This background ensured that enrichment was evaluated relative to genes testable in our experiment and reduced potential bias.

To control for potential bias arising from gene expression levels, we performed an expression-matched permutation analysis. For each gene, expression level was defined as the median TPM across samples within the expressed gene universe. Genes were stratified into quantile-based expression bins (10 bins by default). For the observed gene set generating novel isoforms, we recorded the number of genes in each expression bin and generated matched random gene sets with the same expression-level distribution and repeated this procedure across 200 permutations. Pathway enrichment analysis was then performed for each randomized gene set. For each pathway significantly enriched in the genes with novel isoforms, an empirical P-value (p_emp) was calculated as the fraction of expression-matched permuted gene sets whose P-value (p_perm) was less than or equal to the observed P-value (p_obs) obtained for the gene set with novel isoforms:

DTE analysis

RSEM expected read counts were used as input for downstream differential expression analysis. DTE between NPC tumor samples and immortalized nasopharyngeal cell lines was assessed in DESeq2 (v1.38.3) (Love et al. 2014) using the Wald test under a negative binomial generalized linear model, with Benjamini–Hochberg correction for multiple testing. Transcripts with |log₂FC| > 2 and adjusted P < 0.05 were considered significantly differentially expressed.

Identification of transcripts encoding putative cell surface transmembrane proteins

An isoform was classified as encoding a putative cell surface transmembrane protein if its ORF contained at least one transmembrane helix predicted by DeepTMHMM (Hallgren et al. 2022) and was predicted by DeepLoc (Ødum et al. 2024) to localize to the cell membrane. For external validation, we used the Cancer Surfaceome Atlas (TCSA) as a gene-level reference set of genes encoding cell surface proteins. TCSA integrates nine independent resources to define a curated cancer surfaceome. Because currently available databases generally lack isoform-resolved subcellular localization annotations, we could not reliably assign cell surface localization to a single canonical isoform for every gene. Therefore, the gain or loss of predicted cell surface localization was assessed at the isoform level relative to the annotated protein products of the corresponding gene.

Identification of splicing-derived neoantigens by TS-SNAD

Tumor specificity was assessed by integrating long-read and short-read RNA-seq data, as the sequencing depth of our long-read data set was insufficient for robust transcript-level abundance estimation. Long-read RNA-seq was used primarily to construct a LR-seq-derived transcriptome by resolving full-length transcript structures and identifying both annotated and previously unannotated isoforms, whereas short-read RNA-seq was used to quantify transcript expression based on the LR-seq-derived transcriptome. Tumor-specific transcripts were defined as transcripts with TPM ≥ 1 in at least one tumor sample and no detection in any of four immortalized NP cell lines in our cohort. Splice junction coverage was quantified using SJ.out.tab files generated by STAR alignment of RNA-seq data to the GRCh38 genome, with the LR-seq-derived transcriptome annotation provided as a reference. For neojunction prioritization, all novel splice junctions derived from tumor-specific novel transcripts were subjected to further filtering. Candidate neojunctions were required to meet the following criteria: (1) the average junction read coverage in tumor samples was higher than that in immortalized nonmalignant NP cell lines; (2) fewer than five reads support the novel splice junctions in each immortalized NP cell line; and (3) fewer than five reads support the novel splice junctions in each of 11 nonmalignant nasopharyngeal tissue samples from three public NPC RNA-seq data sets (GEO: GSE118719, GSE134886, and GSE68799); and (4) the novel splice junctions were required to be absent from all GTEx data sets, including short-read RNA-seq from normal tissues and ONT long-read data from 88 GTEx tissues and cell lines.

For neoantigen prediction, 9-mer peptides spanning neojunctions were extracted. NetMHCpan (version 4.0) (Jurtz et al. 2017) was then used with the default parameters to predict HLA binding affinity, using neojunction-derived 9-mer peptides and patient-specific HLA genotypes as input. Only candidate neoantigen peptides with strong predicted HLA-I binding were retained, defined as those with a binding affinity percentile rank <0.5 and predicted IC₅₀ < 500 nM.

HLA-I haplotype prediction

HLA-I genotypes for 70 patients from our local cohort were determined from our previous whole-genome sequencing data and subsequently used for neoantigen binding affinity prediction. For public data sets, HLA-I genotypes were predicted from RNA-seq FASTQ files using OptiType (version 1.3.3) (Szolek et al. 2014). We selected OptiType based on its demonstrated high accuracy (>99%) in predicting HLA-I alleles from RNA-seq data.

IFNG ELISpot assay

The immunogenicity of each neoantigen peptide was evaluated using an IFNG ELISpot assay. Briefly, PBMCs from donors carrying the HLA-B*40:01 allele were isolated from whole blood using Ficoll-Paque (Cytiva 17144002) density gradient centrifugation. The collected PBMCs were then resuspended in RPMI-1640 (Gibco 11875093) complete medium supplemented with 10% fetal bovine serum (Gibco A5256701) and incubated at 37°C in a humidified incubator with 5% CO₂ to allow cell recovery and stabilization. Subsequently, 3 × 10⁵ PBMCs in 100 µL were seeded into each well of a 96-well ELISpot plate (Mabtech 3420-4AST-2), followed by the addition of 2 µL of peptide stimulator (final concentration 20 µg/mL). Anti-CD3 monoclonal antibody was used as a positive control, and medium containing an equivalent concentration of DMSO served as a negative control. Plates were then incubated in a humidified 5% CO₂ incubator for at least 20 h without disturbance at 37°C. Measures were taken to prevent evaporation during incubation. After incubation, cells were removed, and the plates were washed five times with PBS. A biotinylated anti-IFNG detection antibody (1:1000 dilution) was added to each well and incubated for 2 h at 37°C. After another round of washing, streptavidin-ALP (1:1000) was added and incubated for 1 h at room temperature. Finally, BCIP/NBT substrate solution was applied to develop colorimetric spots representing IFNG-secreting cells. Spot-forming units (SFUs) were quantified using an automated ELISpot reader (AID iSpot Spectrum). Responses were considered positive if spot counts exceeded background by at least threefold.

Cell cultures, plasmids, and transfection

NP69 cells were cultured in keratinocyte SFM (1×) kit (Gibco 17005042), supplemented with 1% penicillin–streptomycin (10,000 U/mL, Gibco 15140122). HEK293T cells were maintained in DMEM (high glucose, Gibco 11965118) supplemented with 10% FBS (Sigma-Aldrich F2442-500ML) and 1% penicillin–streptomycin. All cell lines were cultured with 5% CO₂ at 37°C and were routinely tested and confirmed to be free of mycoplasma contamination.

ORFs of the novel transcripts, followed by a 3×FLAG tag, were cloned into the pcDNA3.1(+) vector, and all plasmid constructs were verified by sequencing. HEK293T cells were transiently transfected using Lipofectamine 3000 reagent (Invitrogen L3000001) according to the manufacturer's instructions. NP69 cells were transfected with FuGENE HD (Promega E2311) according to the manufacturer's instructions.

Reagents and antibodies for immunoblotting

The following antibodies were obtained from the indicated sources and used for immunoblotting: anti-FLAG M2 mouse mAb (Sigma-Aldrich F1804) and anti-β-actin rabbit antibody (ABclonal AC038). Secondary antibodies were obtained from the following sources: HRP-conjugated goat antimouse IgG (H+L) (Invitrogen 31430) and HRP-conjugated goat anti-rabbit IgG (H+L) (ABclonal AS014).

Immunofluorescence

Cells grown on coverslips were fixed with 4% paraformaldehyde for 10 min at room temperature and washed with PBS. Fixed, nonpermeabilized cells were then incubated with wheat germ agglutinin-AF488 (MCE, HY-NP163A) in PBS for 10 min at room temperature to label the plasma membrane. Following WGA staining, cells were washed with PBS and subsequently permeabilized with 0.2% Triton X-100 in PBS for 5 min. Cells were then blocked with 5% BSA in PBS for 60 min at room temperature. Cells were incubated with anti-FLAG M2 mouse monoclonal antibody (Sigma-Aldrich F1804) overnight, followed by incubation with goat antimouse IgG (H+L) secondary antibody (Invitrogen A11004) for 60 min. Nuclei were stained with DAPI. Coverslips were mounted using antifade mounting medium, and images were captured using a Leica Stellaris 8 Falcon confocal microscope.

Proteomic profiling and peptide database search

Public NPC mass spectrometry data were obtained from the iProX partner repository (https://www.iprox.cn/page/SCV017.html) under accession number IPX0001265000. Peptide identification was performed using MS-GF+ (version 2018.10.15) (Kim and Pevzner 2014) against a custom sequence database comprising a total of 42,613 merged ORFs, including 22,017 ORFs derived from long-read novel transcripts and 20,596 reviewed human protein sequences from UniProtKB/Swiss-Prot database. The following parameters were set for database searching: carbamidomethyl (C), iTRAQ4plex (N-term), and iTRAQ4plex (K) were specified as fixed modifications; oxidation (M), deamidated (NQ), acetyl (K), and methyl (K) were specified as variable modifications. The precursor mass tolerance was set to 20 ppm, and the product ion tolerance for MS/MS was 0.05 Da. Trypsin was specified as the digestion enzyme, allowing up to two missed cleavages. mzID profiles identified from the search engine were then pooled using the R/Bioconductor package MSnbase, and peptide-to-spectrum matches (PSMs) satisfying both spectra and peptide FDR cutoffs <1% were retained for further analysis. To validate neojunction-derived neoantigens by mass spectrometry, only unique 9-mer peptides with unambiguous spectral matches were considered.

Survival analysis

Survival analysis was performed to evaluate the prognostic associations of transcript expression levels in NPC patients. Patients were stratified into high- and low-expression groups based on the median transcript levels. Kaplan–Meier survival curves were generated and compared using the log-rank test with the survival (version 3.3) and survminer (version 0.4.9) R packages.

Data sets

Public RNA-seq data sets from NPC cohorts were obtained from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) with the accession numbers GSE118719, GSE134886, GSE68799, and GSE102349. RNA-seq splice junction data from normal human tissues were retrieved from the GTEx project (https://www.gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression). Long-read RNA-seq data for 88 GTEx tissue and cell lines, based on the ONT platform, were also retrieved from the GTEx project (https://www.gtexportal.org/home/downloads/adult-gtex/long_read_data). Public mass spectrometry proteomic data of NPC were downloaded from the iProX partner repository (number IPX0001265000).

Ethics declarations

For the two NPC tumor samples, written patient consents were obtained from the patients according to institutional clinical research approval. The study protocol was approved by the joint Chinese University of Hong Kong-New Territories East Cluster clinical research ethics committee at the Chinese University of Hong Kong, Hong Kong SAR. The collection and use of PBMCs from healthy donors were approved under the protocol “collection of peripheral blood mononuclear cells for new drug research and development” (protocol no. LP202210). Written informed consent was obtained from all PBMC donors before sample collection.

Code availability

All code required to reproduce the data processing and downstream analyses in this study has been provided as Supplemental Code files and archived in a publicly accessible GitHub repository. Within this repository, the analysis scripts are available at GitHub (https://github.com/CityUHK-CompBio/NPC-LR-Seq/tree/main/script), and the TS-SNAD pipeline is available at GitHub (https://github.com/CityUHK-CompBio/NPC-LR-Seq/tree/main/TS-SNAD). The TS-SNAD directory includes the full pipeline implementation, example input and output files, parameter descriptions, and a step-by-step usage guide.