Bioinformatics Analysis Of Gene Expression Changes Conferring Drug Resistance In Tumor Samples From Melanoma Patients With BRAF Mutations

BackgroundMelanoma is a highly malignant tumor of the skin melanocytes that is insensitive to chemo/radiotherapy. Most patients with melanoma have a poor prognosis, and the disease is associated with a high mortality rate. Patients with melanoma frequently exhibit activation of the BRAF gene due to somatic mutations, with up to 50% of patients having BRAF oncogenic mutations. The most common BRAF mutation in melanoma, BRAF V600E, accounts for about 79% of BRAF mutations. BRAF is a member of the RAF kinase family that includes ARAF, BRAF, and CRAF. BRAF mutations can lead to constitutive activation of downstream signaling through the MAPK/MEK-ERK pathway, cell migration and proliferation, and tumor growth.Good clinical outcomes have been obtained with melanoma treatments targeting BRAF mutations and using MEK/ERK inhibitors. As a first-line clinical treatment for melanoma, vemurafenib (PLX4032) is a potent inhibitor of mutated BRAF and a sensitive therapy for advanced melanoma. However, targeted inhibitors usually maintain their efficacy for only 8 to 9 months. Thereafter, the patient can develop resistance to the inhibitor, and the tumor can achieve rapid growth. Thus, controlling drug resistance is a key issue in melanoma treatment.Several studies have attempted to elucidate the mechanisms of drug resistance in melanoma patients. Hepatocyte growth factor (HGF) has been found in the stromal cells of patients carrying BRAF mutations, and a correlation has been demonstrated between HGF-secreting stromal cells and resistance to Raf inhibitors. Another potential cause of resistance development is COT/Tpl2 overexpression by cells Mutations of the NARS gene were found in vemurafenib-resistant cells in vitro and in the lymph node cells of vemurafenib-resistant patients. NARS mutations may be active in the MAPK pathway and may lead to resistance of melanoma cells to targeted inhibitors.The BRAF mutation itself could lead to resistance development. In previous reports, patients with the BRAF V600E mutation had a poor prognosis due to acquired resistance to vemurafenib and trametinib. Other studies showed that the BRAF V600E mutation or MEK inhibitor resistance may be related to EGFR activation in tumor cells. For example, Prahallad et al. found that in a subset of patients, BRAF V600E inhibitors may lead to EGFR activation, which, in turn, may enhance the resistance of cancer cells to BRAF inhibitors. Sun et al. found that EGFR expression enhances the proliferation of melanoma cells in the presence of inhibitors against BRAF or MEK.There is still no effective clinical treatment for patients with EGFR-activating feedback. Moreover, it is unknown whether patients with EGFR-activating BRAF V600E mutations experience gene expression changes between tumor sites before and after disease or treatment. Therefore, the aim of this study was to analyze changes in the expression of genes by melanoma tumors in patients with EGFR-activating BRAF mutations and during the development of drug resistance. The overall goal was to identify potential drug targets for melanoma treatment-resistant patients.Materials and MethodsIdentification of differentially expressed genes (DEGs)from a public databaseWe downloaded sample numbers SRR961663, SRR961664, SRR961665, SRR961666, SRR961667, and SRR961668 of GSE50535 from the Gene Expression Omnibus (GEO). We used TopHat and cufflinks software packages to analyze assembly data and differences in the gene expression profiles, respectively. From the TCGA Cutaneous Melanoma Database, we downloaded RNA-seq data of samples, taken from tumor and nontumor sites, from patients with BRAF V600E mutations of unknown EGFR activation status. Using Edger, we analyzed DEGs between pre-vs. post-treatment data from the GEO site and between pre-treatment data from the GEO site vs. tumor or nontumor sample data from the TCGA site. A P value less than 0.05 was considered statistically significant.Functional enrichment analysis of DEGsWe performed a functional enrichment analysis of DEGs, including geneontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, using the Database for Annotation Visualization and Integrated Discovery (DAVID). GO analysis included the categories of cellular component (CC), biological process (BP), and molecular function (MF). For all analyses, differences with a P value less than 0.05 were regarded as statistically significant.ResultsGene expression changes in tumors before and after treatment for drug-resistant patients with BRAF mutationsWe analyzed six cases of GSE50535 melanoma from three patient samples before and after treatment. This analysis resulted in the identification of 94 significant DEGs (62 up- and 32 down-regulated genes), which we clustered through GO analysis in DAVID.In the CC ontology, most of the significantly enriched genes were associated with the terms plasma membrane (PM) (30 genes,33.25% of all DEGs; P=0.002), PM part (21 genes,22.58%; P= 0.002), and sarcolemma (4 genes,4.30% of all DEGs; P=0.003):BGN, COL6A3, COL6A2, and CACNA1C; P= 0.003). In the BP ontology (Fig. 1(b)), nine genes (PTPRC, CORO1A, IKZF1, CLU, CD4, INPP5D, C1QC, CD74, and HLA-DRA) were enriched in the positive regulation of immune system process (P=1.08×10). Other categories with significant enrichment included immune response (P= 0.003), positive regulation of lymphocyte activation (P= 7.29×10), and lymphocyte differentiation (P= 9.84×10). Eight genes were enriched in the positive regulation of cell differentiation (P= 7.29×10). In MF ontology (Fig. 1(c)), categories with significant enrichment of DEGs included extracellular matrix structural constituent (5 genes:BGN, ELN, COL6A2, CD4, and COL5A1; P= 5.43×10), glycoprotein binding (4 genes; P= 5.61×10), kinase binding (6 genes: PTPRC, CORO1A, ITGB2, CD4, TNNI3, and TOP2B; P= 0.001), and protein kinase binding (5 genes:PTPRC, ITGB2, CD4, TNNI3, and TOP2B; P= 0.003).We also performed a KEGG pathway analysis of the 94 DEGs The most significantly enriched KEGG pathway was cell adhesion molecules (CAMs). Five DEGs were enriched in the antigen processing and presentation pathway.Gene changes in tumor compared to nontumor melanoma samplesWe sought to gain a better understanding of the differences between tumor samples from melanoma patients with BRAF mutations and nontumor samples from melanoma patients. We compared SRR961663, SRR961665, and SRR961667 melanoma samples from three patients before treatment to a nontumor sample from the TCGA Cutaneous Melanoma Database. We identified 274 significant DEGs, which we then clustered through GO analysis in DAVID.In the CC ontology, the most significantly enriched categories were adhesion-related items pertaining to the ribosome, including the cytosolic ribosome (27 genes; P= 3.27×10), ribosome (34 genes; P= 2.59×10), ribosomal subunit (26 genes; P= 2.59×10), cytosolic small ribosomal subunit (15 genes; P= 3.38×10), MHC protein complex (11 genes; P= 1.03×10), and MHC class I protein complex (6 genes; P= 5.25×10). In the BP ontology (Fig.2(b)), most categories with significant DEG enrichment were associated with the translation process, including translational elongation (33 genes; P= 1.28×10) and translation (37 genes; P= 1.50×10). DEGs were also enriched in immune-related categories, such as immune response (39 genes; P= 1.08×10), antigen processing and presentation (15 genes; P=5.67×10), and antigen processing and presentation of peptide antigens (9 genes; P= 2.59×10). In the MF ontology, categories with significant enrichment in DEGs included ribosome structural constituent (31 genes; P= 4.67×10), structural molecule activity (40 genes; P= 1.30×10), RNA binding (27 genes; P= 1.73×10), MHC class II receptor activity (6 genes; P= 3.43x 10), and MHC class I receptor activity (5 genes; P= 5.64x 10).To gain a better understanding of the function of genetic differences, we performed pathway analysis of these 274 DEGs. KEGG pathway terms with significantly enriched DEGs included the ribosome (30 genes; P= 6.20×10), antigen processing and presentation (17 genes; P= 4.64×10), and systemic lupus erythematosus (14 genes; P= 4.09×10).DEGs between tumor and nontumor samples from melanoma patientsNine genes(C1QC, CADPS, CD74,CLU, CORO1A, FMN1, HLA-DPA1, HLA-DRA, and LSP1) were differentially expressed in both DEG analyses described above (before vs. after treatment and tumor vs. nontumor samples). Greatest enrichments in DGEs were found in the GO category of immune process (P-1.21×10) and the KEGG pathway of antigen processing and presentation (P= 7.81×10).DEGs between BRAF V600E mutation patients with and without EGFR activationFinally, we analyzed DEGs between samples taken from tumor and nontumor sites of patients with the BRAF V600E mutation of unknown EGFR activation status from the TCGA Cutaneous Melanoma Database. Twenty-seven DGEs changed in both groups. GO analysis by DAVID revealed that nine genes were enriched in the immune response category (P= 1.22X10), and eight genes were enriched in the defense response category (P= 7.96×10).ConclusionWe analyzed gene expression changes in the tumors of patients with EGFR-activating BRAF V600E mutations before and after BRAF inhibitor treatment. We also compared samples from taken from tumor vs. nontumor sites in patients with BRAF V600E mutations of unknown EGFR activation status. We analyzed the functions of DEGs by GO annotation and KEGG pathway enrichment analyses. Through this study, we hope to provide information to guide the development of new therapeutic strategies for melanoma patients with EGFR activation who are resistant to typical drugs.RNA-seq technology provides a powerful tool for analyzing gene expression. Using RNA-seq data for melanoma patients with EGFR activation, we identified 94 genes that were differentially expressed between before and after treatment, including 62 up-and 32 down-regulated genes. Gene functional annotation revealed 30 genes related to membranes in CC oncology, which provides a possible direction for future studies.Ion channels on the membrane are involved in numerous tumor cell activities, including cell proliferation, differentiation, secretion, survival, and others. We found significant upregulation of CACNA1C(P=5.00×10), which encodes the a-1 subunit of a voltage-dependent calcium channel located on the PM. This gene not readily detectable in normal tissues. Calcium channel proteins have been associated with primary tumors of the colon, lung, and skin. Some drugs targeting CACNA1C, such as magnesium sulfate and nicardipine, have been reported. Magnesium sulfate can be used to inhibit the action potential of muscle cells, thereby reducing the frequency and strength of contractions. Nicardipine is a strong calcium channel inhibitor with important vasodilatatory and antihypertensive characteristics, which can be used to enhance the efficacy of certain antitumor agents.We also found a significant downregulation of SLC4A10 (P=5.005×10), which belongs to a small family of sodium-coupled bicarbonate transporters that regulate the intracellular pH of neurons. In addition to unlimited cell proliferation, cancer is characterized by an altered cellular environment that promotes tumor cell proliferation and metastasis. The interior pH of cancer cells is alkaline (pH≥7.4), whereas tumor cells are more acidic than normal cells in the extracellular environment. Downregulation of SLC4A10 leading to decreased Cl/HCO3 transport may be related to environmental alterations for melanoma growth.We identified 274 genes that were differentially expressed between untreated EGFR-activated melanoma samples and nontumor samples from the TCGA Database. GO and pathway analyses revealed the enrichment of many genes in immune-related processes, especially antigen processing and presentation processes. One enriched gene was CD74, which encodes a protein associated with the class Ⅱ major histocompatibility complex (MHC) and is an important chaperone that regulates antigen presentation. CD74 serves as a cell-surface receptor for the cytokine macrophage migration inhibitory factor, which, when bound to the encoded protein, initiates survival pathways and cell proliferation. Prior studies have indicated that CD74 is only expressed in melanoma cells and not in benign melanocytes Milatuzumab is a drug used for the treatment of tumors expressing the CD74 antigen.Antigen presentation plays a key role in the development of melanoma vaccines. Broadly speaking, a tumor cell is also an antigen-presenting cell. Tumor cells form a complex with MHC class Ⅰ molecules by the cytosolic processing pathway, and tumor antigens on the cell surface are recognized by CD8T cells. Alternatively, tumor cells can form MHC class Ⅰ or Ⅱ molecules by lysosomal processing of tumor antigens from dendritic cells or specialized antigen-presenting macrophages. These tumor cells are then recognized by CD8/CD4T cells.We found significant enrichment in many genes related to ribosomal processes. Ribosome synthesis and translational control are essential processes for cells. Several tumor suppressor genes and proto-oncogenes can affect the formation or modification of ribosomes. However, the mechanisms by which these genes affect ribosomes remain unclear at present, and further experiments are needed.We analyzed tumor and nontumor samples from patients with BRAF V600E mutations of unknown EGFR status from the TCGA Database. We compared these results with DEGs between tumor samples from patients with EGFR-activating BRAF V600E mutations before and after treatment. Twenty-seven genes were changed in both comparisons. These genes were mainly enriched in categories related to immune response, which suggests that using immunotherapy in the early stages of melanoma may be a useful therapeutic approach. Nine genes were changed between the comparisons of tumor vs. nontumor tissues and treated vs. untreated tumors. Most of these genes had different changes in different stages; for example, a gene that was downregulated in the tumor area compared to the nontumor area might be upregulated upon acquiring resistance after treatment. DEGs that were significantly enriched in the immune process category by GO analysis included CLU, C1QC, CD74, and HLA-DRA. The C1QC gene is the target of many drugs, as the C1QC protein is an important component of human complement. Studies have associated a lack of C1QA with lupus and glomerulonephritis. Several C1QC-targeting drugs, such as tositumomab, palivizumab, and cetuxima, are used clinically to treat cancer. CLU has been shown to be related to some biological events, such as apoptosis and tumor development, as well as neurodegenerative diseases.However, our study has many deficiencies. For example, the sample size was relatively small; e.g., there were only three samples with EGFR activation and only one nontumor sample. The small sample size may lead inaccuracies when comparing the tumor vs. nontumor samples. Therefore, to explore the most common genetic changes in melanoma patients, we not only compared SRR961663, SRR961665, and SRR961667 melanoma samples from three patients before treatment to a nontumor sample from the TCGA Cutaneous Melanoma Database, but also analyzed tumor and nontumor samples from patients with BRAF V600E mutations of unknown EGFR status from the TCGA Database. We compared these results with DEGs between tumor samples from patients with EGFR-activating BRAF V600E mutations before and after treatment. We tried to improve the accuracy of these findings and will continue to collect samples with the selected mutations for further research.The treatment of melanoma is a complex process. We found significant changes in genes related to the PM of BRAF inhibitor-resistant melanoma patients with EGFR-activated tumors before and after treatment. Significant changes in immune process-related genes were also found in melanoma patients between tumor and nontumor samples. Although these findings provide potential directions for clinical melanoma-specific therapy, follow-up studies of melanoma are needed.