We compared DDX56 RNA expression among different tissues using RNA sequencing (RNA-seq) datasets from The Cancer Genome Atlas (TCGA). RNA-seq data (TPM) and related clinical data were downloaded from UCSC Xena (http://xena.ucsc.edu/) (Goldman et al., 2020). The data corresponded to 31 solid tumor types: adrenocortical carcinoma (ACC), bladder urothelial carcinoma (BLCA), breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), glioblastoma multiforme (GBM), head and neck squamous cell carcinoma (HNSC), kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), brain lower grade glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), mesothelioma (MESO), ovarian serous cystadenocarcinoma (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), uveal melanoma (UVM). Using the “Gene” module of TISCH (http://tisch.comp-genomics.org/search-gene/), we determined the RNA expression of DDX56 in multiple cell types based on single-cell RNA-seq data (Sun et al., 2021). UALCAN carried out a comparative analysis of protein expression (http://ualcan.path.uab.edu/index.html) (Chandrashekar et al., 2017). Only tumors with matched normal tissue data were used for differential analysis. We defined clinical stages as follows: early stages, TNM stages I/II; advanced stages, TNM stages III/IV.

Univariate Cox regression was used to assess the prognostic significance of DDX56 across cancer types. Multivariate Cox regression was used to identify independent prognostic factors. The surv_cutpoint function (from R package “survminer”, https://github.com/kassambara/survminer) was used to determine the optimal cutoff values of DDX56 expression level. Only solid tumor data with complete clinical information were included in the survival analysis.

Genes co-expressed with DDX56 were screened by calculating the Spearman correlation coefficient between DDX56 and all other genes in each cancer. The Benjamini–Hochberg method was used to decrease the false discovery rate (adjusted p-value). To gain functional and mechanistic insights regarding DDX56, we performed a pre-ranked gene set enrichment analysis (GSEA) (Subramanian et al., 2005) based on Hallmarker’s gene sets, which were downloaded from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/) (Liberzon et al., 2015). The genes with significant correlation coefficients (adjusted p < 0.05) were sorted according to Spearman correlation coefficient and then used in GSEA.

The CRISPR pooled libraries consist of thousands of plasmids, each containing multiple gRNAs for each target gene (Shalem et al., 2014). In a CRISPR screening experiment, target cells are treated with the pooled library to create a population of mutant cells that are then screened for a phenotype of interest. Essential genes for specific phenotypic screens can be curated by comparing sgRNA abundance. Therefore, the database includes research information on whether a particular gene is essential for a certain phenotype in a particular tumor cell line (Joung et al., 2017). Compiled CRISPR screen data were obtained from the Biological General Repository for Interaction Datasets (BioGRID) (https://orcs.thebiogrid.org/) (Oughtred et al., 2019). Data from proliferation-based CRISPR screens in solid tumor cancer cell lines were selected. These datasets were used to verify the mitogenic activity of DDX56. Resistance-related CRISPR screening evidence was also retrieved from the BioGRID database. The methods were as previously described except that the studies chosen focused on the “response to chemicals” phenotype.

Half-maximal inhibitory concentration (IC50) data and associated RNA-seq data from multiple cancer cell lines were obtained from the CellMiner database (https://discover.nci.nih.gov/cellminer/) (Reinhold et al., 2012). We conducted Spearman correlation analysis between DDX56 expression and drug IC50 in order to investigate the relationship between DDX56 expression and drug sensitivity.

Proportions of TILs were estimated using the CIBERSORT function (https://cibersort.stanford.edu/). (Newman et al., 2015; Chen et al., 2018). Based on the Spearman correlation coefficient, we assessed the statistical correlation between DDX56 expression level and the proportion of TILs.

We evaluated the value of DDX56 in predicting immunotherapy response in an immune checkpoint blockade therapy cohort. The transcriptomic data and relevant clinical data from this cohort were obtained from dbGaP (phs000452) (Liu et al., 2019a). Univariate Cox regression was used to assess the prognostic significance of DDX56.

We conducted gene mutation analysis on cBioportal (https://www.cbioportal.org/) based on TCGA Pan-Cancer Atlas data (Cerami et al., 2012). The correlation of DDX56 RNA expression level with copy number variation (CNV) was evaluated by MEXPRESS (https://mexpress.be/) (Koch et al., 2019). The correlation between the expression data and the DNA methylation data was estimated by using MEXPRESS (https://mexpress.be/). Oncodb (http://oncodb.org/) was used to compare methylation profiles between tumor tissues and normal tissues (Koch et al., 2019; Tang et al., 2022). We used the Toolkit for CistromeDB (http://dbtoolkit.cistrome.org/) to predict which transcription factors had the greatest potential to enhance expression of DDX56 (Zheng et al., 2019).

We compared non-normally distributed continuous variables using Wilcoxon rank-sum test and compared categorical variables using chi-square test between two groups. Kaplan–Meier and log-rank tests were conducted for survival analysis. Unless otherwise stated, all data analysis was performed in R (version 4.1.0). Unless otherwise specified, p < 0.05 was considered to indicate statistical significance in all analyses. The tumors involved in each analysis are recorded in Supplementary Table S1.

We compared DDX56 RNA expression between different cancer tissues and matched normal tissues based on TCGA data. The pan-cancer analysis showed that the DDX56 RNA expression level was higher in 16 solid cancers than in their matched normal tissues (Wilcoxon test, p < 0.05) (Figure 1A), consistent with findings in lung squamous cell carcinoma (LUSC), osteosarcoma (OV), and colon cancer (COAD) (Zhu et al., 2020; Pryszlak et al., 2021; Wu et al., 2021). UALCAN was used to determine the protein expression of DDX56 in different types of cancer. According to significant unique analyses (provided by UALCAN, Student’s t-test), DDX56 protein was significantly over-expressed in nine tumor types (p < 0.05, Figure 1B). An effective tumor biomarker and drug target requires specific and high expression in tumor cells compared with other components of the tumor microenvironment. Therefore, we analyzed pan-cancer single-cell RNA-seq data using the TISCH database and observed that DDX56 was mainly expressed in tumor cells rather than immune cells, stromal cells, and others (Figure 1C), indicating the potential of DDX56 as a tumor biomarker and drug target.

Next, we assessed the correlation between DDX56 expression and the clinical features of patients. Although few significant differences were observed between groups in age or gender (Supplementary Figure S1), the expression level of DDX56 was consistently higher in patients in the advanced stages than those in the early stages of diseases including ACC, BRCA, COAD, HNSC, KIRC, and LUAD (Wilcoxon test, p < 0.05, Figure 2A and Supplementary Figure S1). Further, we investigated the prognostic value of DDX56 at a pan-cancer level. The results showed that high DDX56 expression was correlated with worse outcomes in 11 cancer types (univariate Cox regression, p < 0.05, hazard ratio >1) (Figure 2B and Supplementary Figure S2). Then, we implemented multivariate Cox regression analysis based on DDX56 expression level and other clinical factors and found that DDX56 was an independent predictor for cancer prognosis in nine tumor types; moreover, worse clinical outcomes in patients with higher expression of DDX56 were observed in eight tumor types, comprising ACC, HNSC, KICH, KIRC, LIHC, LUAD, THCA, and UVM (Figure 2C). These results suggest that DDX56 high expression may be a significant predictor of prognosis and function as a promotor in multiple tumor types.

To explore the molecular mechanism of DDX56 in carcinogenesis, a network of genes co-expressed with DDX56 was first conducted in each tumor dataset, and then pathway enrichment analysis was performed using GSEA. The results revealed that DDX56 mainly participated in cell-proliferation-related signaling pathways, including G2/M checkpoint, MYC targets v1, MYC targets v2, and E2F targets (Figure 3A) (Liberzon et al., 2015). Concomitantly, we observed enrichment of tumor metabolism-related pathways including oxidative phosphorylation and unfolded protein response. In order to validate the possible pro-proliferation function of DDX56, we collected and analyzed CRISPR screens data from the BioGRID Open Respository of CRISPR Screens. CRISPR-based genetic screening is a powerful tool to identify the genes required for specific functions such as cell viability and chemical resistance. We found 512 CRISPR screens studies in 19 tumor types, which confirmed the effects of DDX56 on the proliferation of tumor cells at the pan-cancer level (Figure 3B and Supplementary Table S2).

As the GSEA results had shown enrichment of unfolded protein response (Figure 3A), a chemoresistance-related pathway, (Bahar et al., 2019; Sisinni et al., 2019), we speculated that DDX56 was associated with chemotherapy resistance. To investigate whether DDX56 was relevant to drug resistance, we searched the CRISPR screens evidence in the BioGRID database. We found that DDX56 was an essential gene required for drug resistance in five different tumors (Figure 4A and Supplementary Table S3). Moreover, we obtained IC50 levels of all available drugs for tumor cell lines from the CellMiner database and calculated the Spearman correlation coefficient between the IC50 values and DDX56 RNA expression in matched tumor cell lines. Our results revealed that high expression of DDX56 was correlated with increased resistance to cladribine, entinostat, amuvatinib, triethylenemelamine, irofulven, IWR-1, JZL-195, XL-147, and Cpd-401 (p < 0.05, Figure 4B). On the other hand, high expression of DDX56 could lead to enhanced sensitivity to a geldanamycin analog and alectinib (p < 0.05, Figure 4B).

The GSEA results suggested that DDX56 may have a critical role in inhibiting cell apoptosis and suppressing antitumor immunity (IL-2 STAT5 signaling, IL-6 JAK STAT3 signaling, and interferon gamma response, Figure 5A) (Liberzon et al., 2015). To further explore the potential relationships between DDX56 and immune cells, we examined the correlations between DDX56 and several immune cell markers including immune cells (PTPRC), T cells (CD3D, CD4, CD8A), B cells (CD19), and MHC class II molecules (Supplementary Figure S3). The results suggested that the expression of DDX56 was associated with immune infiltration in the tumor microenvironment. Next, we estimated the proportions of immune cells in each TCGA tumor type by using CIBERSORT. We observed that DDX56 expression was negatively correlated with immune infiltration levels of plasma cells, resting dendritic cells (DC), and CD4⁺ memory T cells (Figure 5B). Given that DDX56 participates in the immune infiltration process, it may have crucial biological functions in immunotherapy. Using transcriptome and clinical data from a recently published immunotherapy study (Liu et al., 2019b), we found that anti-PD-1-treated patients with high DDX56 expression had shorter progression-free survival than those with low DDX56 expression (Figure 5C).

The occurrence and development of cancer is related to gene alterations. To determine the genomic characteristics of DDX56 in cancers, comparative analysis of DDX56 was performed using cBioPortal. We observed that the main alterations of DDX56 were amplifications and mutations in multiple tumors (Figure 6A). The most common mutation type was missense mutation, which could lead to alteration of protein structure and functions (Figure 6B). As CNV can lead to higher gene mRNA levels (Haraksingh and Snyder, 2013), we analyzed the correlation between RNA expression and CNV at the pan-cancer levels. The results showed that the CNV and RNA expression of DDX56 DNA in tumor tissues were significantly associated (Pearson coefficient >0.2, p < 0.05, Figure 6C). In addition, increased RNA transcription may result from low levels of DNA methylation (Bradner et al., 2017). We compared DDX56 methylation status between various tumors and normal tissues and observed that DDX56 methylation in five tumor types was lower compared with that in paired normal tissues (Wilcoxon test, p < 0.05) (Figure 6D). Furthermore, we estimated the correlation between expression and DNA methylation data using MEXPRESS. We observed that methylation of some CpG dinucleotides and CpG islands, including cg25257687 and cg01998345, was significantly negatively correlated with the RNA expression of DDX56 (Pearson correlation coefficient <0 and p < 0.05, Supplementary Figure S4). Our results suggested that lower methylation levels of DDX56 DNA may result in high expression levels of DDX56 RNA in some cancers.

Finally, screening was performed to search candidate transcription factors possibly regulating DDX56 expression. According to our results, LYL1, MECP2, PHF2, ETV7, and CDK9 were the top five transcription factors with the potential to alter the expression level of DDX56 (Figure 6E). Furthermore, we estimated the co-expression relation of these top five transcription factors and DDX56 (Spearman correlation, Supplementary Figure S5). We compared the expression levels of the top five transcription factors in tumor tissues and normal tissues (Wilcoxon test, Supplementary Figure S6). For example, we observed that the expression of CDK9 was highly related to the expression of DDX56 in THCA (Spearman correlation coefficient 0.384, p < 0.05, Supplementary Figure S5 and Supplementary Table S4). Simultaneously, we observed that the expression levels of CDK9 and DDX56 were both high in THCA tumor tissue (Wilcoxon test, p < 0.05, Supplementary Figure S6 and Figure 1A). We thus consider that CDK9 might be responsible for the high expression of DDX56 in THCA.