High expression of COMMD7 is an adverse prognostic factor in acute myeloid leukemia

RNA-sequencing data and bioinformatics analysis

The pan-cancer RNA-seq data of TCGA and GTEx with toil processed uniformly were downloaded from UCSC XENA (https://xenabrowser.net/datapages/) [10–13]. Level 3 HTSeq-FPKM and HTSeq-Count data of the AML samples were obtained from the TCGA website (https://portal.gdc.cancer.gov/repository) for further analysis. This study was in full compliance with the published guidelines of TCGA and GTEx.

Differentially expressed gene (DEG) analysis

The DESeq2 R package was adopted to compare expression data of low- and high-expression of COMMD7 (cut-off value of 50%) in AML samples (HTseq-Count) to identify DEGs [14]. The top 10 DEGs were performed by heat map.

Functional enrichment analysis

DEGs with the threshold for |logFC| >1.5 and padj <0.05 were applied for functional enrichment analysis. Gene Ontology (GO) functional analysis comprising cellular component (CC), molecular function (MF), and biological process (BP), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, were implemented using the ClusteProfiler package in R [15].

Gene set enrichment analysis (GSEA)

R package ClusteProfiler (3.14.3) was used for GSEA to elucidate the functional and pathway differences between the high- and low-expression groups of COMMD7 [15]. The gene set was permutated 1,000 times for each analysis. Adjusted P-value < 0.05 and FDR q-value < 0.25 were considered to be statistically significant.

Immune infiltration analysis by single-sample Gene Set Enrichment Analysis (ssGSEA)

Immune infiltration analysis of COMMD7 was conducted by ssGSEA using GSVA package in R (3.6.3). A total of 24 types of infiltrating immune cells were obtained as previously described [16]. Spearman correction was used to analyze the correlation between COMMD7 and the enrichment scores of 24 types of immune cells. Wilcoxon rank-sum test was used to analyze the enrichment scores of high- and low-COMMD7 expression groups.

PPI network

The PPI network of DEGs was predicted using the Search Tool for the Retrieval of Interacting Genes (STRING) database [17]. The interaction score threshold of 0.4 was set as the cut-off criterion. The PPI network was mapped using Cytoscape (version 3.7.1) [18], and the most significant modules in the PPI network were identified using MCODE (version 1.6.1) [19]. Selection criteria were as follows: MCODE scores >5, degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2. Metascape (https://metascape.org/gp/index.htm) was used to conduct the pathway and process enrichment analysis.

Prognostic model generation and prediction

In order to individualize the prediction of overall survival (OS) and event-free survival (EFS) in AML patients, a nomogram was generated using the RMS R package (version 5.1-3), which included prominent clinical characteristics and calibration plots. The calibration curves were evaluated graphically by mapping the nomogram-predicted probabilities against the observed rates, and the 45°line represented the best predictive values. Concordance index (C-index) was used to determine the discrimination of the nomogram, and the bootstrap approach was used to calculate 1000 resamples. In addition, C-index and receiver operating characteristic (ROC) were used to analyze and compare the predictive accuracy of the nomogram and separate prognostic factors. All statistical tests were double-tailed with 0.05 as the statistical significance level.

Statistical analysis

All statistical analyses and graphs were analyzed and displayed by R (3.6.2) [20]. The expression of COMMD7 in unpaired samples was analyzed by Wilcoxon rank-sum test, with Wilcoxon signed-rank test used in paired samples. Kruskal-Wallis test, Wilcoxon signed-rank test, and logistic regression analysis were used to evaluate the relationship between clinical/cytogenetic characteristics and COMMD7 expression. Cox regression analysis and Kaplan-Meier method were used to evaluate the prognostic factors. Multivariate Cox analysis was adopted to compare the impact of COMMD7 expression on survival along with other clinical features. The median COMMD7 expression was regarded as the cut-off value. In all tests, P value < 0.05 was considered statistically significant. Moreover, ROC analysis was performed on the pROC package to assess the effectiveness of the transcriptional expression of COMMD7 in distinguishing AML from healthy samples. The computed area under the curve (AUC) value ranging from 0.5 to 1.0 indicated 50-100% discrimination ability.

COMMD7 expression in pan-cancers and AML

RNA-seq data from UCSC XENA (https://xenabrowser.net/datapages/) was downloaded in TCGA and GTEx formats processed uniformly through the toil process. By comparing the expression of COMMD7 normal samples in TCGA and GTEX database and corresponding tumor samples in TCGA database, COMMD7 was found significantly high expressed in 28 types of cancer (Figure 1A), including acute myelogenous leukemia (LAML) (Figure 1B).

Identification of DEGs in AML samples with low- and high-expressed COMMD7

The high- and low-expression groups' gene expression profiles were analyzed for differences in the median mRNA expression. A total of 529 DEGs from gene expression RNA-seq-HTSeq-Counts, including 92 up-regulated and 437 down-regulated, were identified statistically significant between COMMD7 high- and low-expressed groups (|log fold change (logFC)| > 1.5, P < 0.05) (Figure 2A). The top five up-regulated DEGs and top five down-regulated DEGs between COMMD7 high- and low-expressed groups were illustrated by the heat map (Figure 2B).

Functional enrichment analysis of DEGs

To better understand the functional implication of 529 DEGs between high- and low- expression of COMMD7 in AML, GO and KEGG functional enrichment analysis was performed by clusterProfiler package (Supplementary Table 1, Figure 3). The association with the biological process (BP) included pattern specification process, regionalization, and mesenchyme development; cellular components (CC) included collagen-containing extracellular matrix, ion channel complex, and basement membrane; molecular function (MF) included receptor ligand activity, DNA-binding transcription activator activity/RNA polymerase II-specific, extracellular matrix structural constituent. KEGG included PI3K-Akt signaling pathway, focal adhesion, and ECM-receptor interaction.

GSEA analysis was conducted to gain further insight into the biologic pathways involved in AML with different COMMD7 expression levels. GSEA was performed between low- and high-COMMD7 expression datasets to identify critical signaling pathways involved in AML. Significant differences (FDR < 0.05, ADJ P < 0.05) were observed in the enrichment of MSigDB Collection (C2.all.v7.0.symbols.gmt) of these pathways (Supplementary Table 2 and Figure 4). Gene mutations or fusions with a good prognosis of AML, such as PML-RARa fusion, NPM1 mutation, AML-ETO fusion, and CBFB-MYH11 fusion, were enriched in COMMD7 low- expression phenotype based on NES, with adjusted P value <0.05 and FDR value <0.05 (Figure 4A–4D). On the contrary, in the high expression of COMMD7 phenotypes, factors with poor prognosis in AML, such as FLT3-ITD fusion and MLL fusion, presented significantly enriched (Figure 4I–4J). So did the pathways involved in AML and other tumor development, such as MAPK, RAS, Hedgehog, and Wnt pathways (Figure 4E–4H). Other genetic variants, such as phosphorylated TP53 targets and MYC targets, were also significantly enriched in such phenotype (Figure 4K–4L).

Immune infiltration analysis in AML

Spearman correlation analysis showed that the expression level of COMMD7 in the AML microenvironment was correlated with the immune cell infiltration level quantified by SSGSEA. Specifically, COMMD7 was positively associated with NK CD56bright cells and active dendritic cells (aDCs) (Figure 5).

PPI enrichment analysis in AML

The network of COMMD7 and its potential co-expressed genes in COMMD7-related DEGs was constructed by STRING, with a threshold of 0.4 (Supplementary Table 3). A total of 529 DEGs were screened out ( |log fold change (logFC)| >1.5, P < 0.05). The PPI network with 238 nodes and 367 edges was displayed by Cytoscape-MCODE (Figure 6A). The most significant module with a MCODE score of 7.317 contained 42 nodes and 150 edges (Figure 6B). Meantime, Metascape-MCODE was used to identify densely connected PPI network components of COMMD7, shown in Supplementary Figure 1. The three best-scoring GO terms by p-value as the functional description of the corresponding components were shown in Supplementary Table 4.

Association between COMMD7 expression and clinical features and cytogenetic risks

The main clinical characteristics of AML in TCGA was shown in Table 1. A total of 151 cases (68 females and 83 males) were analyzed in this study, with an average age of 56.7 years. Among them, COMMD7 expression was low in 76 (50.3%) AML patients and high in the remaining 75 (49.3%) cases. The median COMMD7 expression (log2(TPM+1)), which is 5.783, was regarded as the cut-off value. Correlation analysis suggested that COMMD7 expression was significantly correlated with cytogenetic risk and white blood cell (WBC) count (×10⁹/L) (P < 0.001). In addition, COMMD7 expression was significantly associated with other factors including bone marrow (BM) (P = 0.007), peripheral blood (PB) blasts (%) (P = 0.005), FAB classification (P = 0.036), FLT3 mutation (P = 0.004), IDH1 R132 mutation (P = 0.046), and NPM1 mutation (P = 0.014).

Logistic analysis was applied to further verify the relationship between AML clinicopathological factors and the COMMD7 high-low dichotomy. As a result, high expression of COMMD7 showed a significant positive correlation with high WBC count (>20 × 10⁹/L) (odds ratio [OR], 3.16; P < 0.001) and high PB blasts (>70%) (OR, 2.89; P = 0.002), whereas negatively correlated with FLT3 mutation (OR, 0.32; P = 0.002) and NPM1 mutation (OR, 0.34; P = 0.01) (Table 2). What is more, the potential value of COMMD7 in differentiating AML patients from healthy individuals was examined by ROC curve analysis, with the AUC of 0.760, revealing that COMMD7 had potential as a biomarker (Figure 7A). Besides, the Wilcoxon Rank SUM test was used to compare the expression of COMMD7 in patients with different clinicopathological features. The result showed that COMMD7 was significantly high-expressed in the patients with BM blasts (>20%; P = 0.014), WBC counts (>20 × 10⁹/L; P = 0.002), FAB classification (non-M3 type; P = 0.019), cytogenetics risk (intermediate/poor; P < 0.001), NPM1 mutation (negative; P = 0.005), FLT3 mutation (negative; P = 0.004), and IDH1 R132 mutation (negative; P = 0.019) (Figure 7B–7H).

High COMMD7 impacted the prognosis of AML in patients with different clinicopathological status

The relationship between COMMD7 expression and prognosis was analyzed in AML patients by using Kaplan-Meier. As seen in Figure 8A, patients with high expression of COMMD7 had a strongly worse prognosis than those with low COMMD7 expression (hazard ratio [HR], 1.91(1.25-2.93); P = 0.003). Kaplan-Meier analysis presented that high expressed COMMD7 correlated with poor prognosis in the subgroups of BM blasts≥ 20% (P = 0.024), PB blasts ≤ 70% (P = 0.007), age >60 (P = 0.009), FLT3 mutation negative (P = 0.009), IDH1 R132 mutation positive (P = 0.001), R140 mutation positive (P = 0.002), R172 mutation positive (P = 0.001), NPM1 mutation positive (P < 0.001), RAS mutation positive (P = 0.002), RUX1 mutation negative (P = 0.004), and DNMT3A mutation negative (P = 0.019) (Figure 8B–8L).

Likewise, the forest plot illustrated the prognostic value of COMMD7 in various AML subtypes using univariate Cox regression, with a conclusion consistent with the above results (Figure 9).

Hereafter, univariate Cox proportional hazards regression was used to assess the factors influencing OS, disclosing that COMMD7 (high- vs. low-, P = 0.003) was a predictive factor for worse OS, so did cytogenetic risk (poor & intermediate vs. favorable, P < 0.001) and age (>60 vs. ≤60, P < 0.001) (Table 3). Cytogenetic risk, age and COMMD7 were then included in multivariate Cox regression, suggesting that age > 60 (P < 0.001) and high expression of COMMD7 (P = 0.01) were independent prognostic factors for worse OS (P < 0.05).

Prognostic model of COMMD7 in AML

To better predict AML patients' prognosis, a nomogram was constructed based on the Cox regression analysis results using the RMS R package (Figure 10A). Three independent prognostic factor variables, age, cytogenetic risk, and COMMD7 expression, were included in the model, selected into the prediction model at a statistical significance level of 0.2. Based on multivariate Cox analysis, a point scale was used to assign points to these variables. The straight line was drawn upward to determine the points of the variables, and the sum of the points assigned to each variable was rescaled to a range of 0–100. The points of each variable were accumulated and recorded as the total points. The probability of AML patient survival at 1-, 3-, and 5-year was determined by drawing a line from the total point axis straight down to the outcome axis. The 1-year survival probability was determined by drawing a vertical line downward on the total point axis along the 162-direction ending axis, suggesting the probability of 1-year survival < 20%, both of the probability of 3- and 5-year < 10%. The prediction results of the nomogram calibration curve of OS were consistent with all patients' observation results (Figure 10B).