The cohort contains molecular profiles from the Cancer Genome Atlas (TCGA) PDAC dataset containing transcriptome data from 185 patient samples (14). We then retrieved for TCGA cohort with multi-omics profiles, including transcriptome expression, somatic mutations, copy number alterations (CNAs), clinicopathological features and clinical outcomes. From UCSC Xena (http://xena.ucsc.edu/), we downloaded gene expression data for 165 primary PDAC cases. Transcripts per kilobase million (TPM) values were then converted from FPKM to improve comparability between samples and show greater similarity to microarray results (15). GENCODE22 mapping from UCSC Xena was used to convert Ensembl IDs for transcriptomes into gene symbols. Furthermore, a copy number segment file was obtained from FireBrowse (http://firebrowse.org/). By using cBioPortal, we collected data on somatic mutations, clinicopathological features, progression-free survival (PFS), and overall survival (OS). (https://www.cbioportal.org/).

To test the reproducibility of our identified molecular subtypes, we collected a total of five external independent cohorts, including PACA_AU (n = 269) which was downloaded from International Cancer Genome Consortium (ICGC) (16), E-MTAB-6134 containing 309 consecutive patients with PDAC from a multi-center consortium (17), and three cohorts from the Gene Expression Omnibus (GEO) database; the GEO datasets included GSE71729 (n = 125) (18), GSE62452 (n = 66) (19), and GSE57495 (n = 63) (20). From the corresponding archives, we retrieved clinical outcomes only for PDAC patients and excluded those with neuroendocrine or acinar cell carcinomas. When multiple probe IDs were associated with a gene symbol in three microarray datasets from the GEO database, the median value was taken into account. The three independent GEO datasets were combined considering the relatively small sample size. Potential cross-dataset batch effects were eliminated by the R package sva under an empirical Bayesian framework and batch effects were further investigated using principal component analysis (21).

To many molecular components, NAT10 was known as an acetyltransferase. A recent study discovered NAT10’s ability to acetylate mRNA, providing new insight into the epitranscriptome. The literature yielded a list of 2,156 genes that are altered by NAT10-mediated ac4C (22).

The R package MCPcounter was used to assess TME composition in PDAC patients (23). An analysis of transcriptomic markers determines the scores of a given cell population, which are transcriptomic features that are strongly, specifically, and consistently expressed. Based on these scores, large cohorts and inter-sample comparisons can be made between tumor populations (24). We calculated infiltrating immune/stromal cells and tumor purity in tumor tissue using the R package “estimate” (25). In practice, it is frequently difficult to obtain cancer samples with sufficient tumor purity, particularly for PDAC. In this context, we kept those tumor samples with tumor purity of greater than 70% to ensure sufficient tumor cells and microenvironment cells for analytic purposes.

The functional orientation of the TME was determined using literature signatures (26). Immune suppression (CXCL12, TGFB1, TGFB3, and LGALS1), T cell activation (CXCL9, CXCL10, CXCL16, IFNG, and IL15), T cell survival (CD70 and CD27), regulatory T cells (FOXP3 and TNFRSF18), major histocompatibility complex class I (HLA-A, HLA-B, HLA-C, HLA-E, HLA-F, HLA (CXCL13). The geometric mean signature expression was used to compute scores for each signature. Additionally, DNA methylation scores were calculated using protocols described in the literature for tumor-infiltrating lymphocytes (MeTILs) (27).

Unsupervised hierarchical clustering of samples in the discovery TCGA cohort was performed using Euclidean distance and Ward linkage criteria based on polygenic abundance scores from MCPcounter included cohorts (28). Based on the dendrograms and the relatively small sample size, three clusters were empirically determined.

Using single-sample gene set enrichment analysis (ssGSEA), we identified 21 replication stress signatures from the literature (29). Hierarchical clustering identified two replication stress subtypes.

A differential expression analysis was performed using the R package “limma” (30). We created a pre-ranked gene list based on log2FoldChange values derived from differential expression analysis for gene set enrichment analysis (GSEA). The Hallmark pathway was used to determine functional enrichment using the R package clusterProfiler (31, 32).

We performed comprehensive characterization on subtypes identified from multi-omics aspect using the R package MOVICS (33), including prognosis, mutational landscape, chromosomal stability and clinicopathological features. Specifically, for mutational landscape, significantly, differentially, and frequently mutated genes among subtypes were identified (mutational rate > 3% with P < 0.05). The individual fraction of copy number-altered genome (FGA) was calculated to represent chromosomal instability based on copy number segment data using the threshold of 0.2 (33). In addition, we detected and localized recurrent focal somatic CNAs by GISTIC2.0 through GenePattern with criteria by default as the following: threshold for amplification or deletion: 0.1, q-value: 0.25 with 75% confidence level (https://www.genepattern.org/) (34).

In order to classify genes based on gene expression, nearest template prediction was used (NTP). This prediction algorithm provides a convenient model-free approach that uses only feature gene lists and test datasets for single-sample class prediction and is flexible and beneficial in external cohort applications (35, 36).

We used the R package to analyze chemical compounds. Drug sensitivity for PDAC was predicted for the TCGA cohort using the R ‘pRRophetic’ package based on GDSC drug sensitivity and phenotypic data (https://www.cancerrxgene.org/). In the R package, ridge regression was used to estimate the AUCs for each sample treated with a particular chemotherapeutic agent (37).

For immunotherapy especially to immune checkpoint inhibitors, we collected a transcriptome profile and clinical response of 47 patients with melanoma who treated with either anti-PD1 or anti-CTLA4 immunotherapies (38). A comparison of transcriptome profiles enabled us to predict the clinical response to immune checkpoint inhibitors using subclass mapping (39).

A total of 179 PDAC and 332 standard pancreatic tissue specimens from the TCGA and GTEx databases, as well as OS, PFS, DFS, and Disease Special Survival (DSS), were selected to assess the correlation between differential expression of NAT10 and prognosis. The human protein atlas: This database, also known as the human protein expression atlas, is a study of the immunohistochemical staining status of proteins in normal human tissue, cancer tissue and cancer cell lines, allowing for differential protein expression analysis in tumor tissues, and for Gene and tumor survival analysis. In this study, immunohistochemical analysis of NAT10 protein in PDAC tissues was performed using The human protein atlas database.

we analyzed the expression of NAT10 in tumor and normal tissues, human Pancreatic ductal adenocarcinoma cell lines (Capan-1, CFPAC, PANC-1, MIAPaCa-2 and BXPC-3) and pancreatic duct epithelial tissue cell line (HPDE6-C7) were obtained from BeNa Culture Collection (Beijing, China). Capan-1 was maintained in Iscove’s Modified Dulbecco’s Medium (IMEM) supplemented with 20% Fetal Bovine Serum (FBS) and HPDE6-C7 was maintained in Dulbecco’s Modified Eagle Medium (DMEM) supplemented with 10% FBS, at 37°C in 5% CO₂.

All samples used in this study were approved by the Ethics Committee of Shanghai Outdo Biotechnology Co., Ltd (No. HPanA120Su02). The operation time was from January 2009 to August 2013, and the follow-up time was from November 2014, with a follow-up period of 1.2 to 5.8 years. A total of 66 cases of cancerous tissue and 54 cases of paracancerous tissue were obtained for immunohistochemistry. Pancreatic ductal adenocarcinoma tissues collected were fixed with 4% paraformaldehyde, dehydrated, embedded into paraffin blocks and sliced into 4μm. The sections were hydrated after dewaxing and incubated with 3% H₂O₂. After antigen repair and being blocked, slides were blocked and incubated overnight at 4°C with rabbit antibody specific for NAT10 (1:100) (Abcam, USA). For 20 minutes, the sections were incubated with a streptavidin-biotin-peroxidase complex, followed by 30 minutes with the secondary antibody. We stained the slide with 3, 3-diaminobenzidine (DAB) substrate kit for a peroxidase reaction and counterstained with hematoxylin. The slides were analyzed under an optical microscope with a digital camera.

The total RNA was extracted using the Trizol kit according to the instructions, and its concentration and purity were checked using the Nano-300 microspectrophotometer. Cells were used for subsequent experiments after qRT-PCR was used to measure the expression levels of target RELATIVE NAT10 mRNA. The used specific primer sequences were as follows: NAT10-forward primer (GAGACAGACCCCGAATGACC), reverse primer (GGAGAGCAAGGCTAGGAACC) and GAPDH-forward primer (GGAGCGAGATCCCTCCAAAAT), reverse primer (GGCTGTTGTCATACTTCTCATGG).

Cells were collected and lysed on ice with a mixture of RIPA lysate(Thermo Fisher, USA) and PMSF(Beyotime, ST506, China), total protein was extracted from each group of cells, and the BCA kit(Thermo Fisher, USA) was used to determine the protein content of each group. 40μg of the sampled proteins were subjected to SDS-PAGE, 250 mA constant flow membrane transfer (0.22 μm PVDF), and closed with TBST blocking solution (containing 50 g/L skimmed milk powder) for 2 h. Then the following primary antibodies were incubated overnight at 4°C: rabbit monoclonal anti-NAT10 (1 : 1000, ab194297, Abcam), rabbit monoclonal anti-β-actin (1 : 1000, 4970, Cell Signaling Technology), the next day the cells were rinsed with TBST and the secondary antibody was added for 2 h at room temperature. Using an ECL kit(Thermo Fisher, USA), the expression of the target bands was detected. Use the ImageJ software (version 1.50b; National Institutes of Health) for grayscale analysis. Each experiment was repeated three times to minimize experimental error.

Small interfering RNA (siRNA) was transfected into Capan-1 cells by Lipofectamine 3000. Cells were collected, and NAT10 protein expression levels were determined using Western blot. siRNAs were obtained from Shanghai Genepharma. Sequences of siRNA were as following: siRNA1 NAT10:sense 5′-GGCAGACUAUUCAGUAUAUTT-3′, anti-sense 5′-AUAUACUGAAUAGUCUGCCTT-3′; siRNA2 NAT10:sense 5′-GGCCAAAGCUGUCUUGAAATT-3′, anti-sense 5′-UUUCAAGACAGCUUUGGCCTT-3′; Negative control: sense 5′-UUCUCCGAACGUGUCACGUTT-3′, anti-sense 5′-ACGUGACACGUUCGGAAATT-3′.

Cells from the above groups were digested 48 h after transfection, collected in serum-free DMEM medium and counted, and diluted to 2.5×10⁴ cells/100 μL. 700 μL of complete medium (DMEM with 100 mL/LFBS) was added to the lower chamber, 8 μm diameter TranswellTM cells were placed in a 24-well plate, and 200 μL of the serum-free medium was added to the upper chamber. After incubation for 24h, the chambers were removed and fixed in 40g/L paraformaldehyde for 30min; 1g/L crystal violet was used for staining for 30min.

Using a 6 cm plate with 1000 cells inoculated and incubated at 37°C, clone formation was assessed by number of colonies formed. According to the characteristics of each cell line, the medium was changed regularly after 10-15 days. We washed the adherent cells twice in PBS, fixed them with 4% paraformaldehyde, and stained them with 0.1% crystal violet (Sigma-Aldrich, NY, USA). The total number of colonies (50 cells per colony) was counted using Image J software.

R was used to conduct statistical analyses (v4.0.2). P values or adjusted P values less than 0.05 were considered significant for all statistical comparisons. Mann-Whitney U tests were used to compare two conditions, and Kruskal-Wallis tests were used to compare three conditions. The independence of categorical variables was determined using Fisher’s exact tests. To analyze Kaplan-Meier survival curves and log-rank test calculations, two-sided log-rank P values were calculated using the R package survminer. Each experiment was repeated three times to minimize experimental error. Data are shown as mean ± standard error of the mean (SE). The 2^{-ΔΔ Ct} method was used to analyze the results of real-time PCR in all of the experiments.

To narrow down the gene list and search for the prognostication-relevant genes that were associated with NAT10, we conducted univariate analysis using Cox proportional hazard regression model in PDAC cases with sufficient tumor purity (n = 56), taking into account both PFS and OS. We identified 695 genes (Table S1) and 574 genes (Table S2) related to PFS and OS (P < 0.05), relatively, and a total of 458 genes were shared (Figure 1A). We then performed unsupervised clustering using these 458 NAT10-relegent prognostic genes and revealed three subtypes among PDAC cases in TCGA cohort (Figure 1B). The identified subtypes were tightly associated with PFS (P < 0.001; Figure 1C) and OS (P < 0.001; Figure 1D); notably, C3 showed the most inferior outcome as compared to other cases. We then investigated the prognostic value of our proposed NAT10-relavnt PDAC subtype, which we found remains an independent prognostic factor after adjusting for other clinical characteristics regarding PFS and OS (Figure 1E).

In order to determine whether our newly identified NAT10-relevant prognostic subtypes can be reproducible in external cohorts of patients with PDAC, we tested their reproducibility. In this manner, we searched for the top 30 NAT10-relevant genes that were significantly and uniquely upregulated (log2FoldChange > 0, FDR< 0.05) in each of the three subtypes from TCGA cohort, ending up with a 88-gene classifier in which 30 genes for C1, 28 genes for C2, and 30 genes for C3 (Table S3). We then deployed the 88-gene classifier to the discovery TCGA cohort using NTP algorithm (Figure 2A); we demonstrated superior predictive accuracy with only five misclassified (Kappa = 0.857, P < 0.001; Figure 2B), which made us confident of using such classifier in external cohorts.

In this manner, we first deployed the 88-gene classifier to the combined GEO cohort with sufficient tumor purity (n = 55) where potential batch effect was removed and investigated by PCA (Figure 2C, D). NTP revealed three NAT10-relevant subtypes in GEO cohort (Figure 2E), which showed significantly separated OS rate and the predicted CS3 presented with the most unfavourable prognosis (Figure 2F). Using the same strategy, we reproduced three subtypes in PACA_AU (n = 48) and E-MTAB-6134 cohorts (n = 156), respectively (Figures 2G, H). Likewise, the reproduced subtypes demonstrated prognostication-relevance (P = 0.052 in PACA_AU cohort and P < 0.001 in E-MTAB-6134 cohort) and CS3 showed inferior outcome as compared to other subtypes (Figures 2I, J). C1, C2, C3 were the molecular typing of the training set, and CS1,CS2,CS3 are the molecular typing of the validation set.

To further understand the biology-relevance behind each of the three subtypes, we performed differential expression analysis and GSEA using Hallmark pathways. Significantly and uniquely upregulated Hallmark pathways (NES > 0, FDR < 0.05) were identified for each of the three subtypes,we found that the pathway represented by Pancreas beta cells was up-regulated in the CS1 subtype, and the pathway represented by MYC Targets V2 was up-regulated in the CS2 subtype. And the pathway represented by Epithelial mesenchymal translation is up-regulated in CS3 subtype. (Figure 3A; Table S4). Of note, we found immune-related pathways were significantly activated in the aggressive CS3 subtype as compared to other cases; these pathways included inflammatory response (NES: 1.86, FDR = 0.002), interferon-gamma response (NES: 1.90, FDR = 0.002), and epithelial mesenchymal transition (EMT; NES: 2.20, FDR = 0.002). In addition, cell cycle-relevant pathways (i.e., E2F targets [NES: 1.64, FDR = 0.002] and G2M checkpoint [NES: 1.68, FDR = 0.002]) were significantly upregulated in CS3 comparing to other subtypes (Figure 3B).

As we have shown the immune dysfunction of among three subtypes, we then decided to in-depth profile the TME of PDAC. Therefore, we quantified global enrichment of immune/stromal cells and tumor-infiltrating lymphocyte based on DNA methylation data. We further investigated 10 TME cell infiltration levels in TCGA cohort, seven gene signatures for the functional orientation, and surveyed the PDAC samples for the expression of immune checkpoints (Figure 4A). We found that the aggressive C3 showed highly infiltrated with immune (Figure 4B), stromal cells (Figure 4C) and tumor-infiltrating lymphocytes (Figure 4D), and enriched in both immune activation and suppression factors. Moreover, C3 showed relatively higher expression of nearly all immunotherapy targets, including CD274, PDCD1, PDCD1LG2, CTLA4, and HAVCR2. As our findings demonstrated that C3 presented with highly infiltrated TME and activated interferon-gamma signaling pathway, we therefore hypothesized that C3 may have higher likelihood of responding to immune checkpoint inhibitor as compared to other subtypes. In order to accomplish this, we mapped subclasses in the TCGA cohort. According to the results, the transcriptome profiles of C3 and melanoma patients who responded to anti-CTLA4 immunotherapy were highly similar (P = 0.006, Bonferroni adjusted P = 0.072, FDR adjusted P = 0.072; Figure 4E). Our findings add to the literature that molecular classification using NAT10-relevant genes may be capable of identifying ideal candidates for immunotherapy in patients with PDAC.

Furthermore, the genomic landscape plays an integral role in anti-tumor immunity as well as molecular features. Several factors can trigger T-cell responses, including tumor mutational burden (TMB) and the presence of neoantigens (38, 40), while aneuploidy may result in immune evasion and reduced response to immunotherapy (41). To the end, we first explored the genetic differences among three subtypes. We dissected the mutational landscape across all the samples and calculated the TMB for each tumors; we found that CS1 harboured significantly less TMB as compared to other cases (P = 0.019; Figure 5A). We then identified a total of six mutations that showed significantly differential rate of mutational frequency among three subtypes (P < 0.05); these mutations occurred in at least 3% of the PDAC in TCGA cohort (Table S5). We then investigated if these genes included any mutations that were previously identified as driver mutations for cutaneous melanoma (42), and three mutations were identified, including KRAS (68%), TP53 (52%), and CDKN2A (27%). Strikingly, none of these three mutations occurred in C1 (Figure 5B).

We then investigated chromosomal instability by first profiling the broad-level CNA across the whole human gene in the entire cohort, and three subtypes we proposed, respectively (Figure 5C). Consistently, we found that CS1 showed generally higher chromosomal instability as compared to other subtypes by calculating the individual FGA values as well as fraction genome gained and lost (FGG and FGL) (Figure 5D). Next, we profiled the focal-level CNA for PDAC (Figure 5E) and we found that C1 showed significantly higher focal-level CN amplifications (P = 0.065; Figure 5F); while C3 enriched in CN deletions (P = 0.009; Figure 5G) comparing to other two subtypes.

We investigated the potential of targeting RS (Replication stress) as a PDAC subgroup treatment regimen because we found that cell cycle-related oncogenic pathways are activated in aggressive C3 and activation of cell cycle pathways activates cell cycle checkpoint regulatory proteins involved in RS such as ATR and WEE1, which have been described as closely related proteins associated with DNA damage responses leading to cisplatin resistance (29). Unsupervised hierarchical clustering using 21 replication stress markers resulted in two RS isoforms (RS-High and RS-Low) in the TCGA, GEO, PACA AU and E-MTAB-6134 cohorts (Figures 6A–D). There was a significant enrichment of aggressive C3 in RS-High (all, P 0.05), indicating chemoresistance based on cisplatin. We estimated AUC for three ATR inhibitors (VE-821, VE-822, and AZD6738) and two WEE1 inhibitors (WEE1 inhibitor and MK-1775) provided by GDSC using a model-based prediction strategy similar to the literature (43). Cell cycle checkpoint inhibitors were generally more effective against the RS-High subtype in all four cohorts than against the RS-low (Figures 6E–H).

Gemcitabine is the standard first-line treatment for unresectable locally advanced or metastatic pancreatic ductal adenocarcinoma since 1997 (44). In this context, we further investigated the association between NAT10-relevant subtypes and sensitivity of responding to gemcitabine. Firstly, we profiled the distribution of NAT10 in three PDAC subtypes, and we found that the single-sample level enrichment score of NAT10 tended to increase from C1 to C3 in four cohorts (Figures 7A–D). Secondly, we performed correlation analysis and we found that the predicted AUC of gemcitabine was positively correlated with the NAT10 enrichment score (Figures 7E–H), which indicated that PDAC patients with increasing NAT10 expression level may experience increasing resistance to treatment of gemcitabine.

Using the GEPIA website, NAT10 expression levels were analyzed between the TCGA and GTEx databases in pancreatic ductal adenocarcinoma tissues. Pancreatic ductal adenocarcinoma tissues expressed significantly higher levels of NAT10 than normal pancreatic tissues, and the difference was statistically significant (Figure 8A). A search from the Human Protein Atlas database showed that NAT10 expression was significantly higher in pancreatic ductal adenocarcinoma tissues (Figure 8B). An analysis of pancreatic ductal adenocarcinoma prognoses from the TCGA and GTEx databases can be found on the GEPIA website. According to the survival curve analysis, NAT10 expression levels were negatively correlated with OS, PFS, and DSS. Patients with abnormally high NAT10 expression had poor prognoses, with statistically significant differences (Figure 8C). AKT activation and cell proliferation characteristics were positively correlated (Figures 8D, E), suggesting that the abnormal expression of NAT10 may promote the malignant proliferation of pancreatic ductal adenocarcinoma by activating the PI3K-AKT pathway. On the other hand, NAT10 was positively correlated with TGF-β and angiogenic ability (Figures 8F, G).

The expression of the NAT10 gene in in tumor and normal tissues, and the results are shown in (Figure 9A). There was a significant increase in NAT10 expression in the Capan-1 cell line in comparison to normal pancreatic ductal epithelial cells (P<0.05) (Figures 9B, C). The role of NAT10 in tumorigenicity was studied more thoroughly by transfecting Capan-1cells with a specific siRNA targeting it. We designed two siRNAs to knock down NAT10. Western blot analysis confirmed the effectiveness of NAT10 knockdown, and both siRNAs significantly knocked down NAT10, and the protein level of NAT10 were barely detectable after transfection with siRNA#2 (Figure 9D). The data demonstrate that the siRNA targeting NAT10 has high specificity and transfection efficiency. In both siRNA-transfected cells, cell proliferation was inhibited with increasing time. According to the Capan-1 cell line colony formation assay, knockdown of NAT10 led to a reduction in colonies (Figure 9E). By knocking down NAT10, pancreatic ductal adenocarcinoma cells are inhibited from proliferating and forming colonies. Transwell analysis also revealed that specific siRNA inhibited migration of cells knockdown of NAT10 (Figure 9F). In addition, the siRNA targeting NAT10 reduced the IC50 of gemcitabine in Capan-1 cells from 76.28μM to 23.34μM (Figure 9G), indicating that NAT10 knockdown significantly reduced the resistance of pancreatic ductal adenocarcinoma cells to gemcitabine.