Molecular Characteristics, Oncogenic Roles, and Relevant Immune and Pharmacogenomic Features of EVA1B in Colorectal Cancer
Colorectal cancer (CRC) involving the colon and rectum ranks the third leading cause of cancer mortality globally, with over 1.85 million cases and 850,000 deaths each year . Recently, CRC presents increasing incidence among younger cases . This disease is a heterogeneous disease with distinct pathogenesis mechanisms, involving somatic mutation, genetic fusion, genetic instability, as well as epigenetic alteration . Among newly diagnosed CRC, 20% of cases have occurring metastasis at diagnosis as well as 25% will metastasize following localized disease . Surgical resection is the major therapeutic regimen of CRC. Nevertheless, there are very few treatment regimens for metastatic patients. Although chemotherapy is usually recommended, merely few targeted therapies appropriate for cases who have specific mutation profiling are available, such as EGFR inhibitor and VEGF inhibitor . Hence, it is urgently required to develop novel molecule targets against CRC. Immunotherapy is a novel alternative against CRC treatment, which utilizes cases’ own immune system to combat tumor cells . A minority of CRC cases present microsatellite instability (MSI), a molecule predictor of defective DNA mismatch repair, and a predictive biomarker of immunotherapeutic response, but most of them are microsatellite-stable (MSS) . Enhanced tumor mutational burden (TMB) and neoantigen load in MSI cancers sustain the infiltrations of immune effector subpopulations, and anti-cancer immune response is strongly in relation to MSS counterparts . Nevertheless, a few MSS cancers present increased TMB as well as infiltrating immune cells to respond to immunotherapy. Hence, novel predictors of immunotherapeutic responses will be required. EVA1B is a protein coding gene, which is a critical paralog of EVA1A gene. It is an endoplasmic reticulum and lysosome-relevant protein involving autophagy and apoptosis . Previously, it triggers papillary thyroid carcinogenesis and epithelial-mesenchymal transition (EMT) through Hippo signaling . Flubendazole exerts an antitumor role through modulating EVA1A-mediated autophagy and apoptosis in breast carcinoma . MiR-125b relieves oxaliplatin resistance in liver carcinoma through negatively modulating EVA1A-mediated autophagy . EVA1B possesses high sequence similar to EVA1A gene, and EVA1B protein presents the similar domain to EVA1A protein . Nevertheless, the role of EVA1B in CRC remains indistinct. This study was conducted to investigate the correlation of EVA1B expression with prognosis, tumor immune and pharmacogenomic features in CRC. 

Pan-Cancer Analyses RNA sequencing data (FPKM), somatic mutational data [Mutation Annotation Format (MAF)], copy number variation (CNV) data, and clinical information of 33 cancer 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); lymphoid neoplasm diffuse large B-cell lymphoma (DLBC); 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); acute myeloid leukemia (LAML); 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)] were curated from The Cancer Genome Atlas (TCGA) through the GDC data portal ( https://portal.gdc.cancer.gov/ ). Tumor Immune Estimation Resource (TIMER; version 2.0) web tool ( http://timer.cistrome.org/ ) was adopted to analyze the expression of EVA1B across 33 cancer types . For TCGA dataset, FPKM value was transformed to TPM value. The expression profiling of four major DNA methyltransferases (DNMT1, DNMT2, DNMT3A, and DNMT3B), DNA mismatch repair genes (MLH1, MSH2, MSH6, PMS2, and EPCAM), and m 6 A regulators (CBLL1, VIRMA, METTL3/14, RBM15/15B, WTAP, ZC3H13, ALKBH5, FTO, ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1/2/3, LRPPRC, YTHDC1/2, and YTHDF1/2/3) was extracted from pan-cancer specimens. TMB and MSI of pan-cancer were also harvested from TCGA project. Spearman correlation analysis was conducted to investigate the correlation of EVA1B with above factors across pan-cancer. Estimation of Immunological Features The abundance of tumor-infiltrating immune cells was inferred with single-sample gene set enrichment analysis (ssGSEA) derived from gene set variation analysis (GSVA) package . In total, 122 immunomodulatory factors containing MHCs, receptors, chemokines, and immune stimulators were collected from Charoentong et al. . Additionally, immune checkpoint molecules were retrieved from the study of Auslander et al. . Overall infiltration of stromal and immune cells was estimated in CRC tissues with estimation of stromal and immune cells in malignant tumors using expression data (ESTIMATE) algorithm in accordance with mRNA expression data . Cancer immunity cycle was curated from previous research , and the activities of all steps were estimated with ssGSEA . Collection of CRC Datasets The GSE14333 , GSE39582 , and GSE87211  datasets were collected from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/ ) project. Batch effects were adjusted with ComBat algorithm derived from sva package . Somatic mutational data were visualized with maftools package . GISTIC2.0 was adopted to analyze amplification and deletion using CNV data . Raw “CEL” files of microarray profiling from Affymetrix platform were downloaded, followed by background correction and quantile normalization via robust multiarray averaging algorithm utilizing affy and simpleaffy packages . Meanwhile, normalized matrix files of microarray profiling from other platforms were directly curated. Curation of Gene Sets of Known Biological Processes Gene sets of known biological processes containing CD8+ T effector, DNA damage repair, pan-fibroblast TGF-β response signature (pan-F-TBRS), antigen processing machinery, immune checkpoints, EMT1-3, FGFR3-relevant gene signatures, angiogenesis, KEGG discovered histones, Fanconi anemia, cell cycle, DNA replication, nucleotide excision repair, homologous recombination, mismatch repair, WNT target, as well as cell cycle regulators were curated from previous research . The activities of biological processes were quantified with ssGSEA algorithm. Functional Enrichment Analysis GSEA was conducted to investigate the differences in signaling pathways activated in two subpopulations, with the gene set “c2.cp.kegg.v6.2.symbols.gmt” as the reference. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of EVA1B-relevant genes was conducted with clusterProfiler package . GO categories contained biological process (BP), cellular component (CC), and molecular function (MF). Hallmark gene sets were harvested from the Molecular Signatures Database , and their activities were quantified with ssGSEA algorithm. Screening EVA1B-Relevant Genes To identify EVA1B-relevant genes, CRC patients were classified into high and low EVA1B expression subpopulations with the mean value of EVA1B expression. Through empirical Bayesian method from limma package , differentially expressed genes between subpopulations were determined. The significance criteria of EVA1B-relevant genes were set as |fold-change| >1.5 and adjusted p -value <0.01. Exploitation of an EVA1B-Derived Genomic Model Through univariate Cox regression models, interaction of EVA1B-relevant genes with CRC prognosis was estimated in TCGA cohort. EVA1B-relevant genes with p -value <0.05 were determined as prognostic factors, which were input into random survival forest analysis . The number of Monte Carlo iterations was set as 100, and the number of steps forward was 5. Thereafter, their relative importance was ranked. In accordance with relative importance >0.28, characteristic EVA1B-relevant genes were determined, which were input into a multivariate Cox regression model. The formula of EVA1B-relevant genomic model was conducted as follows: risk score = Σ k − 1 n Exp i ∗ e HR i , in which n indicated the number of characteristic EVA1B-relevant genes, Exp i indicated the expression of characteristic genes, and e HR i indicated the regression coefficients of genes derived from the multivariate Cox regression analysis. In accordance with this formula, risk score of each CRC patient was calculated. Thereafter, CRC patients were clarified into high- and low-risk subpopulations following the mean value of risk score. The risk score distribution, survival state, and heatmap of the expression of characteristic EVA1B-relevant genes were drawn following patients’ risk score. Kaplan–Meier curves of overall survival (OS), disease-specific survival (DSS), and progression-free survival (PFS) were depicted between high- and low-risk subpopulations using survival and survminer packages, followed by log-rank test. Thereafter, receiver operating characteristic (ROC) curves of 1-, 3-, and 5-year survival were established utilizing survivalROC package, and area under the curve (AUC) was calculated to assess the discrimination . Additionally, prognostic value of the EVA1B-derived genomic model was externally verified in the GSE14333, GSE39582, and GSE87211 cohorts. Assessment of Drug Sensitivity The Genomics of Drug Sensitivity in Cancer (GDSC) project ( https://www.cancerrxgene.org/ ) offers drug sensitivity information of 138 anticancer agents across approximately 75,000 experiments in 700 cancer cell lines . The half-maximal inhibitory concentration (IC50) that represented drug response was estimated with pRRophetic package . Nomogram Establishment Uni- and multivariate Cox regression models were established to evaluate the associations of clinical indicators (age, gender, stage, T, N, and M) and EVA1B-relevant risk score with CRC patients’ OS. ROC curves were conducted to evaluate the prognostic value of this nomogram. Additionally, calibration curves were utilized to assess the consistency between the actual and nomogram-predicted survival probabilities. Clinical Specimens and Data Source Formalin-fixed, paraffin-embedded specimens, including primary carcinoma specimens ( n = 131), corresponding noncancerous normal tissues ( n = 19) and liver metastasis specimens ( n = 19) for immunohistochemistry (IHC) were obtained from 131 CRC patients who underwent surgery between 2014 and 2015. CRC primary carcinomas were assessed according to the 8th edition of the American Joint Committee on Cancer (AJCC) staging system. CRC patient data and tissue samples were obtained from the Liaoning Cancer Hospital. Patients included in the current study did not receive preoperative radiotherapy or chemotherapy prior to the study. The current study (20210804GP) was approved by the Ethics Committee of the Liaoning Cancer Hospital, Shenyang, China. Immunohistochemistry EVA1B expression levels in samples were determined by IHC. Tissue sections were incubated with EVA1B antibody (Manufacturer: Atlas Antibodies; Catalog number: HPA043537). The level of EVA1B expression was determined by counting the percentages of positively stained immunoreactive cells and evaluating cell staining intensity. EVA1B IHC staining was scored as “−” ~ “+” and “++” ~ “+++,” which represented negative and positive, respectively. All samples were reviewed by two independent, experienced pathologists in who were blinded to the identity of the samples. Statistical Analysis Statistical analysis was implemented with R software (v3.4.1; https://www.r-project.org/ ) and its appropriate packages. Comparison between groups was conducted utilizing Student’s t - or Wilcoxon rank-sum test. Spearman correlation test was adopted to determine the interactions between variables. Chi-square test was performed to analyze the correlation between EVA1B and clinicopathologic characteristics. Kaplan–Meier curves and the log-rank test were used to compare the OS between EVA1B expression level. The prognostic ability of the predictors for OS was evaluated by ROC curves and the AUC values. Univariate and multivariate Cox regression analyses were utilized to evaluate the independent prognostic value of EVA1B regarding OS. p -value <0.05 was set as the threshold. 

RNA sequencing data (FPKM), somatic mutational data [Mutation Annotation Format (MAF)], copy number variation (CNV) data, and clinical information of 33 cancer 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); lymphoid neoplasm diffuse large B-cell lymphoma (DLBC); 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); acute myeloid leukemia (LAML); 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)] were curated from The Cancer Genome Atlas (TCGA) through the GDC data portal ( https://portal.gdc.cancer.gov/ ). Tumor Immune Estimation Resource (TIMER; version 2.0) web tool ( http://timer.cistrome.org/ ) was adopted to analyze the expression of EVA1B across 33 cancer types . For TCGA dataset, FPKM value was transformed to TPM value. The expression profiling of four major DNA methyltransferases (DNMT1, DNMT2, DNMT3A, and DNMT3B), DNA mismatch repair genes (MLH1, MSH2, MSH6, PMS2, and EPCAM), and m 6 A regulators (CBLL1, VIRMA, METTL3/14, RBM15/15B, WTAP, ZC3H13, ALKBH5, FTO, ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1/2/3, LRPPRC, YTHDC1/2, and YTHDF1/2/3) was extracted from pan-cancer specimens. TMB and MSI of pan-cancer were also harvested from TCGA project. Spearman correlation analysis was conducted to investigate the correlation of EVA1B with above factors across pan-cancer. 

The abundance of tumor-infiltrating immune cells was inferred with single-sample gene set enrichment analysis (ssGSEA) derived from gene set variation analysis (GSVA) package . In total, 122 immunomodulatory factors containing MHCs, receptors, chemokines, and immune stimulators were collected from Charoentong et al. . Additionally, immune checkpoint molecules were retrieved from the study of Auslander et al. . Overall infiltration of stromal and immune cells was estimated in CRC tissues with estimation of stromal and immune cells in malignant tumors using expression data (ESTIMATE) algorithm in accordance with mRNA expression data . Cancer immunity cycle was curated from previous research , and the activities of all steps were estimated with ssGSEA . 

The GSE14333 , GSE39582 , and GSE87211  datasets were collected from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/ ) project. Batch effects were adjusted with ComBat algorithm derived from sva package . Somatic mutational data were visualized with maftools package . GISTIC2.0 was adopted to analyze amplification and deletion using CNV data . Raw “CEL” files of microarray profiling from Affymetrix platform were downloaded, followed by background correction and quantile normalization via robust multiarray averaging algorithm utilizing affy and simpleaffy packages . Meanwhile, normalized matrix files of microarray profiling from other platforms were directly curated. 

Gene sets of known biological processes containing CD8+ T effector, DNA damage repair, pan-fibroblast TGF-β response signature (pan-F-TBRS), antigen processing machinery, immune checkpoints, EMT1-3, FGFR3-relevant gene signatures, angiogenesis, KEGG discovered histones, Fanconi anemia, cell cycle, DNA replication, nucleotide excision repair, homologous recombination, mismatch repair, WNT target, as well as cell cycle regulators were curated from previous research . The activities of biological processes were quantified with ssGSEA algorithm. 

GSEA was conducted to investigate the differences in signaling pathways activated in two subpopulations, with the gene set “c2.cp.kegg.v6.2.symbols.gmt” as the reference. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of EVA1B-relevant genes was conducted with clusterProfiler package . GO categories contained biological process (BP), cellular component (CC), and molecular function (MF). Hallmark gene sets were harvested from the Molecular Signatures Database , and their activities were quantified with ssGSEA algorithm. 

To identify EVA1B-relevant genes, CRC patients were classified into high and low EVA1B expression subpopulations with the mean value of EVA1B expression. Through empirical Bayesian method from limma package , differentially expressed genes between subpopulations were determined. The significance criteria of EVA1B-relevant genes were set as |fold-change| >1.5 and adjusted p -value <0.01. 

Through univariate Cox regression models, interaction of EVA1B-relevant genes with CRC prognosis was estimated in TCGA cohort. EVA1B-relevant genes with p -value <0.05 were determined as prognostic factors, which were input into random survival forest analysis . The number of Monte Carlo iterations was set as 100, and the number of steps forward was 5. Thereafter, their relative importance was ranked. In accordance with relative importance >0.28, characteristic EVA1B-relevant genes were determined, which were input into a multivariate Cox regression model. The formula of EVA1B-relevant genomic model was conducted as follows: risk score = Σ k − 1 n Exp i ∗ e HR i , in which n indicated the number of characteristic EVA1B-relevant genes, Exp i indicated the expression of characteristic genes, and e HR i indicated the regression coefficients of genes derived from the multivariate Cox regression analysis. In accordance with this formula, risk score of each CRC patient was calculated. Thereafter, CRC patients were clarified into high- and low-risk subpopulations following the mean value of risk score. The risk score distribution, survival state, and heatmap of the expression of characteristic EVA1B-relevant genes were drawn following patients’ risk score. Kaplan–Meier curves of overall survival (OS), disease-specific survival (DSS), and progression-free survival (PFS) were depicted between high- and low-risk subpopulations using survival and survminer packages, followed by log-rank test. Thereafter, receiver operating characteristic (ROC) curves of 1-, 3-, and 5-year survival were established utilizing survivalROC package, and area under the curve (AUC) was calculated to assess the discrimination . Additionally, prognostic value of the EVA1B-derived genomic model was externally verified in the GSE14333, GSE39582, and GSE87211 cohorts. 

The Genomics of Drug Sensitivity in Cancer (GDSC) project ( https://www.cancerrxgene.org/ ) offers drug sensitivity information of 138 anticancer agents across approximately 75,000 experiments in 700 cancer cell lines . The half-maximal inhibitory concentration (IC50) that represented drug response was estimated with pRRophetic package . 

Uni- and multivariate Cox regression models were established to evaluate the associations of clinical indicators (age, gender, stage, T, N, and M) and EVA1B-relevant risk score with CRC patients’ OS. ROC curves were conducted to evaluate the prognostic value of this nomogram. Additionally, calibration curves were utilized to assess the consistency between the actual and nomogram-predicted survival probabilities. 

Formalin-fixed, paraffin-embedded specimens, including primary carcinoma specimens ( n = 131), corresponding noncancerous normal tissues ( n = 19) and liver metastasis specimens ( n = 19) for immunohistochemistry (IHC) were obtained from 131 CRC patients who underwent surgery between 2014 and 2015. CRC primary carcinomas were assessed according to the 8th edition of the American Joint Committee on Cancer (AJCC) staging system. CRC patient data and tissue samples were obtained from the Liaoning Cancer Hospital. Patients included in the current study did not receive preoperative radiotherapy or chemotherapy prior to the study. The current study (20210804GP) was approved by the Ethics Committee of the Liaoning Cancer Hospital, Shenyang, China. 

EVA1B expression levels in samples were determined by IHC. Tissue sections were incubated with EVA1B antibody (Manufacturer: Atlas Antibodies; Catalog number: HPA043537). The level of EVA1B expression was determined by counting the percentages of positively stained immunoreactive cells and evaluating cell staining intensity. EVA1B IHC staining was scored as “−” ~ “+” and “++” ~ “+++,” which represented negative and positive, respectively. All samples were reviewed by two independent, experienced pathologists in who were blinded to the identity of the samples. 

Statistical analysis was implemented with R software (v3.4.1; https://www.r-project.org/ ) and its appropriate packages. Comparison between groups was conducted utilizing Student’s t - or Wilcoxon rank-sum test. Spearman correlation test was adopted to determine the interactions between variables. Chi-square test was performed to analyze the correlation between EVA1B and clinicopathologic characteristics. Kaplan–Meier curves and the log-rank test were used to compare the OS between EVA1B expression level. The prognostic ability of the predictors for OS was evaluated by ROC curves and the AUC values. Univariate and multivariate Cox regression analyses were utilized to evaluate the independent prognostic value of EVA1B regarding OS. p -value <0.05 was set as the threshold. 

Analysis of Expression, Genetic and Epigenetic Alterations, and Immunological Characteristics of EVA1B Across Pan-Cancer depicts the workflow of this study. EVA1B, as known as FAM176B, displayed remarkably increased expression in most cancer types in comparison with corresponding normal tissues, containing BLCA, BRCA, CHOL, COAD, ESCA, GBM, HNSC, KIRC, LIHC, PRAD, READ, STAD, and THCA . Differently, downregulated EVA1B was investigated in CESC, KICH, KIRP, and UCEC. DNA methyltransferases exert critical roles in changing chromatin structure and gene expression. We noted that EVA1B presented negative interactions with four major DNA methyltransferases containing DNMT1, DNMT2, DNMT3A, and DNMT3B in most cancer types . Additionally, EVA1B was negatively correlated to DNA mismatch repair genes (MLH1, MSH2, MSH6, PMS2, and EPCAM) across cancer types . M 6 A methylation represents the most common RNA modification, which affects the complexity of cancer progression . Especially, we focused on the interactions of EVA1B with m 6 A regulators in CRC. With the mean value of EVA1B expression, CRC patients were clustered into high and low EVA1B expression subpopulations. As depicted in , downregulated EVA1B was correlated to most m 6 A regulators like CBLL1, VIRMA, METTL14, and RBM15. TMB becomes a reliable biomarker of sensitivity to immune checkpoint blockage . In , EVA1B presented remarkably positive interaction with TBM across THYM, and LGG. Oppositely, it was negatively correlated with TMB across BRCA, CESC, ESCA, HNSC, STAD, and UCEC. MSI is a hypermutation phenotype caused by frequent polymorphisms in short repetitive DNA sequence as well as single nucleotide substitution due to DNA MMR defects . We noted the positive interaction of EVA1B with MSI in THCA, BLCA, DLBC, HNSC, KIRC, and PRAD but negative interaction with MSI in COAD, READ, STAD, and UCEC . Thereafter, we analyzed the immunological characteristics of EVA1B across pan-cancer. Our data suggested that EVA1B was remarkably correlated to immune cell infiltration , MCH molecules , chemokines , immunostimulatory factors , receptors , and immune checkpoint molecules  across pan-cancer. Prognostic Significance of EVA1B and Its Association With Drug Sensitivity in CRC We also focused on the biological role and significance of EVA1B in CRC. We combined the COAD and READ datasets from TCGA project and removed batch effects . Survival analysis was indicative that patients with high EVA1B expression possessed more undesirable OS outcomes than those with low EVA1B expression . Targeted therapies play a key role in CRC management, and genetic alterations that attribute to CRC heterogeneity is correlated to the responses to targeted therapies . Reliable predictive biomarkers for targeted therapies remain greatly lacking. Thus, we evaluated whether EVA1B expression was correlated to drug sensitivity. As a result, EVA1B presented positive association with sensitivity to bleomycin, lenvatinib, zoledronate, floxuridine, mitoxantrone, simvastatin, topotecan, and dasatinib . Differently, EVA1B expression was negatively correlated to sensitivity to nilotinib, tamoxifen, panobinostat, docetaxel, vinorelbine, crizotinib, ethinyl estradiol, and paclitaxel in CRC. Therefore, EVA1B might be a predictive marker of above agents. Immunological and Biological Significance of EVA1B in CRC Our ssGSEA results demonstrated the remarkably increased infiltration levels of most immune cells in high EVA1B expression subpopulation . Additionally, high EVA1B expression was correlated to increased stromal and immune score as well as reduced tumor purity in CRC . There were prominently enhanced activities of most steps within cancer immunity cycle in high EVA1B expression subpopulation . We also investigated the increased activities of immune activation processes (like CD8+ T effector, and immune checkpoint) and stromal activation processes (like EMT, pan-F-TBRS, FGFR3-related genes, and angiogenesis) in high EVA1B expression subpopulation . Additionally, the interactions of EVA1B expression with cancer immunity cycle and known biological processes were evaluated across CRC specimens. In , EVA1B expression presented positive correlations to most steps within cancer immunity cycle as well as immune and stromal activation processes. Association of EVA1B With Mutational Landscape in CRC In , we evaluated the prevalence of somatic mutation in high and low EVA1B expression subpopulations. APC, TP53, and TTN ranked the first three mutational genes. Nevertheless, no prominent difference in somatic mutation was investigated in high and low EVA1B expression subpopulations. The GISTIC2.0 results demonstrated that amplification and deletion displayed higher frequencies in high EVA1B expression subpopulation  compared with low EVA1B expression subpopulation . Additionally, we calculated the G-score in accordance with the amplitude of the aberrations and the frequency of their incidence in CRC specimens. Following comparison, mutations occurred in more regions in high EVA1B expression subpopulation . Signaling Pathways Associated With EVA1B We further analyzed the signaling pathways involving EVA1B via GSEA. In , high EVA1B expression was positively correlated to cell adhesion molecules (CAMs), chemokine signaling pathways, cytokine–cytokine receptor interaction, Leishmania interaction, leukocyte trans-endothelial migration, and VEGF signaling pathways. Meanwhile, low EVA1B expression was negatively correlated to cell cycle, homologous recombination, mismatch repair, nucleotide excision repair, oocyte meiosis, and spliceosome . The above data showed that EVA1B was involved in several oncogenic pathways, indicating the oncogenic role of EVA1B. Identification of EVA1B-Derived Genes and Their Biological Significance Through the cutoffs of |fold-change| >1.5 and adjusted p -value <0.01, we determined 602 EVA1B-derived genes in CRC individuals . Their biological significance was further investigated. In , EVA1B-derived genes were mainly correlated to immune response in accordance with GO annotation results. Additionally, there were remarkable interactions of EVA1B-derived genes with tumorigenic pathways (such as PI3K-Akt signaling pathway, proteoglycans in cancer, and ECM-receptor interaction) and immune activation pathways (such as IL-17 signaling pathway, Th17-cell differentiation, cytokine–cytokine receptor interaction, Th1- and Th2-cell differentiation, antigen processing and presentation, intestinal immune network for IgA production, complement and coagulation cascades; ). Overall, EVA1B-derived genes might exert remarkable roles in CRC progression. The EVA1B-Derived Genomic Model for CRC Prognosis In total, 66 EVA1B-derived genes were prominently correlated to CRC patients’ prognosis utilizing univariate-Cox regression models . With random survival forest analysis, we ranked the relative importance of above prognostic EVA1B-derived genes . In accordance with relative importance >0.28, six characteristic EVA1B-derived genes were determined, containing IGFBP3, EGFL7, SULT1B1, MMP1, CYP2W1, as well as ISYNA1 . On the basis of them, a multivariate-Cox regression model was conducted and risk scoring system was developed following the formula: risk score = (−0.117800092) * ISYNA1 expression + 0.067856101 * CYP2W1 expression + (−0.145078003) * MMP1 expression + (−0.107319254) * SULT1B1 expression + 0.427836868 * EGFL7 expression + 0.188926758 * IGFBP3 expression. With the mean value, we classified CRC patients into high- and low-risk subpopulations . Survival analysis uncovered that low-risk CRC individuals displayed a remarkable survival advantage . Nevertheless, no distinct difference in survival status was noted between high- and low-risk subpopulations . Additionally, the heterogeneity in expression of IGFBP3, EGFL7, SULT1B1, MMP1, CYP2W1, as well as ISYNA1 was displayed in two subpopulations . ROC curves confirmed that the EVA1B-derived genomic model possessed the well potency in estimating 1-, 3-, and 5-year OS probabilities . Verification of Prognostic Significance of EVA1B-Derived Genomic Model in CRC Except for OS, our analysis demonstrated that high-risk subpopulation was indicative of more undesirable DSS and PFS outcomes in comparison with low-risk subpopulation in TCGA cohort . The prognostic significance of the EVA1B-derived genomic model was verified in external cohorts. Our data confirmed that high-risk score was predictive of unfavorable OS outcome in the GSE14333, GSE39582, and GSE87211 datasets . The EVA1B-Derived Genomic Model Predicts Drug Responses of CRC Patients Further analysis was conducted to uncover whether this EVA1B-derived genomic model was predictive of drug responses of CRC individuals. In , we noted that high-risk group had significantly lower estimated IC50 of sunitinib and docetaxel relative to low-risk group, indicating that high-risk subpopulations presented higher sensitivity to sunitinib and docetaxel. Meanwhile, the low-risk group showed significantly lower estimated IC50 of sorafenib and gemcitabine, indicating that low-risk subpopulations were more likely to respond to sorafenib and gemcitabine. We also evaluated the interactions of EVA1B-derived genes with drug responses. As a result, MMP1 was negatively correlated to IC50 of mithramycin, actinomycin D, depsipeptide, and homoharringtonine; EGFL7 presented positive correlation to IC50 of raltitrexed, gemcitabine, cytarabine, and TFdU; IGFBP3 was positively correlated to IC50 of lenvatinib and JNJ-42756493 but negatively correlated to IC50 of LEE-011, cobimetinib (isomer 1), and selumetinib; CYP2W1 displayed positive relationship to IC50 of tegafur and fluorouracil . These data indicated that EVA1B-derived genes might be correlated to sensitivity to above agents. Establishment of a Reliable Nomogram for Prediction of CRC Prognosis Univariate-Cox regression analysis uncovered the remarkable association of age, stage, T, N, M, and EVA1B-derived risk score with CRC prognosis . Among them, age, T, and EVA1B-derived risk score served as independent prognostic indicators of CRC . By integrating above three independent prognostic indicators, this study exploited a nomogram to estimate CRC patients’ survival outcomes . This EVA1B-derived risk score occupied the most contribution to the prediction of 1-, 3-, and 5-year OS duration. ROC curves confirmed the desirable efficacy in predicting patients’ survival outcomes . Additionally, we evaluated the prediction performance of the nomogram with calibration curves. Our data demonstrated that 1-, 3-, and 5-year predicted by this nomogram was close to the actual survival duration . Above data indicated the superior predictive capacity of this nomogram. Molecular Mechanisms Involving the EVA1B-Derived Genomic Signature We further investigated the molecular mechanisms involving the EVA1B-derived genomic signature. As depicted in , unfolded protein response, E2F targets, G2M checkpoint, mTORC1 signaling, MYC targets, DNA repair, and oxidative phosphorylation presented higher activities in low-risk subpopulation. Meanwhile, tumorigenic pathways were remarkably activated in high-risk subpopulation, such as P53 pathway, TGF-beta signaling, apoptosis, Hedgehog signaling, EMT, Notch signaling, and angiogenesis, indicative of high-risk subpopulation’s unfavorable survival outcomes. Moreover, immune activation pathways such as allograft rejection, inflammatory response, and complement were significantly activated in high-risk subpopulation. Prognostic Significance of EVA1B-Derived Genes and Their Interaction With m 6 A Methylation and Immune Cell Infiltration demonstrates that TRADD, MXD4, NR1H2, EVA1B, and C11orf68 acted as risk factors of CRC outcomes. Additionally, we investigated that above EVA1B-derived genes presented positive correlations to m 6 A regulator ALKBH5 and RBM15B but negative correlations to other m 6 A regulators in CRC . In , these EVA1B-derived genes were negatively correlated to activated CD4+ T cell, effector memory CD4+ T cell, memory B cell, and type 2 T helper cell but positively correlated to other tumor-infiltrating immune cells in CRC. shows the association among drug sensitivity, predictive immune responses, and risk status of CRC individuals. EVA1B Overexpression Is Correlated With CRC Tumorigenesis and Poor Outcomes in CRC Patients EVA1B IHC was performed on 19 pairs of primary CRC, corresponding noncancerous and matching liver metastasis specimens from 19 liver metastasis CRC patients. The findings showed that EVA1B was overexpressed in primary CRC tissues compared with the corresponding noncancerous normal tissues ( p < 0.01). In addition, expression level of EVA1B was similar in the liver metastases samples to the expression level in matching primary CRC tissues without statistical significance . In addition, the relationship between EVA1B expression and clinicopathologic characteristics of 131 patients with CRC was explored . Expression of EVA1B was highly correlated with pN stage ( p < 0.001), pM stage ( p = 0.005), AJCC stage ( p = 0.001), CEA ( p = 0.003), and CA199 ( p = 0.002). However, the finding showed no significant correlation between EVA1B expression level and age, gender, tumor location, differentiation, and pT stage of CRC patients . Kaplan–Meier curve was generated to explore the prognosis value of EVA1B in these 131 CRC patients. The findings showed that CRC patients with high expression levels of EVA1B presented with poor OS compared with those with low expression levels of EVA1B ( p < 0.001, ). AUC values for EVA1B expression (0.630, 0.679, and 0.720 for the 2-, 4-, and 6-year OS, respectively) showed that EVA1B accurately discriminated between the high- and low-expression groups . Univariate and multivariate Cox analyses were performed, indicating that EVA1B was an independent as a prognostic factor by adjusting for stage, gender, location, differentiation, stage, CEA, and CA199 . 

depicts the workflow of this study. EVA1B, as known as FAM176B, displayed remarkably increased expression in most cancer types in comparison with corresponding normal tissues, containing BLCA, BRCA, CHOL, COAD, ESCA, GBM, HNSC, KIRC, LIHC, PRAD, READ, STAD, and THCA . Differently, downregulated EVA1B was investigated in CESC, KICH, KIRP, and UCEC. DNA methyltransferases exert critical roles in changing chromatin structure and gene expression. We noted that EVA1B presented negative interactions with four major DNA methyltransferases containing DNMT1, DNMT2, DNMT3A, and DNMT3B in most cancer types . Additionally, EVA1B was negatively correlated to DNA mismatch repair genes (MLH1, MSH2, MSH6, PMS2, and EPCAM) across cancer types . M 6 A methylation represents the most common RNA modification, which affects the complexity of cancer progression . Especially, we focused on the interactions of EVA1B with m 6 A regulators in CRC. With the mean value of EVA1B expression, CRC patients were clustered into high and low EVA1B expression subpopulations. As depicted in , downregulated EVA1B was correlated to most m 6 A regulators like CBLL1, VIRMA, METTL14, and RBM15. TMB becomes a reliable biomarker of sensitivity to immune checkpoint blockage . In , EVA1B presented remarkably positive interaction with TBM across THYM, and LGG. Oppositely, it was negatively correlated with TMB across BRCA, CESC, ESCA, HNSC, STAD, and UCEC. MSI is a hypermutation phenotype caused by frequent polymorphisms in short repetitive DNA sequence as well as single nucleotide substitution due to DNA MMR defects . We noted the positive interaction of EVA1B with MSI in THCA, BLCA, DLBC, HNSC, KIRC, and PRAD but negative interaction with MSI in COAD, READ, STAD, and UCEC . Thereafter, we analyzed the immunological characteristics of EVA1B across pan-cancer. Our data suggested that EVA1B was remarkably correlated to immune cell infiltration , MCH molecules , chemokines , immunostimulatory factors , receptors , and immune checkpoint molecules  across pan-cancer. 

We also focused on the biological role and significance of EVA1B in CRC. We combined the COAD and READ datasets from TCGA project and removed batch effects . Survival analysis was indicative that patients with high EVA1B expression possessed more undesirable OS outcomes than those with low EVA1B expression . Targeted therapies play a key role in CRC management, and genetic alterations that attribute to CRC heterogeneity is correlated to the responses to targeted therapies . Reliable predictive biomarkers for targeted therapies remain greatly lacking. Thus, we evaluated whether EVA1B expression was correlated to drug sensitivity. As a result, EVA1B presented positive association with sensitivity to bleomycin, lenvatinib, zoledronate, floxuridine, mitoxantrone, simvastatin, topotecan, and dasatinib . Differently, EVA1B expression was negatively correlated to sensitivity to nilotinib, tamoxifen, panobinostat, docetaxel, vinorelbine, crizotinib, ethinyl estradiol, and paclitaxel in CRC. Therefore, EVA1B might be a predictive marker of above agents. 

Our ssGSEA results demonstrated the remarkably increased infiltration levels of most immune cells in high EVA1B expression subpopulation . Additionally, high EVA1B expression was correlated to increased stromal and immune score as well as reduced tumor purity in CRC . There were prominently enhanced activities of most steps within cancer immunity cycle in high EVA1B expression subpopulation . We also investigated the increased activities of immune activation processes (like CD8+ T effector, and immune checkpoint) and stromal activation processes (like EMT, pan-F-TBRS, FGFR3-related genes, and angiogenesis) in high EVA1B expression subpopulation . Additionally, the interactions of EVA1B expression with cancer immunity cycle and known biological processes were evaluated across CRC specimens. In , EVA1B expression presented positive correlations to most steps within cancer immunity cycle as well as immune and stromal activation processes. 

In , we evaluated the prevalence of somatic mutation in high and low EVA1B expression subpopulations. APC, TP53, and TTN ranked the first three mutational genes. Nevertheless, no prominent difference in somatic mutation was investigated in high and low EVA1B expression subpopulations. The GISTIC2.0 results demonstrated that amplification and deletion displayed higher frequencies in high EVA1B expression subpopulation  compared with low EVA1B expression subpopulation . Additionally, we calculated the G-score in accordance with the amplitude of the aberrations and the frequency of their incidence in CRC specimens. Following comparison, mutations occurred in more regions in high EVA1B expression subpopulation . 

We further analyzed the signaling pathways involving EVA1B via GSEA. In , high EVA1B expression was positively correlated to cell adhesion molecules (CAMs), chemokine signaling pathways, cytokine–cytokine receptor interaction, Leishmania interaction, leukocyte trans-endothelial migration, and VEGF signaling pathways. Meanwhile, low EVA1B expression was negatively correlated to cell cycle, homologous recombination, mismatch repair, nucleotide excision repair, oocyte meiosis, and spliceosome . The above data showed that EVA1B was involved in several oncogenic pathways, indicating the oncogenic role of EVA1B. 

Through the cutoffs of |fold-change| >1.5 and adjusted p -value <0.01, we determined 602 EVA1B-derived genes in CRC individuals . Their biological significance was further investigated. In , EVA1B-derived genes were mainly correlated to immune response in accordance with GO annotation results. Additionally, there were remarkable interactions of EVA1B-derived genes with tumorigenic pathways (such as PI3K-Akt signaling pathway, proteoglycans in cancer, and ECM-receptor interaction) and immune activation pathways (such as IL-17 signaling pathway, Th17-cell differentiation, cytokine–cytokine receptor interaction, Th1- and Th2-cell differentiation, antigen processing and presentation, intestinal immune network for IgA production, complement and coagulation cascades; ). Overall, EVA1B-derived genes might exert remarkable roles in CRC progression. 

In total, 66 EVA1B-derived genes were prominently correlated to CRC patients’ prognosis utilizing univariate-Cox regression models . With random survival forest analysis, we ranked the relative importance of above prognostic EVA1B-derived genes . In accordance with relative importance >0.28, six characteristic EVA1B-derived genes were determined, containing IGFBP3, EGFL7, SULT1B1, MMP1, CYP2W1, as well as ISYNA1 . On the basis of them, a multivariate-Cox regression model was conducted and risk scoring system was developed following the formula: risk score = (−0.117800092) * ISYNA1 expression + 0.067856101 * CYP2W1 expression + (−0.145078003) * MMP1 expression + (−0.107319254) * SULT1B1 expression + 0.427836868 * EGFL7 expression + 0.188926758 * IGFBP3 expression. With the mean value, we classified CRC patients into high- and low-risk subpopulations . Survival analysis uncovered that low-risk CRC individuals displayed a remarkable survival advantage . Nevertheless, no distinct difference in survival status was noted between high- and low-risk subpopulations . Additionally, the heterogeneity in expression of IGFBP3, EGFL7, SULT1B1, MMP1, CYP2W1, as well as ISYNA1 was displayed in two subpopulations . ROC curves confirmed that the EVA1B-derived genomic model possessed the well potency in estimating 1-, 3-, and 5-year OS probabilities . 

Except for OS, our analysis demonstrated that high-risk subpopulation was indicative of more undesirable DSS and PFS outcomes in comparison with low-risk subpopulation in TCGA cohort . The prognostic significance of the EVA1B-derived genomic model was verified in external cohorts. Our data confirmed that high-risk score was predictive of unfavorable OS outcome in the GSE14333, GSE39582, and GSE87211 datasets . 

Further analysis was conducted to uncover whether this EVA1B-derived genomic model was predictive of drug responses of CRC individuals. In , we noted that high-risk group had significantly lower estimated IC50 of sunitinib and docetaxel relative to low-risk group, indicating that high-risk subpopulations presented higher sensitivity to sunitinib and docetaxel. Meanwhile, the low-risk group showed significantly lower estimated IC50 of sorafenib and gemcitabine, indicating that low-risk subpopulations were more likely to respond to sorafenib and gemcitabine. We also evaluated the interactions of EVA1B-derived genes with drug responses. As a result, MMP1 was negatively correlated to IC50 of mithramycin, actinomycin D, depsipeptide, and homoharringtonine; EGFL7 presented positive correlation to IC50 of raltitrexed, gemcitabine, cytarabine, and TFdU; IGFBP3 was positively correlated to IC50 of lenvatinib and JNJ-42756493 but negatively correlated to IC50 of LEE-011, cobimetinib (isomer 1), and selumetinib; CYP2W1 displayed positive relationship to IC50 of tegafur and fluorouracil . These data indicated that EVA1B-derived genes might be correlated to sensitivity to above agents. 

Univariate-Cox regression analysis uncovered the remarkable association of age, stage, T, N, M, and EVA1B-derived risk score with CRC prognosis . Among them, age, T, and EVA1B-derived risk score served as independent prognostic indicators of CRC . By integrating above three independent prognostic indicators, this study exploited a nomogram to estimate CRC patients’ survival outcomes . This EVA1B-derived risk score occupied the most contribution to the prediction of 1-, 3-, and 5-year OS duration. ROC curves confirmed the desirable efficacy in predicting patients’ survival outcomes . Additionally, we evaluated the prediction performance of the nomogram with calibration curves. Our data demonstrated that 1-, 3-, and 5-year predicted by this nomogram was close to the actual survival duration . Above data indicated the superior predictive capacity of this nomogram. 

We further investigated the molecular mechanisms involving the EVA1B-derived genomic signature. As depicted in , unfolded protein response, E2F targets, G2M checkpoint, mTORC1 signaling, MYC targets, DNA repair, and oxidative phosphorylation presented higher activities in low-risk subpopulation. Meanwhile, tumorigenic pathways were remarkably activated in high-risk subpopulation, such as P53 pathway, TGF-beta signaling, apoptosis, Hedgehog signaling, EMT, Notch signaling, and angiogenesis, indicative of high-risk subpopulation’s unfavorable survival outcomes. Moreover, immune activation pathways such as allograft rejection, inflammatory response, and complement were significantly activated in high-risk subpopulation. 

6 A Methylation and Immune Cell Infiltration demonstrates that TRADD, MXD4, NR1H2, EVA1B, and C11orf68 acted as risk factors of CRC outcomes. Additionally, we investigated that above EVA1B-derived genes presented positive correlations to m 6 A regulator ALKBH5 and RBM15B but negative correlations to other m 6 A regulators in CRC . In , these EVA1B-derived genes were negatively correlated to activated CD4+ T cell, effector memory CD4+ T cell, memory B cell, and type 2 T helper cell but positively correlated to other tumor-infiltrating immune cells in CRC. shows the association among drug sensitivity, predictive immune responses, and risk status of CRC individuals. 

EVA1B IHC was performed on 19 pairs of primary CRC, corresponding noncancerous and matching liver metastasis specimens from 19 liver metastasis CRC patients. The findings showed that EVA1B was overexpressed in primary CRC tissues compared with the corresponding noncancerous normal tissues ( p < 0.01). In addition, expression level of EVA1B was similar in the liver metastases samples to the expression level in matching primary CRC tissues without statistical significance . In addition, the relationship between EVA1B expression and clinicopathologic characteristics of 131 patients with CRC was explored . Expression of EVA1B was highly correlated with pN stage ( p < 0.001), pM stage ( p = 0.005), AJCC stage ( p = 0.001), CEA ( p = 0.003), and CA199 ( p = 0.002). However, the finding showed no significant correlation between EVA1B expression level and age, gender, tumor location, differentiation, and pT stage of CRC patients . Kaplan–Meier curve was generated to explore the prognosis value of EVA1B in these 131 CRC patients. The findings showed that CRC patients with high expression levels of EVA1B presented with poor OS compared with those with low expression levels of EVA1B ( p < 0.001, ). AUC values for EVA1B expression (0.630, 0.679, and 0.720 for the 2-, 4-, and 6-year OS, respectively) showed that EVA1B accurately discriminated between the high- and low-expression groups . Univariate and multivariate Cox analyses were performed, indicating that EVA1B was an independent as a prognostic factor by adjusting for stage, gender, location, differentiation, stage, CEA, and CA199 . 

This study presented an integrative analysis of molecular characteristics, oncogenic roles, and relevant immune and pharmacogenomic features of EVA1B across pan-cancer, especially CRC. Our data demonstrated the remarkable upregulation of EVA1B expression in most cancer types. Additionally, its upregulation was significantly correlated to DNA methyltransferases, DNA mismatch repair genes, m 6 A regulators, TMB, and MSI across pan-cancer. Tumorigenesis is a multistep process where normal cells acquire genetic and epigenetic alterations contributing to tumor initiation and progression . The pathological mechanisms of CRC initiation and progression include chromosomal instability, high MSI and methylation, contributing to oncogene, tumor suppressor gene, and mismatch repair-relevant gene mutations . Our data were indicative of the critical role of EVA1B in rectal and colon carcinogenesis. High EVA1B expression was indicative of undesirable CRC patients’ clinical outcomes. Nevertheless, we did not note the distinct difference in survival status between high- and low-risk subpopulations due to this retrospective cohort. The role of EVA1B in survival status of CRC should be investigated in a prospective cohort in our future studies. Previous research proposed the upregulation of EVA1B in glioma . At present, the therapeutic regimen is surgical resection plus chemotherapy in more advanced or inoperable patients . Nevertheless, immunotherapy to supplement curative and palliative therapy has been undergoing clinical trials . Our further analysis demonstrated that EVA1B upregulation was correlated to enhanced immune cell infiltration, increased stromal and immune activation, and elevated activities of cancer immunity cycle. Additionally, higher frequencies of amplification and deletion were noted in high EVA1B expression subpopulation. Our data were indicative that EVA1B might be utilized as a promising immunotherapeutic predictor. We further determined 602 EVA1B-derived genes in CRC that were mainly correlated to immune response, tumorigenic pathways (such as PI3K-Akt signaling pathway, proteoglycans in cancer, and ECM-receptor interaction), and immune activation pathways (such as IL-17 signaling pathway, Th17 cell differentiation, cytokine-cytokine receptor interaction, Th1 and Th2 cell differentiation, antigen processing and presentation, intestinal immune network for IgA production, complement and coagulation cascades), indicative of their remarkable roles in CRC progression. Utilizing random survival forest analysis, we conducted an EVA1B-derived genomic model, composed of IGFBP3, EGFL7, SULT1B1, MMP1, CYP2W1, and ISYNA1. Following verification in diverse external datasets, the genomic model reliably and independently predicted patients’ prognosis and relapse. Additionally, it possessed the potential in estimating drug responses, like sunitinib, sorafenib, gemcitabine, and docetaxel. There was remarkable correlation of characteristic EVA1B-derived genes with responses to small molecular compounds. TRADD, MXD4, NR1H2, EVA1B, and C11orf68 served as remarkable risk factors of CRC outcomes. Previously, modulation of TRADD restores cellular homeostasis as well as alleviates apoptosis . It triggers tumor inhibition through regulation of ULF-dependent p19Arf ubiquitylation . CDK4/6 inhibitor upregulates MXD4 expression that negatively modulates MYC in CD8+ T cells . NR1H2 modulates cholesterol homeostasis within human cells, controlling fitness and function of activated T cells . Nevertheless, more experimental evidence should be conducted to uncover the functions of above EVA1B-relevant genes in CRC. A nomogram acts as a powerful tool to quantify an individual’s risk in a clinical setting through integrating diverse risk factors . Herein, this study utilized the nomogram in prediction of 1-, 3-, and 5-year OS probabilities via incorporating age, T stage, as well as EVA1B-derived genomic model for CRC individuals. Each indicator was assigned a score on the basis of its contribution to survival risk. Thereafter, ROC curves demonstrated that the nomogram presented the favorable efficiency in predicting an individual patient’ OS outcome. Additionally, calibrated curves confirmed the actual survival duration agreed with the nomogram-estimated survival duration. A few limitations of our research will be pointed out. The functional role of EVA1B in tumor immunity requires in-depth experimental verification. Moreover, prognostic value of EVA1B needs to be verified in larger CRC cohorts. 

Overall, our study conducted an integrative analysis to uncover molecular characteristics, oncogenic roles, and relevant immune and pharmacogenomic features of EVA1B in CRC. Our findings indicated that EVA1B acted as a convincing prognostic marker as well as a predictor of therapeutic responses. 

Publicly available datasets were analyzed in this study. These data can be found here: https://portal.gdc.cancer.gov/ and GEO, https://www.ncbi.nlm.nih.gov/geo/ , under accession numbers GSE14333, GSE39582, and GSE87211. 

The studies involving human participants were reviewed and approved by the human ethics committees of the Liaoning Cancer Hospital. The patients/participants provided their written informed consent to participate in this study. 

BM and YML conceived and designed this study. YML and KW performed the bioinformatics analyses and visualization. YL collected data and performed the statistical analysis. BM, YML and QM wrote the original draft. BM and YML revised the manuscripts. All authors revised and approved the final manuscript. 

This study was supported by the National Natural Science Foundation of China (81902383); the Revitalizing Liaoning Talents Program (XLYC1907004); the Young and Middle-aged Scientific and Technological Innovation Talent Support Plan of Shenyang City (RC200223); Cultivation Program of National Science Foundation of Liaoning Cancer Hospital (2021-ZLLH-03); Cultivation Program of National Science Foundation of Liaoning Cancer Hospital (2020-ZLLH-44); Natural Guiding Plan Foundation of Liaoning Province (2019-ZD-0584); and Beijing Xisike Clinical Oncology Research Foundation (Y-QL202101-0039). 

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. 

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher. 

