Biological Function and Clinical Value of VPS13A in Pan-Cancer Based on Bioinformatics Analysis

Introduction

Given the complexity of tumor genesis and progression, it is important to analyze the expression of related genes and evaluate their correlation with clinical prognosis and potential molecular mechanisms. The TCGA (The Cancer Genome Atlas) and GEO (Gene Expression Omnibus) databases contain different tumor functional genomics datasets so that we can perform pan-cancer analysis.^1,2 Vacuolar protein sorting 13 homolog A (VPS13A) located at chr9:77177353–77421541 (GRCh38). The protein encoded by this gene controls the steps of protein circulation from the trans-Golgi network to the endosome, lysosome and plasma membrane.³ Mutations in this gene cause an autosomal recessive disease known as chorea-acanthosis. Selective splicing of this gene can lead to multiple transcriptional variations.⁴

Studies have shown that VPS13A is associated with rhabdomyosarcoma, gastric cancer, and ovarian cancer.^5,6 We summarized the distribution relationship between VPS13A and various tumors in the TCGA database, as shown in Figure 1A. In this study, for the first time, we used the TCGA project and GEO database to conduct VPS13A pan-cancer analysis (including gene expression, survival status, DNA methylation, gene mutations, immune osmosis and related cellular pathways), which can further understand the biological function of VPS13A in ubiquitin-cancer and its value in clinical diagnosis and prognosis of cancer patients.

Figure 1 The expression of VPS13A differed in different tumors and normal tissues. (A) The distribution relationship between VPS13A and various tumors in the TCGA database; (B) The expression differences of VPS13A in different tumors and normal tissues from TIMER2 database; (C) VPS13A was highly expressed in CHOL, DLBC, LAML and THYM tumors, Low expression in COAD and UCEC tumors; (D) The correlation between VPS13A expression and different tumor Pathological stages. Notes: The statistical significance computed by the Wilcoxon test is annotated by the number of stars (*p-value < 0.05; **p-value <0.01; ***p-value <0.001).

Materials and Methods

Gene Expression Analysis

In Tumor Immune Estimation Resource version 2 (TIMER2, http://timer.cistrome.org/) “Gene_DE” module, VPS13A was input to observe the difference in gene expression between tumor and adjacent normal tissue in different tumor or specific tumor subtype of TCGA project, gray columns indicate normal data are available.⁷ Some tumors do not have normal tissue or are highly restricted to normal tissues, such as adrenal cortical carcinoma (ACC), T diffuse large B-cell lymphoma (DLBC), acute myeloid leukemia(LAML), serous cystadenocarcinoma of the ovary (OV), thymic carcinoma (THYM), and so on.

We use the “analysis-box” module of Gene Expression Profiling Interactive Analysis Version 2 website (GEPIA2, http://gepia2.cancer-pku.cn/) to obtain a box map of tumor Tissue and corresponding normal Tissue Expression differences in the Genotype-Tissue Expression (GTEX) database.⁸ Set p value cutoff = 0.01, log2FC (fold change) cutoff = 1, “Match TCGA normal and GTEx data”. In addition, through the ”Pathological Stage Plot” module of GEPIA2, we obtained the violin plots of VPS13A expression in all TCGA tumors at different Pathological stages (Stage I, II, III and IV).⁹ The log2 [TPM (Transcripts per million) +1] transformed expression data were applied for the box or violin plots.

Protein Expression and DNA Methylation Analysis

UALCAN website (http://ualcan.path.uab.edu/analysis-prot.html), an interactive web resource for analyzing cancer Omics data,¹⁰ allowed us to perform protein expression analysis on the Clinical Proteomics Consortium for Cancer Analysis (CPTAC) dataset.¹¹ In this study, we input “VPS13A” and selected available datasets of six kinds of tumors, to investigate the expression level of total protein of VPS13A in primary tumor and normal tissue, respectively. The z-value represents the standard deviation of the median sample for different cancer types.

DNA methylation in primary tumors and normal tissues was analyzed in the TCGA dataset on the UALCAN website,¹² Input “VPS13A”, select available datasets of 11 kinds of tumors: infiltrating breast cancer, cervical squamous cell cancer, adenocarcinoma, cholangiocarcinoma, colorectal adenocarcinoma, esophageal cancer, squamous cell carcinoma of head and neck, clear cell carcinoma of kidney, hepatocellular carcinoma, lung adenocarcinoma, pancreatic cancer, and endometrial cancer. Beta values indicate that DNA methylation levels go from 0(unmethylated) to 1(fully methylated), [Beta-value: 0.7–0.5] for hypermethylation, and [Beta-value: 0.3–0.25] for hypomethylation. P < 0.05 indicated a statistically significant difference in DNA methylation between the primary tumor and normal tissue.

Analysis of Survival Prognosis

We used the “Survival Map” module of GEPIA2 to obtain the OS (overall survival) and DFS (disease-free survival) significance map data of VPS13A in all TCGA tumors.¹³ The threshold value of cutoff (50%) was used as the expression threshold to separate the high expression and low-expression queues. In the survival analysis module, the significance level was set to 0.05, and the log-rank hypothesis test was used to obtain the survival curve, in which P < 0.05 suggested that the difference in gene expression was related to the prognosis of patients.

Analysis of Genetic Changes

We use cBioPortal website (https://www.cbioportal.org/),^14–16 in the “Quick select” part chose “TCGA Pan Cancer Atlas Studies”, Enter “VPS13A” to query the genetic alteration characteristics of VPS13A. Change frequency, mutation type and copy number variation (CNV) of all TCGA tumors can be observed in the “Cancer Types Summary” module. The mutation site information on VPS13A can be displayed in the schematic map of protein structure or three-dimensional structure through the “Mutations” module. We also used the “Comparison/Survival” module to obtain data on the difference in overall, disease-free, progression-free, and disease-free survival of TCGA cancer cases with and without VPS13A gene alteration. Kaplan–Meier graphs with log-rank and p-value are also generated.

Immune Infiltration Analysis

We used the “Immune Genes” module of the TIMER2 Web server to explore the relationship between VPS13A expression and immune infiltration in all TCGA tumors.¹⁷ CD8+ T cells and cancer-related fibroblasts were selected, TIMER, CIBERSORT, CIBERSORT- abs, QUANTISEQ, XCELL, MCPCOUNTER and EPIC algorithms were used to test the purity Spearman rank correlation. The heat map and scatter map containing P value and partial correlation (COR) value of immune infiltration were obtained.¹⁸ COR ranges from +1 to −1 denote perfect positive to perfect negative correlations, sometimes denoted by R.

Gene Enrichment Analysis

Search the STRING website (https://string-db.org/), using a single protein name (“VPS13A”) and an organism (“homo sapiens”),¹⁹ then set the following main parameters: Minimum required interaction score [“low confidence (0.150)”], meaning of network edge meaning (“evidence-line color indicates the type of interaction evidence”), Max number of Interactors to show maximum number of interactors (“no more than 50 interactors”) and active interaction sources (“experiment”). Active interaction sources include Neighborhood in the Genome (green), Gene Fusion (purple), Co-occurrence Across Genomes (blue), Co-Expression (black), Experimental/Biochemical Data (red), Association in Curated Databases (light blue), and Co-Mentioned in PubMed Abstracts (light green).

Finally, available experimentally determined VPS13A binding proteins were obtained. We used the “similar gene” module of GEPIA2 to obtain the top 100 VPS13A-related targeted genes based on the data sets of all TCGA tumors and normal tissues.²⁰ “Correlation Analysis” module was also applied to conduct paired Pearson gene correlation analysis of VPS13A and the selected genes, a scatter plot with log2 TPM, P value and correlation coefficient (R) was obtained. The larger the R, the higher the correlation was.²¹

In addition, the “Gene_Corr” module of TIMER2 was used to input VPS13A-related genes, and purity adjusted Spearman rank correlation test was selected to perform gene correlation analysis on VPS13A and the selected genes, and heat maps and scatter maps with partial correlation (COR) and p values in different tumors were obtained. The larger R was, the higher the correlation was.²²

We used Jvenn, an interactive Venn graph viewer,²³ to perform cross analysis to compare binding proteins and related genes. KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was conducted by combining the two sets of data. Click the “Search & Color” Pathway on the “KEGG Mapper” page, and in the Color input box for the related proteins and markers, set the other parameters: Sponsorship-Specific defaults to HSA, select Include Aliases, and finally click on “Exec”. The pathways with more proteins involved were copied, and histograms were drawn using GraphPad Prism 8. Similarly, 100 Uniprots with similar genes were obtained from the UniProt website, and their colors and markers were input into the KEGG Mapper box. We upload the list of genes to David (Database for annotation, visualization, and integrated discovery) by selecting the identifier (“official”) and species (“ Homo sapiens ”) Settings. And integrated discovery), and obtained the data of functional annotation map.^24,25 Go (Gene Ontology) enrichment^26,27 was performed on these data, and copies of pathways with a large number of genes were involved. GraphPad Prism 8 was used to plot the histogram.

Results

Gene Expression Data Analysis

The expression differences of VPS13A in different tumors and normal tissues were obtained in TIMER2, as shown in Figure 1B. VPS13A was observed in Cholangiocarcinoma (CHOL), Lymphoid Neoplasm Diffuse Large B-cell Lymphoma (DLBC), Acute Myeloid Leukemia (LAML), Thymoma (THYM), Pancreatic cancer Adenocarcinoma (PAAD), Pheochromocytoma and Paraganglioma (PCPG) were highly expressed. VPS13A is low expressed in breast invasive carcinoma (BRCA), cervical squamous cell carcinoma and endocervical carcinoma (CESC), uterine corpus endometrial carcinoma (UCEC), ovarian serous cystadenoma (OV), stomach adenocarcinoma (STAD) and colon adenocarcinoma (COAD). The statistical significance computed by the Wilcoxon test is annotated by the number of stars (*: p-value <0.05; **: p-value <0.01; ***: p-value <0.001). The box chart of expression difference was obtained from GEPIA2 database. VPS13A was highly expressed (with high median) in CHOL, DLBC, LAML and THYM tumors, as shown in Figure 1C. Low expression (low median) in COAD and UCEC tumors, as shown in Figure 1C. The box chart * indicated P < 0.05, and the difference was statistically significant. We also used the “Pathological Stage Plot” module of GEPIA2 to observe the correlation between VPS13A expression and different tumor pathological stages.

Among them, ACC (P=0.02), BRCA (P=0.007), BLCA (P=0.03), LIHC (P=0.04), LUAD (P=0.03), LUSC (P=0.04), TGCT (P=0.01), THCA (P=0.005), OV (P=0.04), P<0.05 is correlated. Others have no correlation. The higher the pathological stages of ACC, BLCA, LUAD, LUSC and OV were, the higher the expression level of VPS13A was, presenting a positive correlation, while BRCA, LIHC, TGCT and THCA presented a negative correlation, as shown in Figure 1D.

Protein Expression and DNA Methylation Analysis

CPTAC dataset results showed that compared with normal tissues, total VPS13A protein was lower expressed in primary tissues of breast cancer (BRCA), serous cystadenocarcinoma of the ovary, endometrial carcinoma of the uterus, colon cancer, clear cell carcinoma of the kidney, and lung adenocarcinoma (LUAD), as shown in Figure 2A.

Figure 2 VPS13A protein expression and DNA methylation analysis. (A) VPS13A protein was lower expressed in these tumors compared with normal tissues; (B) VPS13A showed higher methylation levels in primary tumors compared with normal tissues.

TCGA dataset results showed DNA methylation analysis in primary tumors and normal tissues, VPS13A showed higher methylation levels in primary tumors of infiltrating breast cancer, cervical squamous cell cancer, adenocarcinoma, cholangiocarcinoma, colorectal adenocarcinoma, esophageal cancer, squamous cell carcinoma of head and neck, clear cell carcinoma of kidney, hepatocellular carcinoma, lung adenocarcinoma, pancreatic cancer, and endometrial cancer, as shown in Figure 2B.

Survival Analysis Data

We divided cancer cases into high and low-expression groups based on VPS13A expression level, mainly using the TCGA and GEO datasets to study the correlation between VPS13A expression and prognosis in different tumor patients with datasets. Kaplan-Meier mapping tool was used to analyze the survival data. The high expression of VPS13A gene was associated with longer overall survival (OS) of CHOL and KIRC (P=0.035, 0.0051, P<0.05). Low expression of VPS13A gene was associated with prolonged overall survival of ACC, LGG, LIHC and UCS (P=0.031, 0.0057, 0.027, 0.047, P<0.05), as shown in Figure 3A. The high expression of VPS13A gene was associated with prolonged disease-free survival (DFS) of TGCT and THCA (P=0.022, 0.025, P<0.05). Low expression of VPS13A gene was associated with prolonged disease-free survival of ACC, LIHC, Prad, and UCS (P=0.0017, 0.019, 0.022, 0.015, P<0.05), as shown in Figure 3B. Finally, high expression of VPS13A in ACC, LIHC and UCS was found to be correlated with poor prognosis.

Figure 3 Survival analysis of VPS13A in different tumors. (A) The high expression of VPS13A gene was associated with prolonged overall survival (OS) of ACC, LGG, LIHC and UCS (P=0.031, 0.0057, 0.027, 0.047, P<0.05); (B) Low expression of VPS13A gene was associated with prolonged disease-free survival (DFS) of ACC, LIHC, Prad, and UCS (P=0.0017, 0.019, 0.022, 0.015, P<0.05).

Analysis Data of Genetic Changes

We observed the genetic alteration of VPS13A in different tumor samples in the TCGA dataset. As shown in the Figure 4A, patients with endometrial cancer with “mutation” as the primary type had the highest frequency of VPS13A changes (>12%). VPS13A copy number deletions were present in all genetically altered mature B-cell tumors (2% of the time). Figure 4B shows the type, locus and number of cases of VPS13A gene mutation. We found that missense mutation was the main type of VPS13A gene mutation. We can observe the K3725fs *2 site and SHR-BD in the three-dimensional structure of SND1 protein. SHR-BD:SHR-binding domain of vacuolar-sorting associated protein 13 (2209–2455) PFAM, see Figure 4B. In addition, we explored the potential association between VPS13A gene alterations and clinical survival outcomes in patients with different types of cancer. In the cBioPortal website, through the analysis of 33 tumor mutations and prognosis, it was found that the mutation and non-mutation of VPS13A were not significantly associated with survival on the whole, which may be related to the low number of cases or lack of mutation data in many tumor samples. Figure 4C data showed that compared with UCEC patients without VPS13A change, UCEC patients with VPS13A change had better overall (P=0.0723), disease-free (P=0.315), progression-free (P=6.798e-3), and disease-specific (P=0.0486) survival outcomes. Compared with ACC patients with VPS13A change, ACC patients without VPS13A change had better overall (P=1.56e-8), progression-free (P=1.11e-16), and disease-specific (P=6.15e-9) survival outcomes.

Figure 4 VPS13A gene mutation analysis and survival analysis. (A) The genetic alteration of VPS13A in different tumor samples in the TCGA dataset; (B) The type, locus and number of cases of VPS13A gene mutation;The three-dimensional structure of SND1 protein. (C) UCEC patients with VPS13A change had better overall compared with UCEC patients without VPS13A change; ACC patients without VPS13A change had better overall compared with ACC patients with VPS13A change; (D) The phosphorylation site and sequence information of VPS13A in the PhosphoNet database.

Protein Phosphorylation Data Analysis

Phosphonet and qPhos databases were also used to analyze phosphorylation levels. The qPhos database result interface contains the protein number in UniProt database, the gene name, the position of phosphorylated amino acids in the protein, sequence information, sample name, sample type, experimental conditions, detection method, phosphorylation difference results and significance.^28,29 qPhos database showed that under PaSC standard, TMT detection method was used to obtain the phosphorylation site and sequence information on VPS13A,³⁰ and the difference was significant (P=0.017). The PhosphoNet database shows that VPS13A has 240 phosphorylated sites, with more phosphorylated sites at S839, as shown in Figure 4D.

Immune Infiltration Analysis Data

Tumor-infiltrating immune cells, as an important part of the tumor microenvironment, are closely related to the occurrence, progression or metastasis of tumors.^31–34 It has been reported that tumor-associated fibroblasts in the tumor microenvironment stroma are involved in regulating the functions of various tumor-infiltrating immune cells.^35,36

We used TIMER, CIBERSORT, CIBERSORT- abs, QUANTISEQ, XCELL, MCPCOUNTER and EPIC algorithms to investigate the potential relationship between different levels of immune cell infiltration and VPS13A gene expression in different types of TCGA cancer. Cor 0–1 is a positive correlation, which is red; Cor 0- (−1) is a negative correlation, which is blue. After a series of analyses, based on all or most of the algorithms, we found that there was a statistically positive correlation between the immune infiltration of CD8+T cells and the expression level of VPS13A in DLBC, LUAD, SKCM, TGCT and UVM metastases (Figure 5A). We observed a positive correlation between VPS13A expression and cancer-associated fibroblasts in BRCA, CESC, LIHC, and THYM, but a negative correlation in BLCA, KIRC, and TGCT (Figure 5B). According to TIDE algorithm, the expression level of VPS13A in ACC is positively correlated with the infiltration level of cancer-related fibroblasts (Cor = 0.262, P=2.49e-02), P < 0.05, and the COR > is positively correlated, which is of statistical significance, as shown in Figure 5C.

Figure 5 Immune infiltration analysis. (A) The correlation between the immunoinfiltration of CD8+T cells and the expression level of VPS13A in DLBC, LUAD, SKCM, TGCT and UVM metastases; (B) The correlation between VPS13A expression and cancer-associated fibroblasts in BRCA, CESC, LIHC, and THYM, but a negative correlation in BLCA, KIRC and TGCT; (C) The expression level of VPS13A in ACC is positively correlated with the infiltration level of cancer-related fibroblasts (Cor =0.262, P=2.49e-02).

Enrichment Analysis of VPS13A Related Complexes

In order to further study the molecular mechanism of VPS13A gene in tumor genesis, we tried to screen the targeted VPS13A binding protein and the genes related to VPS13A expression, and performed a series of pathway enrichment analysis. We used the GEPIA2 tool combined with all tumor expression data of TCGA to obtain the top 100 genes associated with VPS13A expression. Paired Pearson gene correlation analysis was conducted for VPS13A and the selected genes, and a scatter plot of correlation analysis was obtained. The larger R was, the higher the correlation was, as shown in Figure 6A. We used the Gene_Corr module of TIMER2 to provide VPS13A-related genes, and obtained heat maps and scatter maps with partial correlation (COR) and p values for different tumors, as shown in Figure 6B. Among them, GCC2 (R=0.66, P=0), ZnF518A (R=0.64, P=0), INO80D (R=0.63, P=0), NHLRC2 (R=0.63, P=0), SREK1 (R=0.62, P=0) and ZnF782 (R=0.62, P=0) were all correlated with VPS13A in different tumors, and R>0 was indicated as a positive correlation, P<0.05, which was statistically significant. Based on the STRING tool, we obtained a total of 50 binding proteins, all of which were supported by experimental evidence. Figure 6C shows a diagram of these protein network interactions. The line between nodes represents the type of interaction between two proteins, and the red represents Experimental Data. The color of protein node is different from the score value. Among them, the score between PMSA1 and PMSA4 was the highest (score=0.994), the score between ULK3 and VPS13A was the lowest (score=0.157).

Figure 6 Enrichment analysis of VPS13A related complexes. (A) Scatter plot of correlation analysis was conducted for VPS13A and the selected genes; (B) The partial correlation (COR) and p values between VPS13A-related genes in different tumors; (C) Protein- Protein Interaction network.

VPS13A Signaling Pathway Analysis

We combined the two datasets for KEGG enrichment analysis and obtained the signaling pathway diagrams of 50 VPS13A binding proteins and 100 VPS13A associated genes. RAB7A is one of the VPS13 family of acting proteins, can regulate autophagosome-lysosomal fusion and autosomal lysosomal maturation.³⁷ In the neurodegenerative pathway and the amyotrophic lateral sclerosis pathway, protein enrichment is the largest, followed by the proteasome pathway, as shown in Figure 7A. The associated genes involved in RNA degradation, autophagy, cell senescence, cell cycle, apoptosis, TGF-β, cAMP, mTOR, Hippo, cGMP-PKG, Wnt, oxytocin, estrogen, metabolic pathway, cancer and other 66 signaling pathways. GO enrichment analysis was conducted on the annotation data of associated genes obtained by David, and it was found that these genes were related to DNA, RNA and protein synthesis, as well as DNA binding, DNA repair, cell division, protein binding and protein phosphorylation, as shown in Figure 7B. VPS13A has a wide range of effects in vivo and is very important for cell survival. Jvenn was used for cross analysis of related genes and binding proteins, and no obvious common members were found, as shown in Figure 7C.

Figure 7 VPS13A signaling pathway analysis. (A) The signaling pathway diagrams of 50 VPS13A binding proteins; (B) GO enrichment analysis was conducted on the 100 VPS13A associated genes; (C) Cross analysis of related genes and binding proteins.

Discussion

VPS13A has been shown to be associated with rhabdomyosarcoma, gastric cancer and ovarian cancer. VPS13A refers to vacuolar protein sorting-associated protein13A, which is a peripheral membrane protein located in the close contact with the endoplasmic reticulum (ER) and mitochondria (MT), and is associated with mitochondria, endoplasmic reticulum and lipid droplets. VPS13A interacts with the endoplasmic reticulum ER resident protein VAP-A through its FFAT domain and with mitochondria through its C-terminal domain.³⁸ The protein encoded by this gene can regulate intracellular protein transport and induce autophagy by mediating protein endocytosis and lysosomal degradation. Defects in the autophagy pathway can be caused by VPS13A deficiency.³⁹ VPS13A can regulate apoptosis by activating the Akt-PI3K pathway and increases cell survival rate.⁴⁰ VPS13A can affect the morphology and function of mitochondria in cells and regulate apoptosis and autophagy, which is very important for cell survival. However, the pan-cancer analysis of VPS13A is still lacking, so for the first time, we investigated the distribution characteristics, expression and prognosis of VPS13A in 33 tumor datasets of TCGA and GEO.

In some tumors, the expression of VPS13A gene was positively correlated with the level of cancer-associated fibroblasts, indicating that VPS13A was associated with the occurrence of cancer, and the higher the pathological stage of some tumors, the higher the expression of VPS13A gene was. Chorea levels of VPS13A-expressed protein in blood samples can be clinically screened to predict the occurrence of cancer and determine the pathological stage of cancer.

The methylation level of VPS13A is higher in the cells of tumor patients, suggesting that VPS13A methylation may be one of the causes of cancer. Therefore, the level of gene methylation can be detected to screen for cancer, and gene demethylation can be used to prevent the occurrence of cancer through treatment. The expression level of VPS13A gene was correlated with the infiltration level of CD8+T cells and other immune cells, indicating that the increased expression level of VPS13A was associated with increased cancer cells and enhanced immune response. Because VPS13A is highly expressed in most tumors, it can reduce the expression level of VPS13A mRNA or protein in patients’ cells, reduce cancer-related fibroblasts, activate tumor cell apoptosis pathway, reduce autophagy level, reduce immune response, promote cell death, improve the therapeutic effect, and prolong the survival time. For some cancer patients, screening for VPS13A mutations in vivo can predict the prognosis of patients. For example, promoting VPS13A mutations in endometrial cancer cells can improve the prognosis of patients.

Conclusion

This time, we analyzed the changes of VPS13A in different tumors from the aspects of gene expression and clinical prognosis, DNA mutation, DNA methylation, protein phosphorylation, and immune cell infiltration. Most of the data were statistically significant, which helped to understand the role of VPS13A in various tumors and to observe the occurrence of tumors from the perspective of clinical tumor samples. However, the sample size of this study is insufficient, and there may be some deviation in the data. It is necessary to increase the clinical data and conduct clinical trials for further verification.