Bioinformatics Analysis of Differential Gene and MicroRNA Expression in Lung Adenocarcinoma: Genetic Effects on Patient Prognosis, as Indicated by the TCGA Database.

Objective: To investigate the differential expression of genes and microRNAs (miRNAs) in patients with lung adenocarcinoma and the relationship between such changes and patient prognosis.
Methods: We analyzed the expression levels of genes and miRNAs in lung adenocarcinoma tissues and adjacent normal tissues using The Cancer Genome Atlas database (TCGA). We analyzed the function of the differentially expressed genes and miRNAs in a co-expression network. Finally, we performed survival analysis of differential genes and miRNAs in the co-expression network using clinical data from the TCGA database.
Results: We successfully identified 6064 differentially expressed genes: 5324 upregulated genes and 740 downregulated genes. And we identified 161 differentially expressed miRNAs: 126 upregulated miRNAs and 35 downregulated miRNAs. We identified several genes that were related to each other in the co-expression network. Further analysis revealed that the high expression levels of G6PC, APOB, F2, PAQR9, and PAQR9-AS1 genes were associated with poor prognosis. However, there was no significant correlation between the expression of hsa-mir-122 with regards to patient prognosis. Conclusions: Our data showed that hsa-mir-122 and a number of related genes may affect the prognosis of patients with lung adenocarcinoma by regulating the cytoskeleton, thus promoting tumor angiogenesis and the metastasis of tumor cells. The high expression levels of some differentially expressed genes was associated with the low survival rate in patients with lung adenocarcinoma. However, the levels of hsa-mir-122 were not correlated with patient prognosis.

Introduction

Of all forms of malignancy, lung cancer is associated with the highest rates of morbidity and mortality. Non-small cell lung cancer (NSCLC) accounts for 85% of lung cancer patients. The most common form of non-small cell lung cancer is lung adenocarcinoma (LUAD); this particular form accounts for approximately 40% of all lung cancers.[2-5] With the use of advanced screening equipment, and the improvement of techniques that can be used to prevent and diagnose NSCLC at an early stage, the prognosis of patients with this particular form of cancer is gradually improving.[6,7] However, most patients tend to be diagnosed when they are in an advanced stage of the disease. Without the possibility of surgery, the only options for such patients are radiotherapy and chemotherapy. Tumor development is a complex process that involves a multitude of genes and different stages of development. Consequently, it is very important that we investigate oncogenes and miRNAs as these play critical roles in the pathogenesis, targeted therapy, and prognostic evaluation, of lung cancer. The Cancer Genome Atlas (TCGA) is a publicly funded project and aims to provide a comprehensive cancer genome map, to allow researchers to explore cancer-related knowledge, to help us to understand the mechanisms underlying cancer, and to facilitate diagnosis or therapeutic interventions. Based on this database, we investigated the expression levels of key genes and related miRNAs in patients with lung adenocarcinoma to ascertain whether these molecules exerted effect on prognosis. Our ultimate aim was to identify new biomarkers for lung adenocarcinoma that may facilitate clinical diagnosis and treatment.

Methods

The Acquisition of Gene Transcriptome Samples

We used the TCGA database (https://cancergenome.nih.gov/) to acquire transcriptomic data (Level 3) from patients with lung adenocarcinoma. A number of different factors were used to screen the database, as follows: disease type (adenomas and adenocarcinomas; project ID (TCGA-LUAD), workflow types (HTSeq-Counts), data type (transcriptome profiling), and experimental strategy (RNA-Seq).

The Acquisition of MicroRNA Sequences

Sequence data (Level 3) relating to miRNAs were acquired for patients with lung adenocarcinoma from the TCGA (https://cancergenome.nih.gov/). The screening conditions were as follows: disease type (adenomas and adenocarcinomas); Project ID (TCGA-LUAD), data typ (miRNA expression quantification), and experimental strategy (miRNA-Seq).

Statistical Methods

We chose to use the edgeR package (in R3.6.2 software, URL: https://www.R-project.org/) to standardize and analyze the data acquired from the TCGA database (https://cancergenome.nih.gov/). We used strict screening conditions (|logFC| ⩾ 2; corrected P value <.05) to identify differentially expressed genes and miRNAs. We use the gplots software package (version 3.0.1.1; https://CRAN.R-project.org/package=gplots) to create volcanic and heat maps. Next, we calculated suitable soft thresholds in the WGCNA package (version: 1.69; https://CRAN.R-project.org/package=WGCNA) in order to create a co-expression matrix for differentially expressed genes and miRNAs. Next, we created a co-expression network for the differentially expressed genes and miRNAs using Cytoscape 3.7.2 (https://cytoscape.org/). software. Functional enrichment analysis, was carried out using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses and data were presented using the ggplot2 package (Version: 3.3.2; https://CRAN.R-project.org/package=ggplot2). Next, we used TargetScan (http://www.targetscan.org; http://www.targetscan.org/vert_72/), miRDB (http://mirdb.org/), PicTar (https://pictar.mdc-berlin.de/), and miRanda (http://www.microrna.org/microrna/getGeneForm.do) to predict the target genes for the miRNA we chose to analyze. Meanwhile, we used a Venn Diagram to visualize the predictive results arising from the target genes of differentially expressed miRNAs within the co-expression network. Finally, we performed survival analysis of differentially expressed genes and miRNAs within the co-expression network using log-rank test methods and clinical data from patients. P < .05 was considered to indicate statistical significance.

Results

Genes Showing Differential Expression When Compared Between Lung Adenocarcinoma and Adjacent Normal Tissues

We downloaded, collated, and then analyzed, transcriptome data from the TCGA database. A total of 551 transcriptome samples were identified, including 497 adenocarcinoma tissue samples and 54 para-cancerous normal tissue samples. According to our selection conditions (|logFC| ⩾ 2 and a corrected P value <.05), we were able to use the edgeR package in R3.6.2 software to identify differentially expressed genes between cancerous and adjacent normal tissues. We identified 6064 differentially expressed genes, including 5324 upregulated genes and 740 downregulated genes (Figure 1).

Figure 1.

Gene expression in lung adenocarcinoma tissues and adjacent normal tissues: (A) heat map of differentially expressed genes between cancer tissues and adjacent normal tissues. The horizontal axis represents sample identifiers. While the vertical axis refers to gene names. Red represents high gene expression levels, while green indicates low expression levels and (B) volcanic map of gene expression levels in lung adenocarcinoma. The horizontal axis represents −log10(FDR), while the vertical coordinate represents logFC. The red points represent upregulated genes while the green points represent downregulated genes.

miRNAs Showing Differential Expression When Compared Between Lung Adenocarcinoma and Adjacent Normal Tissues

We downloaded, collated, and then analyzed, miRNA data from the TCGA database. We acquired a total of 528 miRNA sequences, including 483 adenocarcinoma tissue samples and 45 para-cancerous normal tissue samples. According to our selection conditions (|logFC| ⩾ 2 and a corrected P value <.05), we were able to use the edgeR package in R3.6.2 software to identify differentially expressed miRNAs between cancer and adjacent normal tissues. Finally, we identified 161 differentially expressed miRNAs, including 126 upregulated miRNAs and 35 downregulated miRNAs (Figure 2).

Figure 2.

miRNA expression in lung adenocarcinoma tissues and adjacent normal tissues: (A) heat map of the differentially expressed miRNAs between lung adenocarcinoma and adjacent normal tissues. The horizontal axis represents the sample identifiers, while the vertical axis represents the miRNAs. Red represents high miRNA expression levels while green represents low levels of expression and (B) volcanic map showing miRNA expression in lung adenocarcinoma. The horizontal axis represents −log10(FDR) while the vertical axis represents logFC. The red points represent upregulated miRNAs while the green points represent downregulated miRNAs.

Construction of a Co-expression Network Relating to Differentially Expressed Genes and miRNAs

According to the levels of differentially expressed genes and miRNAs, we were able to use the WGCNA package and topological network analysis to set a soft threshold for a co-expression matrix between differentially expressed genes and miRNAs. It was important that the link between genes and miRNAs in the network diagram followed the distribution of a non-scaled network as this would mean that the co-expression network would be more meaningful in many biological aspects. To do this, we set a soft threshold of 18 with which to create a co-expression matrix between differentially expressed genes and miRNAs. Based on this co-expression matrix, the co-expression network between differentially expressed genes and miRNAs was created using Cytoscape 3.7.2 software (Figure 3). Connections between nodes in the network represent the presence of an association related to expression.

Figure 3.

A co-expression network map of differentially expressed genes and miRNAs between lung adenocarcinoma and adjacent normal tissues. The green nodes represent differentially expressed genes while the red nodes represent miRNAs. The lines between the genes and miRNAs indicates a direct correlation. (A) A co-expression network diagram showing up-regulated hsa-mir-122 and its associations with differentially expressed genes. The lines in the network indicate a quantitative correlation between the genes and miRNAs. All of the genes shown in Figure 3A were upregulated in lung adenocarcinoma tissue and correlated directly with high expression levels of hsa-mir-122. (B-F) The co-expression network diagrams showing differential microRNAs and genes between lung adenocarcinoma and adjacent normal tissues.

Figure 3 shows a range of differentially expressed genes, including VTN, ALB, PAQR9, AFP, CREB3L3, PRODH2, AC080128.1, HRG, APOB, PAQR9-AS1, LGALS14, and G6PC. Only hsa-mir-122, when upregulated, showed a significant correlation with these differentially expressed genes (Figure 3A). Previous studies have shown that hsa-mir-122 is associated with the development of several different types of cancer, including breast cancer, renal cancer, colorectal cancer,[10,11] and gastric cancer. For instance, the overexpression of hsa-mir-122 in patients with gastric cancer has been shown to inhibit the proliferation, metastasis, and invasion of gastric cancer cells, while the downregulation of hsa-mir-122 promotes the process of cell proliferation. However, the overexpression of hsa-mir-122 in patients with colon cancer has been shown to promote development, metastasis, and invasion of tumors. In order to study the role of hsa-mir-122 and its associated genes in the pathogenesis of lung adenocarcinoma further, we selected hsa-mir-122 and several related genes (Figure 3A) in our co-expression network. This network was created using differentially expressed genes and miRNAs between lung adenocarcinoma and adjacent normal tissues.

Functional Enrichment Analysis of Differentially Expressed Genes Related to hsa-mir-122

We used the DAVID website (https://david.ncifcrf.gov/)[13,14] to perform GO (gene ontology) enrichment analysis of genes in our co-expression network diagram of hsa-mir-122 and its associated differentially expressed genes (Figure 3A). Of the differentially expressed genes shown in Figure 3A, APOB, ALB, and VTN were all significantly enriched in 1 particular biological process: receptor-mediated endocytosis (GO:0006898). We also found that ALB, HRG, and VTN were all significantly enriched in a certain cell component: blood cell microparticles (GO:0072562). APOB, HRG, and VTN genes were significantly enriched in 1 particular molecular function: binding heparin (GO:0008201) (Figure 4A). Next, we used the KOBAS website (kobas.cbi.pku. edu.cn/) to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. This showed that the pathway in which CREB3L3, G6PC, and VTN were most significantly enriched was the PI3K-Akt signaling pathway (Figure 4B).

Figure 4.

Enrichment analysis of differentially expressed genes related to hsa-mir-122: (A) data arising from gene ontology functional enrichment analysis of differentially expressed genes related to hsa-mir-122. The horizontal axis refers to FDR while the vertical axis refers to GOTERM_NAME. Different shapes represent biological processes, cellular components, and molecular functions, respectively. The circles refer to biological processes, the triangles refer to cellular components, while the squares refer to molecular functions. The size of the shapes represents the number of genes enriched on the gene ontology term. Coloration reflects −log10(FDR); the FDR value of red is less than green and (B) results arising from KEGG pathway analysis of differentially expressed genes related to hsa-mir122. The bar chart shows the results of KEGG pathway analysis. The vertical axis refers to the names of pathways that were enriched by differentially expressed genes associated with hsa-mir-122. The length of the bar represents the number of genes enriched in the pathway. Coloration reflects −log10(FDR); the FDR value of blue is less than red.

The Prediction and Functional Analysis of Target Genes for hsa-mir-122

We Used TargetScan, miRDB,[18,19] PicTar,[20,21] miRanda[22,23] to predict the target genes for hsa-mir-122. The miRanda website predicted 1872 target genes, while the PicTar, TargetScan, and miRDB, websites predicted 115, 1129 and 516 target genes, respectively. Analyses showed that 27 target genes were consistently found on all 4 database: ALDOA, CCNG1, GIT1, CS, DDR2, G3BP2, IQGAP1, STX16, DICER1, FOXP2, P4HA1, PAK3, MIPOL1, NPEPPS, OCLN, BACH2, DR1, SMYD4, EPO, MOSPD1, BRPF1, CPEB1, MASP1, MAF1, LRP10, LAMC1, and FUNDC2 genes. These 27 target genes were used for functional analysis to further explore the molecular function of hsa-mir-122 in lung adenocarcinoma. Results showed that these target genes were significantly enriched in certain GO annotated entries. For biological process, there was significant enrichment for the negative regulation of transcription by RNA polymerase II (GO:0000122); this refers to any biological process that prevents or attenuates the frequency, rate, or extent of transcription mediated by RNA polymerase II. For cellular component, there was significant enrichment for the nucleus (GO:0005634). Chromosomes are preserved and copied in the membrane-bound organelles of eukaryotic cells. In most cells, the nucleus contains chromosomes other than organelle chromosomes, in which RNA is synthesized and processed. In some species or special types of cells, the process of RNA metabolism or DNA replication may not occur in the nucleus. For molecular function, there was significant enrichment for protein domain specific binding (GO:0019904). This refers to interactions with specific domains of proteins that occur in a selective and non-covalent manner. Of these target genes, KEGG pathway analyses showed that the most significant enrichment for the GIT1, IQGAP1, and PAK3 genes was the actin cytoskeleton pathway.

Prognostic Analysis of Survival in Patients With Lung Adenocarcinoma

We downloaded clinical data from patients with lung adenocarcinoma from the TCGA database; 486 samples were acquired. We carried out survival analysis for hsa-mir-122, and the differentially expressed genes that are associated with hsa-mir-122 using the survival package in R3.6.2 software. We divided the patients into 2 different groups according to the median expression level. In other words, patients with expression levels that were higher than the median were placed into the high expression group, while those with expression levels that were lower than the median were placed into a low expression group. Comparative analysis of the survival status of the 2 groups showed that the expression levels of G6PC, APOB, F2, PAQR9, and PAQR9-AS1 genes exerted significant effects on the survival prognosis of patients and that low expression levels are associated with a better prognosis (P < .05). However, there was no significant correlation between the expression of hsa-mir-122 and the prognosis of patients (Figure 5).

Figure 5.

Kaplan-Meier survival curve for patients with lung adenocarcinoma. The horizontal axis represents the survival time (years) while the vertical axis represents survival rate (also known as cumulative survival probability or survival function).

Discussion

Lung cancer is considered to be the leading cause of cancer-related death. Targeted therapy to driver gene mutations has improved the clinical outcome of some patients. Over recent years, immunotherapy has become a new therapeutic method for patients without targeting mutated genes. Despite the development of such treatments, the 5-year survival rate for patients with lung adenocarcinoma remains unsatisfactory. Most patients are in the median or advanced stage of disease when diagnosed and therefore the therapeutic effect is not as effective. Therefore, there is a significant need to diagnose patients during the early stages of disease as this would offer a better choice of treatment and is normally associated with a better outcome. In order to detect lung cancer earlier, it is necessary to investigate genes that are related to lung cancer and identify molecular biomarkers that can facilitate the routine diagnosis and treatment of lung cancer.

Research has shown that the dysregulation of genes and non-coding RNA molecules is associated with the development of many types of cancer, of which circRNAs and miRNAs are promising tumor biomarkers. miRNAs are non-coding RNAs of small molecules and can regulate the expression of genes and participate in some biological processes, such as transcription and endonuclease processes in the nucleus or cytoplasm. Studies have shown that the abnormal expression of miRNAs is closely related to various types of cancers and plays a key role in cell proliferation, apoptosis, and metastasis. The regulation of gene and miRNA expression can explain many cancer-related molecular mechanisms, such as those involved in colorectal cancer, lymphatic cancer, and lung cancer.[29,30] Some studies have predicted the survival prognosis of patients by constructing the expression profiles of miRNAs.

In order to identify miRNAs, and their associated genes, in lung adenocarcinoma, and to investigate their roles in the pathogenesis of lung adenocarcinoma, we analyzed differences in expression between lung cancer tissues and adjacent normal tissues by analyzing data from the TCGA database. We identified 6064 genes that were differentially expressed, including 5324 upregulated genes and 740 downregulated genes. We also identified 161 differentially expressed miRNAs, including 126 upregulated miRNAs and 35 downloaded miRNAs. Next, we mapped the co-expression network of differentially expressed genes and differential miRNAs with the WGCNA package and related software (Figure 3). We identified a number of differentially expressed genes that interacted with each other in the network, including VTN, ALB, PAQR9, AFP, CREB3L3, PRODH2, AC080128.1, HRG, APOB, PAQR9-AS1, LGALS14, and G6PC. Only up-regulated levels of hsa-mir-122 were correlated with these particular genes (Figure 3A). Previous studies have shown that hsa-mir-122 is associated with the development of different types of cancer, including breast cancer, renal cancer, colorectal cancer,[10,11] gastric cancer. The overexpression of hsa-mir-122 in gastric cancer patients has also been shown to inhibit the proliferation, metastasis, and invasion, of gastric cancer cells; its down-regulation promoted the proliferation of tumor cells. However, the overexpression of hsa-mir-122 in patients with colon cancer has been shown to promote the development, tumor metastasis and invasion.

In order to study the role of hsa-mir-122 and its associated genes in the pathogenesis of lung adenocarcinoma further, we chose hsa-mir-122 and genes that were associated with this miRNA and constructed a co-expression network using differentially expressed genes and miRNAs between lung adenocarcinoma and adjacent normal tissues (Figure 3A). We used the DAVID website[13,14] to perform GO (gene ontology) enrichment analyses of genes within the co-expression network for hsa-mir-122 and its associated differentially expressed genes (Figure 3A). Of the differentially expressed genes shown in Figure 3A, APOB, ALB, and VTN genes were found to be significantly enriched in terms of a specific biological process: receptor-mediated endocytosis (GO:0006898). This is a specific type of molecular transport that is mediated by cell surface receptors. Specific receptors on the cell surface can bind closely to extra-cellular macromolecules, such as ligands. This binding region contains a receptor-ligand complex for endocytosis; this creates a vesicle containing a receptor-ligand complex. This process usually occurs via clathrin pits and vesicles. ALB, HRG, and VTN genes were also significantly enriched in one particular cell component: blood cell microparticles (GO:0072562). This refers to particles in the blood that are derived from other types of cells, such as platelets, endothelial cells, or other types of cells. These particles express membrane receptors and other specific proteins from parental cells. These particles are also uneven in size and are devoid of nucleic acids. Studies have proposed that blood cell particles play a significant role in the pathophysiological mechanisms underlying inflammation and tumor metastasis. APOB, HRG, and VTN genes were significantly enriched in 1 specific molecular functional: binding heparin (GO:0008201). This refers to the selective and non-covalent binding to heparin. Heparin is a glycosaminoglycan and a key component of mast cells. Glycosaminoglycan is associated with malignant transformation in cells and tumor metastasis (Figure 4A). The results of GO functional analysis for differentially expressed genes related to hsa-mir-122 show that they are related to the development and metastasis of tumors.

We also used the KOBAS website (kobas.cbi.pku. edu.cn/) to perform Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. This technique showed that the pathway in which CREB3L3, G6PC and VTN genes are most significantly enriched was the PI3K-Akt signaling pathway (Figure 4B). The PI3K-Akt signaling pathway is one of the complex regulatory pathways found in human malignant tumors; this pathway is involved in the regulation of tumor development, including cell proliferation, genomic instability and metabolism. This pathway can affect the migration of endothelial cells by regulating recombination of the actin cytoskeleton. Because of the importance of this mechanism in tumors, this particular pathway has become a promising target for the therapy of cancer. A number of drugs that target this pathway are now involved in clinical trials. Many factors affect the production of vascular endothelial growth factor (VEGF) through the PI3K-Akt signaling pathway to promote angiogenesis. This pathway plays an important role in tumor development.[37-40] Therefore, CREB3L3, G6PC, and VTN genes play an important role in tumor angiogenesis and metastasis via the PI3K-Akt signaling pathway.

In order to investigate the role of hsa-mir-122 in the pathogenesis of lung adenocarcinoma further, we used TargetScan, miRDB,[18,19] PicTar[20,21] and miRanda[22,23] to predict the target genes of hsa-mir-122. These analyses identified27 target genes that were associated with hsa-mir-122, including ALDOA, CCNG1, GIT1, CS, DDR2, G3BP2, IQGAP1, STX16, DICER1, FOXP2, P4HA1, PAK3, MIPOL1, NPEPPS, OCLN, BACH2, DR1, SMYD4, EPO, MOSPD1, BRPF1, CPEB1, MASP1, MAF1, LRP10, LAMC1, and FUNDC2 genes. Functional analysis of these predicted target genes showed that they were significantly enriched in the following GO annotated entries: biological process (negative regulation of transcription by RNA polymerase II; GO:0000122; cellular component (nucleus; GO:0005634); molecular function (protein domain specific binding; GO:0019904). Of the predicted target genes, GIT1, IQGAP1, and PAK3, were significantly enriched in the actin cytoskeleton pathway. Cytoskeleton proteins are the driving force of the movement of tumor cells and also the core component of cellular pseudopodia. Their structural characteristics are the key in determining the migration ability of tumor cells.[41,42] Different types of pseudopodia play different roles in cellular metastasis and invasion. Of these, the lamellar pseudopodia play an important role in cellular migration. Invasive pseudopodia can help tumor cells degrade the extracellular matrix in order to access blood vessels and complete cellular invasion.[44,45] The production and increased number of pseudopodia in tumor cells are closely related closely to the metastasis and invasion of tumor cells. We hypothesize that hsa-mir-122 may interact with GIT1, IQGAP1, and PAK3 in the nucleus and bind to specific domains of RNA polymerase II sin a selective and non-covalent manner to inhibit RNA polymerase II-mediated transcription and to regulate cytoskeleton-related pathways. Therefore, hsa-mir-122 plays an important role in the development and metastasis of lung adenocarcinoma. Consequently, has-mir-122, and its related genes, should be investigated with regards to its role in regulating the cytoskeleton.

In order to investigate the effects of hsa-mir-122 and its related genes on the prognosis of patients with lung adenocarcinoma, we acquired 486 samples from patients with lung adenocarcinoma patients from the TCGA database, including patient survival time and status. By combining the expression of hsa-mir-122 and its related genes with patient survival data, we were able to analyze the prognosis of patients with lung adenocarcinoma using the survival package in R3.6.2 software. We found that the expression levels of G6PC, APOB, F2, PAQR9, and PAQR9-AS1 were upregulated in lung adenocarcinoma and also related to the prognosis of patients. The survival rate of patients with lung adenocarcinoma who had lower expression levels of these was higher (Figure 5). hsa-mir-122 was also highly expressed in lung adenocarcinoma. However, there was no significant correlation between the expression of this miRNA and prognosis. These findings suggest that hsa-mir-122 and its associated genes may affect the prognosis of patients by regulating the cytoskeleton, angiogenesis, and the metastasis of tumors.

To conclude, this study identified a specific miRNA molecule (hsa-mir-122) that is associated with lung adenocarcinoma and several genes (VTN, ALB, PAQR9, AFP, CREB3L3, PRODH2, AC080128.1, HRG, APOB, PAQR9-AS1, LGALS14, G6PC genes). Of these, CREB3L3, G6PC, and VTN are related to the PI3K-Akt signaling pathway. We also found that G6PC, APOB, F2, PAQR9, and PAQR9-AS1 were related to patient prognosis. Furthermore, we found that GIT1, IQGAP1, and PAK3 were associated with the actin cytoskeleton pathway. This suggests that hsa-mir-122 and its related genes may affect the prognosis of patients with lung adenocarcinoma by regulating the cytoskeleton, tumor angiogenesis, and tumor cell metastasis. High expression levels of of G6PC, APOB, F2, PAQR9, and PAQR9-AS1, are associated with low survival rates in patients with lung adenocarcinoma although there was no correlation between the expression of hsa-mir-122 and prognosis. This paper provides useful information relating to hsa-mir-122 and its related genes in the occurrence and development of lung adenocarcinoma and proposes new ideas for the targeted therapy of lung cancer, and potentially other forms of tumors.