The PEX13 mRNA expression data for 33 tumors retrieved from The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) database were analyzed using the SangerBox portal (http://SangerBox.com/Tool) (18). These tumors included adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large b-cell lymphoma (DLBC), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck cancer (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney papillary cell carcinoma (KIRP), acute myeloid leukemia (LAML), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian cancer (OV), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarcinoma (PRAD), rectum adenocarcinoma (READ), sarcoma (SARC), skin cutaneous melanoma (SKCM), stomach adenocarcinoma (STAD), testicular germ cell tumors (TGCT), thyroid carcinoma (THCA), thymoma (THYM), uterine corpus endometrial carcinoma (UCEC), uterine carcinosarcoma (UCS) and uveal melanoma (UVM). PEX13 mRNA expression distributions were visualized in violin plots. Subsequently, the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) portal (genepattern.org) was applied to further explore the PEX13 mRNA expression in 33 tumors (19). The immunohistochemical (IHC) staining data for the analysis of PEX13 expression in various tumor tissues and corresponding normal tissues were obtained from The Human Protein Atlas (THPA) portal (https://www.proteinatlas.org/) (20). The analyzed tumors included glioma and breast, colorectal, liver, lung, pancreatic, skin, stomach, thyroid and renal cancer. Survival prognosis analysis was performed using GEPIA2 and SangerBox. Through GEPIA2, the overall survival (OS) map data for PEX13 in multiple tumors of TCGA database was obtained. The low- and high-expression cohorts of PEX13 were obtained using the expression thresholds of cutoff-high (>50%) and cutoff-low (<50%) values. Cox regression analysis was used to analyze the OS, disease-free interval (DFI), progression-free interval (PFI) and disease-specific survival (DSS) data in various cancers using SangerBox portal.

The characteristics of PEX13 genetic alterations were explored in multiple tumors using the cBioPortal database (https://www.cbioportal.org/), including mutation, amplification and deep deletion (21). The structural variant data, mutation data and copy number alteration (CNA) data were obtained from the cBioPortal database. Data of somatic CNAs and somatic mutations of PAAD were obtained from TCGA datasets (gdc.cancer.gov/). The CNAs were related to the PEX13 expression and the threshold CNA peaks were analyzed through GISTIC 2.0 (https://cloud.genepattern.org/) (22). According to the PEX13 expression level, the patients with PAAD were divided into the first 25% PEX13low (n=46) and the last 25% PEX13high (n=46) groups. To download and visualize the somatic mutations of patients with PAAD, the ‘maftools’ package (bioconductor.org/packages/release/bioc/vignettes/maftools/inst/doc/maftools.html#7_Visualization) was applied in R software (4.2.1) (https://www.r-project.org/) (23).

Through the GEPIA2 portal, the top 200 PEX13-related genes were explored in TCGA database. Subsequently, the TIMER2 (http://timer.cistrome.org/) portal was used to generate a heatmap of the top 10 PEX13-related genes in various cancers, including KCMF1, KIAA1841, MPHOSPH10, MRPL19, PPP1CB, PUS10, RAB1A, SLC30A6, SMEK2 and YIPF4 (24). The correlation coefficient (R) and P-value were calculated. The Gene Ontology (GO)-biological process (BP) subontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 200 PEX13-correlated genes in pan-cancers were explored via DAVID (https://david.ncifcrf.gov/home.jsp) and SangerBox portal. Moreover, the top 50 PEX13 negatively related genes and top 50 PEX13 positively related genes were explored and heatmaps were generated in PAAD using LinkedOmics (www.linkedomics.org/login.php) (25). The GO-BP and KEGG enrichment analysis of the top 50 PEX13 negatively related genes and top 50 PEX13 positively related genes in PAAD were explored via DAVID and SangerBox portal. Linkedomics (linkedomics.org) portal was used to analyze the enrichment pathways of PEX13 positively and negatively related genes in PAAD (25).

The hub genes affected by PEX13 expression level were obtained through WGCNA using SangerBox portal. First, the PAAD gene expression profile in TCGA database was used to calculate the median absolute deviation (MAD) of each gene. The outlying genes and samples were removed using the GoodSamplesGenes method of the ‘WGCNA’ R software package (genetics.ucla.edu/labs/horvath/CoexpressionNetwork/Rpackages/WGCNA) (26). WGCNA was then applied to build a scale-free co-expression network. As a soft-thresholding parameter, β could emphasize strong gene-gene associations and penalize weak correlations. To ensure a scale-free network, the power of β=7 (0.85) was selected as the soft-thresholding parameter and the power of β=7 (16.23) was selected as the soft-thresholding parameter in mean connectivity. Moreover, the PAAD samples were clustered using Pearson's correlation and average linkage method. The average linkage hierarchical clustering identified 23 modules and the differences among these were analyzed according to the expression level of PEX13, among which darkgrey module had the strongest association with PEX13 expression. Consequently, the darkgrey module was chosen as the important module for extracting hub genes.

Through the STRING database (https://string-db.org/), 50 PEX13 interacting proteins were identified in pan-cancers. Using Cytoscape (3.10.1; cytoscape.org/download.html), a PEX13 interacting protein network map was generated. The key protein-protein interaction network of PEX13 was extracted in this protein-protein interaction network using Cytoscape. Molecular docking technique was applied to predict the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3. By using the HDOCK server (http://hdock.phys.hust.edu.cn/), the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 were analyzed. Using the PyMOL tool (2.5.5; pymol.org/2/), the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 were visualized (27).

In addition, lncRNA/circRNA-miRNA-PEX13 networks were investigated. Through the niRabel (http://bioinfo.univ-rouen.fr/mirabel/index.php), miRDB (http://mirdb.org/), TargetScan (http://www.targetscan.org/vert_72/) and ENCORI (https://starbase.sysu.edu.cn/) databases, 16 miRNAs that might target PEX13 mRNA were identified. Subsequently, four miRNAs were further screened and selected through miRabel, miRDB and miRWalk scores, including hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p. Using the ENCORI online tool, 105 lncRNAs and 96 circRNAs that targeted hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p were predicted, and the lncRNA-miRNA-PEX13 and circRNA-miRNA-PEX13 regulatory networks were visualized via Cytoscape tool.

The SangerBox portal was used to investigate the relationships between PEX13 mRNA expression and tumor infiltration immune cells (TIICs) in various cancers, including T_cells_CD8, T_cells_CD4_memory_resting, B_cell_memory, Tregs, NK_cells_activated, macrophages and dendritic_cells_activated. Subsequently, the correlations between PEX13 mRNA expression and immune checkpoint (ICP) genes, microsatellite instability (MSI), tumor mutation burden (TMB) and tumor purity in multiple TCGA tumors we investigated using SangerBox. Furthermore, TIMER2 portal was applied to analyze the associations between PEX13 and the infiltration levels of B cell and macrophage in pan-cancers, especially in PAAD via different algorithms, such as XCELL, CIBERSORT and EPIC. Moreover, the relationships between PEX13 and stromal score, immune score and ESTIMATE score in multiple cancers were explored using SangerBox, especially in lung adenocarcinoma (LUAD), stomach and esophageal carcinoma (STES), sarcoma (SARC) and lung squamous cell carcinoma (LUSC).

The human PAAD cell lines of PANC1 and BxPC3 were obtained from the Shanghai Cell Bank of the Chinese Academy of Sciences (Shanghai, China). Comprehensive cell line authentication was performed and cells were periodically checked for mycoplasma contamination. The DMEM (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% FBS (FBS; Satorius AG) was used to culture the PAAD cell lines. The PAAD cells were incubated in a cell incubator at 37°C with 5% CO₂.

siRNAs used to knock down the expression of PEX13 were obtained from Shanghai GenePharma Co., Ltd. A total of 2×10⁵ PAAD cells were seeded in 6-well plates, 4.5 µl Lipofectamine™ RNAiMAX reagent (Thermo Fisher Scientific, Inc.) and 40 pmol siRNA were mixed for 10 min in each well at room temperature. Subsequently, the mixture was added to the cells cultured with FBS-free medium and the cells were replaced with full medium at 24 h after transfection. The PAAD cells were collected after 48 h with TRIzol^® reagent (Takara Bio, Inc.). The GoScript reverse transcription system (Promega Corporation) was applied to generate PEX13 cDNA and the GoTaq^® qPCR Master Mix (Promega Corporation) was used to detect the knockdown efficiency of PEX13 siRNAs through qPCR using the ABI QuantStudio 3 system. The siRNA sequences were as follows: si-PEX13-1 sense, 5′-ACUUGAUUCCACCGUUUUCCUTT-3′ and antisense, 5′-GAAAACGGUGGAAUCAAGUAATT-3′; si-PEX13-2 sense, 5′-UUAUUCUAGUCCUGAAUUCCUTT-3′ and antisense, 5′-GAAUUCAGGACUAGAAUAAAGTT-3′; and siRNA negative control (si-NC) sense, 5′-UUCUCCGAACGUGUCACGUTT-3′ and antisense, 5′-ACGUGACACGUUCGGAGAATT-3′. The thermocycling conditions used for qPCR were: Initial denaturation at 95°C for 10 min; 40 cycles of amplification for 15 sec at 95°C and 45 sec at 60°C, and an extension step at 72°C for 2 min. The relative expression levels of the PEX13 mRNA were calculated using the 2^−ΔΔCq method and normalized to the internal control GAPDH (28). The primer sequences were as follows: GAPDH forward, 5′-GGTGGTCTCCTCTGACTTCAACA-3′ and reverse, 5′-GTTGCTGTAGCCAAATTCGTTGT-3′ reverse; PEX13 forward, 5′-GCCCCACTTTCCAATCTGCT-3′ and reverse, 5′-AGAATAGGTGGGGGCACTCT-3′.

A total of 2×10³ PANC1 or BxPC3 cells were plated into 96-well plates. PEX13 expression was knocked down using the PEX13 siRNAs. Subsequently, the optical density (OD) at 450 wavelengths was measured for 5 consecutive days (at days 0, 1, 2, 3, 4 and 5) after incubation for 1 h at 37°C with the CCK-8 reagent (ApexBio Technology) using a multiplate reader.

A total of 2.5×10⁵ PANC1 or BxPC3 cells were plated in 6-well plates. A 200 µl pipette tip was used to scratch the cell monolayer. After 0, 12 and 24 h incubation in an FBS-free medium, the images of wound closure were captured and further explored via an inverted microscope (Olympus Corporation) and Image J software (1.8.0; National Institutes of Health). For the Transwell assay, 8-µm pore filters coated with Matrigel matrix glue (Corning, Inc.) were used. Matrigel was added into the upper chamber for 40 min at 37°C. A total of 5×10⁴ PANC1 and BxPC3 cells in 300 µl serum-free DMEM were seeded into the upper chamber and 600 µl DMEM with 10% FBS medium was added to the lower chamber. After 24 h, cells were fixed with 4% paraformaldehyde at room temperature and stained with 1% crystal violet at room temperature. The images were captured using a microscope (BX53; Olympus Corporation).

In 6-well plates, 1×10³ PANC1 and BxPC3 cells were seeded in each well with 2.5 ml DMEM containing 10% FBS. The cells were incubated for 2 weeks at 37°C in a humidified chamber with 5% CO₂. The colonies were fixed with 4% paraformaldehyde for 30 min at room temperature before being stained with crystal violet for 10 min at room temperature. A cell mass of >50 cells was considered a colony, and Image J (1.8.0) was used for counting colonies.

All of the experimental data were statistically analyzed using GraphPad Prism 8.0 (Dotmatics). Comparisons between two groups were made using an unpaired t-test and the statistical method used for survival analysis was Kaplan-Meier. Cox regression analysis was performed using logrank test to obtain prognostic significance. The correlation analyses in the manuscript were applied using Pearson's correlation analysis. One-way ANOVA with Bonferroni's post-hoc test was used to compare the differences between multiple groups. Two-way ANOVA with Bonferroni's post-hoc test was used to compare the differences between CCK-8 groups (si-NC and si-PEX13-1/si-PEX13-2). The data are presented as mean ± standard deviation from three individual experiments and each experiment was repeated at least three times. P<0.05 was considered to indicate a statistically significant difference.

Through TCGA and GTEx databases analysis, PEX13 was revealed to be highly expressed in multiple tumors, including GBM, LGG, GBMLGG, BRCA, UCEC, LUAD, CESC, ESCA, COAD, STES, READ, STAD, PRAD, HNSC, LUSC, LIHC, high-risk wilms tumor (WT), BLCA, THCA, PAAD, UCS, LAML, acute lymphoblastic leukemia (ALL) and ACC. In addition, the expression level of PEX13 was downregulated in KIRP, pan-kidney cohort (KIPAN), KIRC, SKCM, TGCT and KICH compared with the corresponding normal tissues (Fig. 1A). In GEPIA2 portal, it was discovered that PEX13 was remarkable upregulated in DLBC, GBM, LGG, LUSC, PAAD, STAD and THYM (Fig. 1B). In addition, the results of IHC showed that the expression of PEX13 protein in breast cancer, colorectal cancer, glioma, liver cancer, lung cancer, pancreatic cancer, skin cancer, stomach cancer, thyroid cancer and renal cancer was higher than that in normal tissues (Fig. S1). The present findings indicated that PEX13 mRNA and protein expression levels are increased in a variety of tumors, suggesting that PEX13 may function as an oncogene in various cancers, including PAAD.

The prognostic value of PEX13 in various tumors was analyzed by comparing the OS of patients with PEX13 expression. Analysis of the pan-cancer cohort using the GEPIA2 database demonstrated that high PEX13 mRNA expression levels was associated with a significantly lower OS in patients with ACC, KICH, LGG, LIHC and PAAD, and high PEX13 mRNA expression levels were associated with a significantly higher OS in patients with KIRC (Fig. 2A). Cox regression analysis in SangerBox database suggested that PEX13 mRNA expression level was related to OS, DFI, DSS and PFI in patients with various tumors. Regarding the OS, the findings revealed that high PEX13 mRNA expression was linked to a shorter OS in GBMLGG, LGG, PAAD, ACC, KICH, BRCA, LIHC and LAML (Fig. 2B). Regarding the DFI, high PEX13 mRNA expression was linked to a shorter DFI in PAAD, KIPAN, KIRP and ACC (Fig. S2A). In PFI, the high PEX13 expression was related to shorter PFI in GNMLGG, ACC, LGG, PAAD, LIHC, KICH and UVM (Fig. S2B). Furthermore, a high PEX13 mRNA expression was correlated with a shorter DSS in GBMLGG, LGG, PAAD, ACC, KICH and LIHC (Fig. S2C). It was also observed that high PEX13 mRNA expression was significantly correlated with shorter OS, DFI, PFI, and DSS in PAAD and ACC. The present results revealed that the level of PEX13 mRNA expression was related to the prognosis of various tumors. Furthermore, high PEX13 mRNA expression was linked to a worse prognosis in patients with PAAD.

Mutations, deletions or amplification of oncogenes or tumor suppressor genes have been linked to the growth and progression of several tumors (29). Hence, the present study analyzed diverse types of genetic alterations in the PEX13 gene using the cBioPortal portal, including amplification, mutation and deep deletion. The most common genetic alteration of PEX13 gene in DLBC (18.75%), LUSC (5.75%), BLCA (2.68%), UCS (1.75%), ovarian serous cystadenocarcinoma (1.54%) and ESCA (1.1%) was amplification. In UCEC (2.46%), STAD (1.36%), COADREAD (1.35%) and KIRP (0.71%), the most common genetic alteration of PEX13 gene was mutation (Fig. 3A). Subsequently, the TCGA-PAAD dataset was interrogated to investigate the correlation between PEX13 and specific genomic features such as copy number variations (CNVs) and somatic mutations. Fig. 3B shows the comparison of the CNV profiles in the PEX13high (n=46) and PEX13low (n=43) groups. In the PEX13high group, deletion peaks were found in the 5p12.3, 6p25.3, 9p21.3, 10p15.3, 10q23.2, 18q12.1, 18q21.2 and 19p13.2 chromosomal locations, while amplification peaks were found in the 1p12, 6p21.31, 7q21.13, 8q24.21, 9p13.3, 12p13.33, 18q11.2 and 19q13.2 chromosomal regions (Fig. 3C-D). In the PEX13low group, deletion peaks were found in the 9p23 and 9p21.3 chromosomal locations, while frequent amplification peaks were found in the 8p11.23 and 19q13.2 chromosomal regions (Fig. 3C and D). These findings suggested that high PEX13 expression led to more deletions and amplifications of multiple sites in the genome of patients with PAAD. The PEX13low group indicated that high frequency of somatic mutations in the FRG1B (72%), KRAS (49%), TP53 (47%), UBBP4 (44%) and TTN (40%) genes and the PEX13high group suggested that high frequency of mutations in the KRAS (89%), FRG1B (76%), TP53 (72%), MAMLD1 (35%), MED12L (35%) and UBBP4 (30%) genes (Fig. 3E and F). These findings revealed that high PEX13 expression may lead to the increased mutation probability of KRAS, FRG1B, TP53 and other genes, and to the decreased mutation probability of UBBP4 gene in PAAD, which may be a potential mechanism leading to the occurrence and development of tumors.

These findings suggested the presence of PEX13 gene mutations, amplification and deletion in multiple tumors. The PAAD tissues revealed different somatic mutations and CNVs depending on the level of PEX13 expression. This suggested that alterations in the PEX13 gene could influence the occurrence and progression of multiple cancers, particularly PAAD.

The PEX13-related genes were analyzed in 33 cancers through GEPIA2 database and 200 PEX13-related genes were obtained. The top 10 genes with the highest correlation among 200 PEX13-related genes were generated correlation heat maps in 33 cancers via TIMER2, including KCMF1, KIAA1841, MPHOSPH10, MRPL19, PPP1CB, PUS10, RAB1A, SLC30A6, SMEK2 and YIPF4 (Fig. 4A). Fig. 4B shows that PEX13 expression was positively correlated with PUS10 (R=0.61), KIAA1841 (R=0.61), RAB1A (R=0.59), SMEK2 (R=0.58), YIPF4 (R=0.57), PPP1CB (R=0.57), KCMF1 (R=0.57), SLC30A6 (R=0.57), MPHOSPH10 (R=0.56) and MRPL19 (R=0.56) in GEPIA2. The GO-BP analysis of 200 PEX13-related genes in 33 tumors showed that PEX13-related genes were mainly enriched in mRNA metabolic process, cellular component biogenesis, cell division, intracellular transport, cytokinesis, mRNA processing, cell cycle and RNA splicing (Fig. 4C). PEX13-related genes were mainly enriched in Endocytosis, mRNA surveillance pathway, nucleocytoplasmic transport, T cell receptor, signaling pathway and hippo signaling pathway based on the results of the KEGG enrichment analysis (Fig. 4D).

To further explore the biological functions of PEX13 in PAAD, the enrichment of PEX13 positively and negatively related genes and their enrichment pathways in PAAD were analyzed using the Linkedomics portal. Fig. S3A shows the heat maps of 50 positively and 50 negatively related genes of PEX13 in PAAD. Fig. S3B-E showed the GO-BP and KEGG enrichment analysis of 100 positively and 100 negatively related genes of PEX13, respectively. The results showed that GO-BP was mainly enriched in positive regulation of drug response, small GTPase mediated signal transduction, regulation of epidermal growth factor-activated receptor activity and ERBB signaling pathway among the positively related genes of PEX13 and in the regulation of glycogen (starch) synthase activity, intracellular protein transport, intracellular protein transport and autophagosome maturation among the negatively related genes of PEX13 (Fig. S3B and C). For the KEGG enrichment analysis, PEX13-positive related genes were mainly enriched in Rap1 signaling pathway, ErbB signaling pathway, EGFR tyrosine kinase inhibitor resistance, Neurotrophin signaling pathway and AMPK signaling pathway and PEX13-negatively related genes were mainly enriched in Vitamin B6 metabolism, Apelin and Insulin signaling pathway (Fig. S3D-E). These findings show that PEX13 could regulate mRNA processing, RNA splicing, cell cycle and the response to drugs in cancers. In PAAD, PEX13 expression was related to protein transport, tumor cell response to drugs, multiple cancer pathways and immune pathways.

The TCGA-PAAD cohort was analyzed via WGCNA to identify hub genes strongly associated with PEX13 expression in PAAD. Using the average linkage approach and Pearson's correlation method, the samples of TCGA-PAAD data were first clustered (Fig. 5A-C). With the help of the average linkage hierarchical clustering, a total of 23 modules were obtained (Fig. 5D). According to the expression level of PEX13, we analyzed the differences among 23 modules, among which darkgrey module had the strongest association with PEX13 expression (Fig. 5E-F). The darkgrey module was selected as the significant module for further analysis. Fig. 5G shows the hub genes of darkgrey module and GO-BP and KEGG enrichment analyses were used to analyze these hub genes. In GO-BP analysis, hub genes mainly enriched in cell migration, regulation of ERK1 and ERK2 cascade, wound healing, ephrin receptor signaling pathway, ERBB2-EGFR signaling pathway, positive regulation of MAP kinase activity and signal transduction (Fig. 5H). In KEGG enrichment analysis, hub genes mainly enriched in PI3K-Akt signaling pathway, extracellular matrix (ECM)-receptor interaction, focal adhesion, pathways in cancer and MAPK signaling pathway (Fig. 5I). These results suggested that the hub genes of darkgrey module with different PEX13 expression were mainly involved in cell migration, wound healing, PI3K-Akt and MAPK signaling pathway.

The 50 PEX13 interacting proteins were identified using STRING database and the PEX13 interacting protein network map was generated in Cytoscape (Fig. 6A). In this protein interaction network, key protein interaction network of PEX13 was extracted using Cytoscape, as shown in Fig. 6B. The PEX13 interacting proteins were scored in Cytoscape and a larger volume indicated higher interacting score of these proteins in the network diagram. Subsequently, four proteins had the highest interaction scores with PEX13, including PEX2, PEX12, PEX14 and ABCD3 (Fig. 6B). Molecular docking technique was applied to forecast the binding sites between PEX13 and PEX2, PEX12, PEX14 and ABCD3 (Fig. 6C-F). The results showed that PEX13 may interact with PEX2, PEX12, PEX14 and ABCD3 to influence the biological functions and acted as a significant role in various tumors.

In addition, numerous studies found that lncRNA and circRNA regulated the expression of downstream mRNA by adsorbing target miRNAs, which acted as a crucial factor in tumorigenesis and the progression of tumors (30,31). Therefore, the present study explored the lncRNA/circRNA-miRNA network that may regulate mRNA expression in multiple cancers. LncRNA-miRNA regulatory networks of PEX13 were constructed to investigate the underlying molecular mechanisms of PEX13 in tumors. Through the miRabel, miRDB, TargetScan and ENCORI databases, 16 miRNAs that might target PEX13 mRNA were identified (Fig. 7A). Subsequently, four miRNAs, including hsa-miR-4465, hsa-miR-16-5p, hsa-miR-195-5p and hsa-miR-6838-5p, were screened and selected using miRabel, miRDB and miRWalk. Using the ENCORI online tool, the lncRNAs and circRNAs that targeted their corresponding miRNAs were predicted and Cytoscape was used to visualize the lncRNA-miRNA and circRNA-miRNA regulatory networks (Fig. 7B and C). These findings suggested the upstream ncRNA-miRNA regulatory networks that may control the aberrant expression of PEX13 in cancers.

The relationships between PEX13 and the tumor immune microenvironment (TIM) were further analyzed in a variety of tumors, and the correlation between PEX13 and the composition of tumor infiltrating immune cells (TIICs) was investigated in multiple cancers. Studies showed that TIICs were a significant component of the tumor microenvironment and regulated the occurrence, development and metastasis of tumors (32,33). The Sangerbox database was applied to investigate the correlation between PEX13 and TIIC levels in diverse tumors. The expression of PEX13 in LAML, PRAD, BRCA, PAAD, STES, HNSC, LIHC, KIRC, THCA, STAD, SKCM, THYM and CESC was remarkably related to various immune cells, including T_cells_CD4_memory_resting, B_cell_memory, T_cells_CD8, macrophages_M1, macrophages_M2, Tregs, NK_cells_activated, dendritic_cells_activated and neutrophils (Fig. 8A). Fig. S4A demonstrated the relationship between PEX13 expression and B cell infiltration levels in different tumors. In PAAD, the expression level of PEX13 under different algorithms was mainly negatively related to the infiltration level of B cells (Fig. S4B). Fig. S4C presents a heat map of the relationship between PEX13 expression levels and macrophage infiltration levels in various tumors. PEX13 expression level was also negatively related to the infiltration level of macrophages under different algorithms in PAAD (Fig. S4D).

Furthermore, the correlation between PEX13 expression level and ICP genes, MSI, TMB and tumor purity were analyzed via SangerBox. The occurrence of MSI is caused by the functional defect of DNA mismatch repair in tumor tissues (34). The presence of MSI and ineffective DNA mismatch repair is a crucial clinical tumor marker (35). TMB has a strong relationship to the effectiveness of PD-1/PD-L1 inhibitors and can forecast their benefits in terms of efficacy (36). The present study discovered that PEX13 expression was positively linked with ICP genes in several tumor types (Fig. 8B). PEX13 expression was positively correlated with MSI in TGCT and KIRC and negatively correlated with MSI in GBMLGG, HNSC, PRAD, KIPAN, COAD, DLBC and THCA (Fig. 8C). In LUAD and BRCA and pheochromocytoma and paraganglioma, PEX13 expression was linked favorably with TMB (Fig. 8D). In tumor purity, PEX13 expression was negatively related to purity in KICH, KIPAN and PRAD and positively related to purity in GBMLGG, SARC, LUAD, LUSC, STES, LGG, STAD, BRCA, ESCA, etc. (Fig. 8E). Moreover, it was evaluated the relationship between the ESTIMATE score and PEX13 expression level in various cancers. ImmuneScore reflects the proportion of infiltrating immune cells in cancer tissues. The proportion of stromal cells in tumor tissue was reflected by the Stromal score (StromalScore). The ESTIMATE score (ESTIMATEScore) is the sum of immune and stromal scores and reflects the state of the TIM and tumor purity (18,37). The present study found significant relationships between PEX13 levels and StromalScore, ImmuneScore and ESTIMATEScore across multiple tumors, especially in LUAD, STES, SARC and LUSC (Fig. S5A-C). These findings indicated that PEX13 expression regulated the sensitivity and resistance of certain cancers to immunotherapy. Thus, PEX13 emerged from this analysis as a potential immunotherapy biomarker and predictor of tumor immunotherapy response.

PANC1 and BxPC3 cell lines were used to verify the biological functions of PEX13. The knockdown of PEX13 expression in PANC1 and BxPC3 cells transfected with PEX13 siRNAs was confirmed using RT-qPCR (Fig. 9A). The CCK-8 assay suggested that PEX13 knockdown inhibited the cell proliferation of PANC1 and BxPC3 cells in vitro (Fig. 9B). In addition, wound healing and Transwell experiments demonstrated that PEX13 knockdown suppressed the cell migration and invasion of PANC1 and BxPC3 cells in vitro (Fig. 9C-D). The colony formation assay suggested that PEX13 knockdown could inhibit the cell colony formation ability of PANC1 and BxPC3 cells in vitro (Fig. 9E). These findings showed that PEX13 could increase cell proliferation, migration, invasion and colony formation abilities of PANC1 and BxPC3 cell lines in vitro.