We obtained expression difference of PCIF1 between tumor and normal tissues based on Cancer Genome Atlas (TCGA) database. For those tumors that lack data for normal tissues, we used genotype-tissue expression (GTEx) database as a Supplementary Material . We also obtained violin plots of the PCIF1 expression difference in different pathological stages of TCGA tumors via the “Pathological Stage Plot” module of GEPIA2 (http://gepia2.cancer-pku.cn/) (15). The log2 (transcripts per million (TPM) +1) were applied as log-scale for violin plots. We also examined specificity of PCIF1 RNA expression in different tissues, blood cell lineages, and single cell types via the consensus dataset composed of Human Protein Atlas (HPA), GTEx, and FANTOM5 datasets on the Human Protein Atlas website (https://www.proteinatlas.org/). We applied p-value as the measurable indicator for statistical hypothesis during the whole study and suggested a significant difference for a result of p < 0.01, a statistical difference for a result of p < 0.05, and no statistical difference for p > 0.05 herein and in the rest of the content (16).

We evaluated the correlation of PCIF1 RNA expression and overall survival (OS), disease-free survival (DFS/RFS) via “Survival Analysis” module of GEPIA2 (15). Cutoff-high (50%) and cutoff-low (50%) values were applied as the expression thresholds.

We utilized the cBioPortal web (https://www.cbioportal.org/) to explore genetic alteration characteristics of SND1 by pressing “quick search” and entering PCIF1 (17, 18). The expected results including a three-dimensional structure diagram of PCIF1, alteration frequency, mutation type, copy number alteration (CNA), and PTM sites via “Cancer Types Summary”, “Plots”, and “Mutations” module. We also obtained the differences of OS, DFS, progression-free survival (PFS), and disease-specific survival (DSS) for pan-cancer with or without PCIF1 alteration via the “Comparison/Survival” module and generated Kaplan-Meier plots with log-rank p-value.

The protein expression and phosphoprotein level (S30, S116, S140, S144, S159, T150, T158) of PCIF1 was obtained from the UALCAN website (http://ualcan.path.uab.edu/analysis-prot.html) on the basis of clinical proteomic tumor analysis consortium (CPTAC) dataset (19). In UALCAN, only data for six tumors were available including ovarian cancer, colon cancer, clear-cell renal cell carcinoma (RCC), uterine corpus endometrial carcinoma (UCEC), breast cancer, and lung adenocarcinoma (LUAD).

The TIMER (http://timer.cistrome.org), CIBERSORT (https://cibersort.stanford.edu/), EPIC (https://gfellerlab.shinyapps.io/EPIC_1-1), MCP-Counter (http://github.com/ebecht/MCPcounter), quanTIseq (http://icbi.i-med.ac.at/software/quantiseq/doc/index.html), and xCell (http://xcell.ucsf.edu/) were applied for evaluation of PCIF1 immune infiltration in 33 TCGA tumors. The purity and infiltration level of CD8⁺, CD4⁺, neutrophil, macrophage, and B cell with log2 TPM as scale were illustrated by the “Immune Estimation” module of TIMER2.0. The Rho and p-values were obtained from Spearman’s correlation test.

We tested the correlation of PCIF1 expression and eight immune checkpoint genes (ICGs) including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, SIGLEC15, and TIGIT.

Additionally, we provided Spearman’s correlation analysis of tumor mutational burden (TMB), microsatellite instability (MSI), and PCIF1 expression. The data were acquired from TCGA database, downloaded from the Genomic Data Commons (GDC) data portal website (https://portal.gdc.cancer.gov/) and managed with R software v4.0.3 for statistical analysis. For reliable immune score evaluation, we use immunedeconv, an R package that integrates six latest algorithms, including TIMER, CIBERSORT, EPIC, MCP-counter, quanTIseq, and xCell (https://grst.github.io/immunedeconv) (20).

To figure out interplays of PCIF1 protein, we carried out a series of gene-enrichment analysis. We first obtained top 100 PCIF1-binding proteins in Homo sapiens via the STRING database (http://string-db.org). Parameters we used were “evidence” for meaning of network edges (“evidence”), “low confluence (0.150)” for minimum required interaction score, “100” for max interactors, and no limitations/”experiment” for active interaction sources to get predicted functional partners of PCIF1. We used “Similar Gene Detection” module of GEPIA2 get top 100 PCIF1-correlated genes. Top 5 correlated genes were selected for further analysis. We applied the “Correlation Analysis” module of GEPIA2 to Pearson’s test of the top 5 PCIF1-correlated genes. Venn diagram was provided to obtain an intersection of 100 PCIF1 correlated genes from GEPIA2 and 100 interacted proteins from STRING with Evenn (http://www.ehbio.com/test/venn/). We also conducted global percentage and relation network diagram to depict pathway activity of PCIF1, its interacted proteins, and correlated genes via GSCALite, a powerful and widely applied tool for analyzing cancer multiomics and drug sensitivity (http://bioinfo.life.hust.edu.cn/web/GSCALite/) (21).

In our study, we first assessed differential expression of PCIF1 among tumor and adjacent normal tissues in 28 types of tumors based on data collected from the TCGA database and GTEx database as Supplementary Materials . Expression of PCIF1 had a significant difference in 25 of 28 tumor types like adrenocortical carcinoma (ACC), breast invasive carcinoma (BRCA), colon adenocarcinoma (COAD), glioblastoma multiforme (GBM), kidney chromophobe (KICH), and prostate adenocarcinoma (PRAD) (p < 0.01), while no statistic difference in the rest of tumor types (p > 0.05) ( Figures 1A, B ). Violin plots showed that the difference expression of PCIF1 had statistical significance in each pathological stages of COAD, KICH, kidney renal clear cell carcinoma (KIRC), liver hepatocellular carcinoma (LIHC), pancreatic adenocarcinoma (PAAD), and skin cutaneous melanoma (SKCM) ( Figure 1C ). Additionally, we conducted total protein expression distribution of the PCIF1 among primary tumor and normal tissues from clinical proteomic tumor analysis consortium (CPTAC) dataset. Six available tumors were chosen including breast cancer, clear cell renal cell carcinoma (RCC), colon cancer, ovarian cancer (OV), lung adenocarcinoma (LUAD), and uterine corpus endometrial carcinoma (UCEC). The PCIF1 protein expression in UCEC and LUAD was much more higher in normal tissue than primary tumor (p < 0.01) on the contrary to that in clear cell RCC (p < 0.01). There was no significant difference in the PCIF1 protein expression between normal and tumor tissues in the rest of the three cancers including breast cancer, colon cancer, and OV (p > 0.05) ( Figure 1D ). PCIF1 RNA specificity was low in different tissues, blood cell lineages, single cell types ( Figure 1E ), and cell lines, which was not enclosed.

The OS and RFS significance map of those who had statistical difference, namely, ACC, KIRC, kidney renal papillary cell carcinoma (KIRP), mesothelioma (MESO), SKCM, and thyroid carcinoma (THCA) were demonstrated ( Figure 2 ). As shown in Figure 2A , the high expression of PCIF1 was correlated with higher OS of ACC (p = 0.031), KIRC (p = 0.0079), KIRP (p = 0.02), and MESO (p = 0.014), while only high expression of PCIF1 was correlated with a better RFS (p = 0.023) in KIRC. Highly expressed PCIF1 was positively linked with RFS of THCA (p = 0.038) but had no statistic difference in OS of THCA (p > 0.05). Additionally, we have analyzed the correlation of PCIF1 expression and PFS and DSS data on TCGA datasets ( Supplementary Figure S1 ). Interestingly, there is a correlation between highly expressed PCIF1 and poorer DSS of cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC) (p = 0.0012) and COAD (p = 0.0482). The above data suggested that PCIF1 expression was differentially associated with the prognosis of specific tumor types and mostly presented a positive correlation with OS and RFS.

Furthermore, genetic alteration status of PCIF1 were tested and shown ( Figure 3 ). COAD possessed the highest gene alteration frequency of PCIF1 (9.43%), with 11 of 594 (1.85%) “mutation” cases and 45 of 594 (7.58%) “amplification” cases. The highest mutation frequency (3.07%) of PCIF1 appeared in endometrial carcinoma. There was a fairly low frequency of structural variant in esophagogastric adenocarcinoma, nonsmall-cell lung cancer (NSCLC) and deep deletion in leukemia, NSCLC, and prostate adenocarcinoma ( Figures 3A, B ). The total somatic mutation frequency was 0.9%, in which missense mutations were in the majority. R265H/C was detected in two cases of UCEC and one case of COAD ( Figure 3C ), while the role of amino acid motifs remained unknown. However, altered PCIF1 had no correlation with OS, RFS, PFS, and DSS of pan-cancer and each tumor type (p > 0.05) ( Figure 3E ).

We also analyzed differences in phosphorylation of PCIF1 among normal tissues and tumor lesions (OV, colon cancer, clear cell RCC, UCEC, breast cancer, and lung adenocarcinoma) based on CPTAC dataset ( Figure 4 ). PCIF1 phosphorylation level in OV had no significant difference among normal tissues and tumor sites ( Figure 4A , p > 0.05). The S116, S159, and S144 loci of PCIF1 possessed a higher phosphorylation level in colon cancer tissues, but the phosphorylation level at the T150 locus was higher in normal tissues with statistical difference ( Figure 4B , p < 0.05). An increased phosphorylation level was also seen in clear cell RCC and breast cancer tissues, with statistical difference, compared with normal tissues ( Figures 4C, E , p < 0.05) while the situation was opposite in UCEC ( Figure 4D , p < 0.05). PCIF1 phosphorylation of the T158 locus was elevated in normal tissues compared with lung adenocarcinoma tissues; PCIF1 phosphorylation of the S116 locus had no statistical difference between the two groups. PCIF1 was predicted to be located intracellularly (nucleoplasm, microtubules, cytokinetic bridge). Whether the phosphorylation of PCIF1 has an effect on its location or its function (e.g., catalytic activity) remains unknown, requiring more investigations.

We applied a series of analytical tool on TIMER2.0 including TIMER, CIBERSORT, EPIC, MCP-counter, quanTIseq, and xCell to investigate immune infiltration of PCIF1 ( Figure 5 and Supplementary Figure S2 ). Figure 5 displays the expression level of PCIF1 in CD8+ T cell, CD4+ T cell, neutrophil, myeloid dendritic cell (MDC), macrophage, and B cell in 33 tumor types, among which PCIF1 expression harbored a significant difference in immune cells in liver hepatocellular carcinoma (LIHC) and prostate adenocarcinoma (PRAD). The xCell had revealed that the significant difference of PCIF1 expression in monocytes and macrophages (“macrophage,” “macrophage M1,” and “macrophage M2”) of GBM and sarcoma (SARC). Furthermore, we created linear regressions in ACC, KIRC, KIRP, and THCA to see how two variables fit. The selected tumors are those who had statistical difference not only in expression among normal tissue and tumor sites but also in OS or RFS ( Figure 5 and Supplementary Figure S3 ). Six immune cell types were chosen according to xCell. Figure 5 only included KIRC and THCA, since there was no significant difference of PCIF1 immune infiltration in ACC and KIRP (p > 0.05). PCIF1 expression in CD4+ T cell, neutrophil, MDC and B cell of KIRC, and CD8+ T cell and macrophage of THCA was positively correlated with their infiltration level, while negatively correlated in CD4+ T cell of THCA.

Immune microenvironment influences therapeutic response of immune checkpoint blockade. Expression profile of ICGs is linked with immune infiltration, response of immunotherapy, and prognosis. Therefore, we analyzed eight classical ICGs including CD274 (PD-L1), CTLA4 (CD152), HAVCR2 (TIM3), LAG3, PDCD1 (PD-1), PDCD1LG2, SIGLEC15, and TIGIT that have an either positive or negative correlations with PCIF1 expression in all 33 TCGA tumors (22) ( Figure 6 ). Figure 6 suggests that PCIF1 expression had a strong correlation with almost all the selected immune checkpoint genes in LIHC, thymoma (THYM), THCA, PRAD, and COAD, consistent with data of PCIF1 immune infiltration. No relationships were found in uveal melanoma (UVM), KICH, and esophageal carcinoma (ESCA).

As TMB (23) and MSI (24) are considered biomarkers for immunotherapy, we also depicted the correlation between PCIF1 expression and TMB or MSI with bubble diagram to further evaluate its immunogenicity ( Figure 7 ). The expression of PCIF1 was positively correlated with TMB of MESO, GBM, LUAD, and brain lower-grade glioma (LGG) but negatively correlated with TMB of COAD, THYM, PRAD, THCA, BRCA, and UCEC (p < 0.05). The expression of PCIF1 was positively correlated with MSI of KICH, KIRC, LIHC, and lung squamous cell carcinoma (LUSC) but negatively correlated with MSI of rectum adenocarcinoma (READ) (p < 0.05). Interestingly, there were no overlaps between TMB and MSI groups.

A few progresses have been made in 2019 revealing the role of PCIF1 role as the only known m6Am methyltransferase, but still blank to cancer so far. Thus, we carried out a series of enrichment analysis of PCIF1 to explore its oncogenic role via STRING, GEPIA2, and GSCALite (21). Meanwhile, we gained a list of top 100 genes that have correlations with PCIF1 via the module of “similar gene detection” of GEPIA2 ( Supplementary Table S1 ). Thirty-four proteins had an interaction with PCIF1 in Homo sapiens based on the STRING database; eight interacted proteins were identified by experiments, namely, PDX1, MNAT1, HDAC3, POLR2A, POLR2I, CAMSAP2, C19of43, and FANCM ( Figures 8A, B and Supplementary Table S2 ). Pearson’s test displayed correlation of PCIF1 expression and top 5 correlated genes including ZSWIM1 (R = 0.55, p < 0.01), UBOX5 (R = 0.5, p < 0.01), SLC35C2 (R = 0.5, p < 0.01), IFT52 (R = 0.48, p < 0.01), and EHMT2 (R = 0.47, p < 0.01) ( Figure 8C ). We also obtained global percentage of PCIF1, the abovementioned eight PCIF1-binding proteins (PDX1, MNAT1, HDAC3, POLR2A, POLR2I, CAMSAP2, C19orf43, FANCM), and top 5 correlated genes (ZSWIM1, UBOX5, SLC35C2, IFT52, EHMT2) of TCGA tumors, as well as interaction map of gene and pathway ( Figures 9A, B ). PCIF1 had an activation on DNA damage response and hormone androgen receptor (AR), while an inhibition on apoptosis ( Figure 9A ). However, the Venn diagram suggested that there was no overlap between the interacted and correlated proteins ( Figure 9 ).

Furthermore, we obtained heatplots of negative or positive differentially expressed genes and Spearman’s correlation test of PCIF1-binding proteins in ACC, KIRC, KIRP, and TYHA ( Figures 10A, B ). These four tumors were those with statistical difference in normal and tumor tissues and OS and RFS as mentioned before. We defined absolute statistic value (r value): strong correlation for 0.7–1.0, moderate correlation for 0.4–0.7, weak correlation for 0.2–0.4, and extreme weak correlation for 0.0–0.2 (25). A positive r value indicated a positive correlation and vice versa. We found that in ACC, PCIF1 had a positive moderate correlation with PHF20 (r = 0.649), C20orf43 (r = 0.640), and CHMP4B (r = 0.638), while a negative moderate correlation with LOC90784 (r = −0.608), DNM1 (r = −0.586), and COPS7B (r = −0.581). In KIRC, PCIF1 had a positive moderate correlation with DVL2 (r = 0.536), XAB2 (r = 0.526), and UNK (r = 0.486), while a negative moderate correlation with ACER3 (r = −0.481), SEC24A (r = −0.463), and HSPA4 (r = −0.459) (p < 0.05). In KIRP, PCIF1 had a positive moderate correlation with CLDN3 (r = 0.549), SPEF1 (r = 0.516), and C17orf107 (r = 0.499), while a negative moderate correlation with FAM40B (r = −0.460), CHAC2 (r = −0.456), ATG4A (r = −0.442) (p < 0.05). In TYHA, PCIF1 had a positive moderate correlation with TBC1D17 (r = 0.594), BCKDHA (r = 0.579), and TUBGCP2 (r = 0.549), while a negative moderate correlation with UBA6 (r = −0.565), CAPZA1 (r = −0.543), and BZW1 (r = −0.518) (p < 0.05). KEGG pathway analysis revealed that PCIF1 was positively correlated with pathways associated with ferroptosis, drug metabolism, necroptosis in ACC, spliceosome, Notch signaling pathway in KIRC, ribosome, and metabolism (e.g., carbon, glycine, serine, threonine) in TYHA while negatively correlated with pathways associated with ribosome, DNA replication, cell cycle and mismatch repair in ACC, protein export, ferroptosis, proteasome in KIRC, natural killer cell-mediated cytotoxicity, antigen processing and presentation in KIRP, TNF, JAK-STAT, and NOD-like receptor signaling pathway in TYHA. However, the interaction mechanisms required more investigations.