3.1. BEND family genes expression level in pan-cancers
We first analyzed RNA sequencing data present in TCGA using R software to evaluate the transcriptional levels of BEND in different human tumors(Fig.1A). According to the results, BEND3, BEND5 and, BEND7 were highly expressed in pan-cancers, BEND3P1, BEND3P3, and BEND6 were moderately expressed, while BEND2, BEND3P2, and BEND4 were low expressed in tumor tissues. Further analysis showed that BEND5 expression is the highest in CHOL. The expression level of BEND7 was highest in CHOL. The expression of BEND3 was highest in CHOL. The expression of BEND4 was highest in prostatic PRAD. The expression of BEND6 was highest in LUSC. The expression of BEND3P3 was highest in CHOL. The expression of BEND3P1 was highest in CHOL(Fig.1B).BEND5 and BEND6 were the two genes with the most significant positive correlation(Correlation coefficient=0.32, Fig.1C), while BEND3 and BEND5 were the two genes with the most significant negative correlation(Correlation coefficient=-0.16, Fig.1C).
We also used R software to analyze the RNA sequencing data in the TCGA database to evaluate the differential expression of BEND family genes in pan-cancers. The 11,075 mRNA expression profiles of 33 cancers from TCGA showed that BEND2 (Fig.2A) is highly expressed in many types of tumors, including breast invasive carcinoma(BRCA), Kidney Chromophobe(KICH), kidney renal clear cell carcinoma(KIRC), kidney renal papillary cell carcinoma(KIRP), Liver hepatocellular carcinoma(LIHC), lung adenocarcinoma(LUAD), lung adenocarcinoma (LUSC), Prostate adenocarcinoma(PRAD) and Thyroid carcinoma(THCA).BEND3(Fig.2B) in bladder urothelial carcinoma(BLCA),BRCA,Cholangiocarcinoma(CHOL), colon adenocarcinoma(COAD),glioblastoma multiforme(GBM),head and neck squamous cell carcinoma(HNSC),LIHC,LUAD,LUSC,PRAD,Rectum adenocarcinoma(READ),Stomach adenocarcinoma(STAD),THCA and uterine corpus endometrial carcinoma (UCEC) had higher expressions;meanwhile,it had lower expression in KICH, KIRC, and KIRP. BEND3P1(Fig.2C) was significantly expressed in BLCA, CHOL, COAD, HNSC, LIHC, READ, STAD and was low in KICH, KIRC, KIRP, THCA, and UCEC. BEND3P2(Fig.2D) was highly expressed in HNSC, STAD, and THCA; while it was low in CHOL, KICH, KIRP, LIHC, LUAD, and LUSC. BEND3P3(Fig.2E) was significantly expressed in CHOL, but in BLCA, BRCA, COAD, GBM, KICH, KIRC, KIRP, LUAD, LUSC, PRAD, READ, THCA, and UCEC had lower expression levels. BEND4(Fig.2F) was significantly expressed in PRAD, but in BLCA, BRCA, COAD, ESCA, GBM, KIRP, LIHC, LUAD, LUSC, READ, THCA, and UCEC was low in expression. BEND5 (Fig.2G) was significantly expressed in CHOL and in BRCA, COAD, ESCA, HNSC, KICH, KIRP, LUAD, LUSC, PRAD, READ, STAD, THCA, and UCEC was low in expression. BEDN6(Fig.2H) was significantly expressed in BRCA, CHOL, ESCA, HNSC, KIRP, LIHC, LUAD, and LUSC, and was low in COAD, GBM, PRAD, READ, UCEC. BEDN7(Fig.2I) was significantly expressed in CHOL and LIHC, which was low in BRCA, HNSC, KICH, KIRC, KIRP, LUAD, LUSC, PRAD, THCA, and UCEC.
Then we used the HPA database to verify the differential expression of the BEND family genes in gastric cancer and normal tissues (Fig.3). According to the results, BEND3 was highly expressed in gastric cancer tissues, while BEND5 was down-regulated in gastric cancer tissues. Due to the lack of BEND3P1, BEND3P2, and BEND3P3 data in the HPA database, further immunohistochemical verification of their expression is needed.
3.2. Prognostic value of BEND family genes in pan-cancers
Next, we analyzed the prognostic value of BEND family genes for pan-cancers through different databases. The Kaplan-Meier survival curve suggested(Fig.4) that the expression of BEND family genes is related to the prognosis of several TCGA cancers.Compared with the low-gene expression group,the survival rate of patients in the high expression group was higher.The expression of BEND2 was protective in PRAD (OS:N=499,P=0.042).The expression of BEND3 in ACC(OS:N=79, P=0.007),KIRC(OS:N=531,P=0.033),LIHC(OS:N=368,P=0.021),MESO(OS:N=84,P=0.015),SARC(OS:N=260,P<0.001),SKCM(OS:N=457,P=0.040) was an unfavorable effect,however,which in LAML(OS:N=132,P=0.001),READ(OS:N=158,P=0.010) and UVM(OS:N=80,P=0.027) was protective.The expression of BEND3P1 played an unfavorable role in COAD(OS:N=448,P=0.047),LAML(OS:N=132,P=0.046) and UCEC (OS:N=544,P=0.035),while played a protective role in BLCA(OS:N=406,P=0.047).The expression of BEND3P2 played an protective effect in DLBC(OS:N=47,P=0.012)、HNSC(OS:N=501P=0.019)、LAML(OS:N=132,P=0.017)、PAAD(OS:N=177,P=0.045).BEND3P3 played an unfavorable role in MESO(OS:N=84,P=0.048),and which in KIRC(OS:N=531,P=0.003)、LGG(OS:N=524,P<0.001) and LUAD(OS:N=513,P=0.048) was protective.The expression of BEND4 had an unfavorable effect in DLBC(OS:N=47,P=0.008) and UVM (OS:N=80,P=0.037),and it played a protective role in SKCM(OS:N=457, P<0.001). The expression of BEND5 played an unfavorable role in ACC(OS:N=79,P=0.008)、READ(OS:N=158,P=0.012),and in BLCA(OS:N=406,P<0.001)、BRCA(OS:N=1082,P=0.005)、CESC(OS:N=293,P=0.048)、KIRP(OS:N=286,P=0.006)、LUAD(OS:N=513,P=0.002)、MESO(OS:N=84,P=0.015)、PAAD(OS:N=177,P=0.014)、PRAD(OS:N=499,P=0.007)、TGCT(OS:N=139,P=0.048)、UVM(OS:N=80,P=0.034) it played an protective role.The expression of BEND6 in BLCA(OS:N=406,P=0.013)、KIRC(OS:N=531,P<0.001)、MESO(OS:N=84,P<0.001)、UCEC(OS:N=544,P=0.006)、UVM(OS:N=80,P=0.009) was unfavorable.The expression of BEND7 in BLCA(OS:N=406,P=0.038)、KIRC(OS:N=531,P=0.036)、LGG(OS:N=524,P=0.002) waws protective,however which was adverse in HNSC(OS:N=501,P=0.016)、UCEC(OS:N=544,P=0.018).
Next, we further conducted a Cox analysis on the prognostic value of BEND family genes expression for pan-cancers. According to the results(Fig.5), the expression of BEND2 was an unfavorable factor for the prognosis of ACC, LICH, OV(HR>1, P<0.05). The expression of BNED3 was a protective factor for the prognosis of LAML, READ, UVM(HR<1, P<0.05), and its expression was a disadvantage factor of ACC, KIRP, LIHC, MESO, SARC, SKCM(HR>1, P<0.05). The expression of BEND3P1 was a protective factor for the prognosis of GBM(HR<1, P<0.05) and adversely affected the prognosis of COAD, ESCA, KIRP, UCEC(HR>1, P<0.05). The expression of BEND3P2 had no obvious correlation with the prognosis of pan-cancers. The expression of BEND3P3 was a protective factor for the prognosis of KIRC and LGG(HR<1, P<0.05), but it was an unfavorable factor for the prognosis of MESO and PCPG(HR>1, P<0.05). The expression of BEND4 was a protective factor for the prognosis of LGG(HR<1, P<0.05), and its expression was an unfavorable factor for the prognosis of DLBC, KIRP, PCPG, SARC, THCA(HR>1, P<0.05). The expression of BEND5 was a protective factor for the prognosis of BLCA, HNSC, KIRC, KIRP, LUAD, MESO, PAAD, PRAD, SKCM, UVM(HR<1, P<0.05), and its expression was an unfavorable factor for the prognosis of STAD(HR>1, P<0.05). The expression of BEND6 was an unfavorable factor for the prognosis of BLCA, KIRC, KIRP, LUAD, MESO, READ, SARC, STAD, UCEC(HR>1, P<0.05). The expression of BEND7 was a protective factor for the prognosis of GBM, KIRC, LGG, READ(HR<1, P<0.05), and its expression was an unfavorable factor for the prognosis of UCEC.
3.3. Relationship between the BEND family genes and tumor microenvironment correlation、tumor stem cells
The BEND family genes expression was positively or negatively correlated with Stromal Score(Fig.6A), Immune Score(Fig.6B), and Estimate Score(Fig.6C)in pan-cancers patients. Similarly, BEND family genes expression was positively or negatively correlated with Tumor Purity, RNAss, and DNAss in pan-cancers.BNED3 showed a significantly negative correlation with Stromal Score, Immune Score, and Estimate Score in various cancers, and showed a significantly positive correlation with Tumor Purity, RNAss, and DNAss, suggesting that BEND3 is involved in the occurrence of various tumors.BEND4, BEND5, and BEND6 were positively correlated with Stromal Score, Immune Score, and Estimate Score in various cancers, and significantly negatively correlated with Tumor Purity, RNAss, and DNAss.
3.4. Drug sensitivity analysis of BEND family genes
To analyze the potential correlation between BEND family genes expression and drug sensitivity in different cancer cell lines from the CellMiner™ database, we performed a correlation analysis. The results showed (Fig.7) that the expression of BEND4 is positively correlated with the sensitivity of Fludarabine, XK-469, Asparaginase, Benddamustine, Nelarabine, Chlorambucil, Chelerythrine, and Cytarabine. The expression of BEND5 was positively correlated with the sensitivity of Nelarabine, XK469, Chelerythrine, and PX-316. The expression of BEND7 was negatively correlated with the sensitivity of Tamoxifen, Pipamperone, Ixabepilone, VinorelbinePipamperone, Ixabepilone, and Vinorelbine.
3.5. Association between BEND family genes expression and Immune subtype in Stomach adenocarcinoma(STAD)
To analyze the potential correlation between BEND family genes expression and immune subtypes in STAD, we performed correlation analysis.Immune subtypes include C1(Wound Healing),C2(IFN-G dominant),C3(inflammatory),C4(lymphocyte depleted),C5(immunologically quiet) and C6(TGF-b dominant).According to the results, the expressions of BEND3, BEND3P3, BEND4, BEND5, and BEND6 were all related to the immune subtypes of STAD(Fig.8).BEND3 was highly expressed in C1, C2, and C4, and moderately expressed in C3 and C6.BEND4 was moderately expressed in C1, C2, and C3, and was lowly expressed in C4 and C6.BEND3P3 was highly expressed in C3 and moderately expressed in other subtypes. The expression of BEND5 was highest in C3, moderate in C1, C2, and C6, and low in C4. The expression of BEND6 was highest in C3 and C6, and moderate in other subtypes.
3.6. Association between BEND family genes expression and Clinical relevance(Age、Gender、Grade、Stage and TNM Stage)in Stomach adenocarcinoma(STAD)
We analyzed the correlation between BEND family genes expression in STAD and clinical relevance, including Age, Gender, Grade, Stage, and TNM stage. According to the results(Fig.9), the expressions of BEND3, BEND5, and BEND6 were correlated with Age. The expression of BEND3 was higher in gastric cancer patients that are > 65 years old, while BEND5 and BEND6 expression were higher in ≤65-years olds. The expressions of BEND4, BEND5, and BEND6 were significantly correlated with Grade. With the gradual increase of Grade, the expression of BEND4 and BEND5 gradually increased, and the expression was the highest in G3. The expression of BEND6 was high in G3 and low in G1 and G2. The expression of BEND4, BEND5, and BEND6 was significantly correlated with Stage.BEND4 was highly expressed in StageⅡ and StageⅢ, BEND5 was highly expressed in StageⅡ, StageⅢ, and StageⅣ, and BEND6 was highly expressed in all stages.BEND3P2, BEND4, BEND5, and BEND6 were correlated with the T stage.BEND3P2 and BEND4 showed low expression in all T stages, while BEND5 and BEND6 showed high expression in T2, T3, and T4 stages. There was no correlation between the BEND family genes expression and gender, N and M stage.
3.7. Correlation between BEND family genes expression and TME, Stemness index in STAD
We further explored the relationship between BEND family genes expression and TME and Stemness scores of gastric cancer and conducted correlation analysis. The results showed that the expression of BEND3 is positively correlated with RNAss(Fig.10), while the expressions of BEND3P3, BEND4, BEND5, and BEND6 are negatively correlated with RNAss. The expression of BEND3 was positively correlated with DNAss, while the expressions of BEND3P3, BEND4, BEND5, and BEND6 were negatively correlated with DNAss. The expressions of BEND3 and BEND7 were negatively correlated with Stromal Score, while the expressions of BEND3P3, BEND4, BEND5, and BEND6 were positively correlated with Stromal Score. The expressions of BEND3 and BEND7 were negatively correlated with Immune Score, while the expressions of BEND4, BEND5, and BEND6 were positively correlated with Immune Score. The expressions of BEND3 and BEND7 were negatively correlated with Estimate Score, while BEND3P3, BEND4, BEND5, and BEND6 were positively correlated with Estimate Score.
3.8. An BEND risk model predicted OS and DFS in patients with STAD
We compared the relationship between the BEND family genes expression level and the risk of tumor occurrence in STAD patients, and the expressions of BEND3P1 and BEND6 were significantly correlated with the risk of tumor occurrence (Fig.11A), so a prognostic model composed of BEND3P1 and BEND6 was constructed. The risk score of each patient was calculated according to the risk score calculation formula: risk score =(-0.272087×BEND3P1 expression)+ (0.6450065×BEND6 expression). Patients were assigned to the high-risk groups (n=175) and the low-risk group(n=175) based on the median risk score. We constructed a heat map to show the expression of two genes in the high-risk and low-risk groups (Fig.10A), where the expression of BEND6 in the high-risk groups was higher than that in the low-risk groups, while the expression of BEND3P1 was lower than that in the low-risk groups.Fig.11B showed the distribution of risk scores of STAD patients. After patients were divided into two groups, the risk scores increased from left to right.Fig.11C showed the distribution of survival status and survival time of patients with different risk scores. K-M curve was used to compare the difference in overall survival rate between high-risk groups and low-risk groups(Fig.11D, P=1.249E-02). The results showed that the OS of gastric cancer patients with high-risk scores is significantly lower than that of patients with low-risk scores. According to the univariate analysis between clinical characteristics, risk score, and OS of gastric cancer patients(Fig.11E, F), we drew a time-dependent ROC curve(Fig.11G) to evaluate the predictive value of various pathological features on the prognosis of gastric cancer.The risk score(AUC=0.619),age(AUC=0.578),gender(AUC=0.540),Grade(AUC=0.562),Stage (AUC=0.602),T stage(AUC=0.563),M stage (AUC=0.528) and N stage (AUC=0.577) both had high sensitivity and specificity. At the same time, we constructed a nomogram to predict OS and DFS using the BEND risk score and clinical case characteristics(Fig.11H). We calculated the score of gastric cancer patients according to the nomogram, then added them to obtain the total score to predict 1-year and 3-year survival rates of gastric cancer patients, which was beneficial to guide clinical treatment. The closer the calibration curve is to the diagonal, the more accurate the prediction result is. The calibration curve of the nomogram showed good accuracy in predicting 1-year and 3-year survival rates(Fig.11I, J). The ROC curves of 1 year(AUC=0.619) and 3 years(AUC=0.555) also showed that the predictive ability of the nomogram is very accurate(Fig.11K).
3.9.The relationship between STAD overall survival rate and clinicopathological characteristics by Univariate and multivariate COX analyses.
Univariate and multivariate Cox proportional hazards regression analysis were used to evaluate whether risk model could be used as an independent predictor of adverse survival in STAD patients.According to the results(Table.1),Age(HR=1.024;95%CI,1.005-1.043;p=0.013),Gender(HR=1.562;95%CI,1.020-2392;p=0.040),Stage(HR=1.532;95%CI,1.213-1.935;p<0.001),Tumor size(HR=1.303;95%CI,1.025-1.656;p=0.030),Lymph node(HR=1.254;95%CI,1.056-1.490;p=0.010),Distant metastasis(HR=2.133;95%CI,1.109-4.103;p=0.023) and BEND risk score(HR=2.342;95%CI,1.195-4.589;p=0.013) were risk factors for prognosis. Multivariate analysis showed Age(HR=1.041;95%CI,1.019-1.062;p<0.001),Gender(HR=1.647;95%CI,1.066-2.544;p=0.025) and BEND risk score(HR=3.276;95%CI,1.607-6.680;p=0.001) are independent risk factors for prognosis.
3.10. Relationship between BEND scores and gene expression profiles
We used cBioportal to analyze the mutation and copy number variation of the BEND family genes.The mutation frequencies of BEND2,BEND3,BEND4,BEND5,BEND6 and BEND7 were 3%,7%,2.9%,1.6%,2.9%,1.8%,respectively(Fig. 12A).There was a missense mutation during BEND2, BEND3, BEND4, BEND5, BEND6, and BEND7. Then we performed survival analysis on the two genes that constructed the risk model and found that the high expression of BEND6 suggests a poor prognosis(Fig.12B, P=0.027), and the AUC curve showed that it has a high diagnostic value(Fig.12C, AUC at 1 years=0.579; AUC at 2 years=0.607; AUC at 1 years=0.653), while the predictive effect of BEND3P1 was poor(Fig.12E, F).In addition, we found multiple cancer-related KEGG pathways through Gene Set Enrichment Analysis(GSEA). The BEND6-related signaling pathways were as follows(Fig.12D): BASAL-CELL-CARCINOMA, DNA-REPLICATION, ECM-RECEPTOR-INTERACTION, MELANOMA, and TGF-BETA-SIGNALING-PATHWAY. The BEND3P1-related signaling pathways were as follows(Fig.12G): BASAL-CELL-CARCINOMA, COLORECTAL-CANCER, DNA-REPLICATION, MELANOMA, NON-SMALL-CELL-LUNG-CANCER, NOTCH-SIGNALING-PATHWAY, and T-CELL-RECEPTOR-SIGNALING-PATHWAY.
3.11. Relationship between BEND family gene expression and Immune infiltration
Tumor-infiltrating lymphocytes have been proven to be independent prognostic factors in cancer patients, so we investigated the relationship between BEND family genes expression and immune infiltration. The bar chart showed the number of immune cells in each sample(Fig.13A). For STAD patients, the expression levels of BEND6 and BEND3P1 were mostly positively correlated with the infiltration of different immune cells. For tissues with high expression of BEND6(Fig.13B), the infiltration level of B cell naïve(P=0.026), T cells follicular helper(P=0.021), Monocytes(P<0.001), Macrophages M2(P=0.029), Mast cells resting(P<0.001) was higher than that of tissues with low expression of BEND6.Plasma cells(P=0.034), T cells CD4 memory activated(P=0.008), and Macrophages MO(P=0.042) showed higher infiltration in tissues with low BEND expression. Correlation test analysis(Fig.13C) showed that the BEND6 expression is significantly positively correlated with the infiltration levels of B cell(R=0.19,P=0.00028),Dendritic cells resting(R=0.12,P=0.017), Macrophages M2(R=0.11,P=0.037),Mast cells resting(R=0.33,P=8.9e-11) and Monocytes (R=0.21,P=4.1e-05),which was negatively correlated with the infiltration of Macrophages MO(R=-0.19,P=0.00027),T cells follicular helper(R=-0.11,P=0.037), Neutrophils(R=-0.13,P=0.012),Plasma cells(R=-0.11,P=0.032) and T cells CD4 memory activated(R=-0.18,P=7e-04).After taking the intersection results of the two tests, it could be concluded that the expression level of BEND6 is significantly related to B cells naïve, Plasma cells, T cells CD4 memory activated, T cells follicular helper, Macrophages M0, Macrophages M2, Dendritic cells resting and Neutrophils(Fig.13D).
For tissues with high expression of BEND3P1(Fig.14A), the infiltration level of T cells follicular helper(P<0.001) and Macrophages M1(P=0.035) was higher.However,the infiltration level of NK cells resting(P=0.019),Dendritic cells activated(P=0.030),Neutrophils(P=0.003) was higher in tissues with low BEND3P1 expression.Correlation test analysis(Fig.14C) showed BEND6 expression level is significantly positively correlated with infiltration level of Macrophages M1 (R=0.1,P=0.05),NK cells activated(R=0.11,P=0.033) and T cells follicular helper (R=0.19,P= 0.00028).However the BEND6 expression level was negatively correlated with infiltration level of Mast cells resting(R=-0.11, P=0.035),Neutrophils(R=-0.17,P=0.00088) and NK cells resting(R=-0.18,P=0.00075). After taking the intersection of the two test results, it could be concluded that the expression level of BEND3P1 is significantly correlated with T cells follicular helper, NK cells resting, NK cells activated, Macrophages M1, Mast cells activated, and Neutrophils(Fig.14B). All in all, BEND family genes play an important role in the immune infiltration of STAD.