Identification of potential therapeutic targets for lung cancer by bioinformatics analysis

Introduction

Lung cancer is one of the most common malignancies and has a significant socioeconomic impact on patients and their families (1). In western countries, the mortality rate of lung cancer is 15% and the worldwide mortality rate for patients with lung cancer is 86% (2). The high mortality of lung cancer is mainly attributable to the lack of effective therapeutic methods and the difficulty of obtaining an early diagnosis. Thus, the development of effective therapeutic targets is urgently required.

Differentially expressed genes (DEGs) have been reported to have important roles in lung cancer, and their identification may aid in the elucidation of its underlying molecular mechanisms as well as the discovery of novel biomarkers and treatments (3). Numerous genes, including p53 (3,4), EGFR (5,6), kRAS (7), PIK3CA (8) and EML4 (9), are known to be associated with lung cancer, while others have remained elusive. Futhermore, SEMA5A and -6A were identified as potential therapeutic targets for lung cancer (10–12). Although tremendous efforts have been made to discover novel targets for lung cancer treatments, the current knowledge is insufficient and requires expansion.

In the present study, DEGs between lung cancer and normal lung tissues were identified. Protein-protein interaction (PPI) and transcription factor (TF) regulatory networks were constructed and key target genes were screened. Through the identification of key genes, the possible underlying molecular mechanisms as well as potential candidate biomarkers and treatment targets for lung cancer were explored.

Materials and methods

Affymetrix microarray data

The gene expression profile dataset no. GSE3268 deposited in the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) by Wachi et al (13) based on the GPL96 platform (HG-U133A; Affymetrix Human Genome U133A Array), was subjected to bioinformatics analysis in the present study. The dataset contained a total of 10 chips, including five squamous cell lung cancer tissues and five paired adjacent normal lung tissues obtained from patients with squamous cell lung cancer.

Furthermore, the gene expression profile dataset GSE19804 based on the platform GPL570 (HG-U133_Plus_2; Affymetrix Human Genome U133 Plus 2.0 Array), which was deposited in the GEO database by Lu et al (14), was used. The dataset contained 120 chips, including 60 samples of non-small cell lung cancer tissues and 60 samples of paired normal lung tissues from female Taiwanese patients.

Identification of DEGs

The raw data were pre-processed using the Affy package (15) in R language. DEGs of GSE3268 (DEG1) and GSE19804 (DEG2) between normal groups and disease groups were respectively analyzed using the limma package in R (16). Fold changes (FCs) in the expression of individual genes were calculated and DEGs with P<0.05 and |log FC| >1 were considered to be significant. DEG1 and DEG2 were then combined and the pooled dataset was referred to as the overlapping DEGs in the present study.

Gene ontology (GO) and pathway enrichment analysis of DEGs

GO analysis is a commonly used approach for functional studies of large-scale transcriptomic data (17). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (18) contains information on networks of molecules or genes. The Database for Annotation, Visualization and Integrated Discovery (DAVID) (19) was used to systematically extract biological information from the large number of genes. GO functions and KEGG pathways of the overlapping DEGs were analyzed using DAVID 6.7 with P<0.05.

Construction of PPI network and screening of modules

The Search Tool for the Retrieval of Interacting Genes (STRING) (20) database was used to retrieve the predicted interactions for the DEGs; version 9.1 of STRING covers 1,133 completely sequenced species. All associations obtained in STRING are provided with a confidence score, which represents a rough estimate of the likelihood of a given association to describe a functional linkage between two proteins (21). The overlapping DEGs with a confidence score >0.4 were selected to construct the PPI network using Cytoscape software (version 3.0; http://cytoscape.org/) (22). Cytoscape allows for the visualization of complex networks and their integration to any type of attribute data. The MCODE (23) plugin in Cytoscape was used to divide the PPI into modules. GO functional analysis of genes in the modules was performed using the BinGo 2.44 plugin in Cytoscape (24) with a threshold of P<0.05 using the hypergeometric test.

Transcriptional regulatory network construction

The University of California at Santa Cruz (UCSC) database (http://genome.ucsc.edu) contains information on TF binding sites and the regulated genes (25). Using information collected from the UCSC database, DEGs were matched with their associated TFs. The TF regulatory network then was constructed using Cytoscape software (26).

Results

GO and pathway enrichment analysis of DEGs

From the GEO datasets, information on the expression of 8,172 genes was obtained. The normalized results showed that the expression median after normalization was in a straight line (Fig. 1). A total of 466 DEGs, including 156 upregulated and 310 downregulated genes, were selected.

Figure 1 Boxplot of normalized expression values for the datasets. The dotted lines in the middle of each box represent the median of each sample, and its distribution among samples indicates the level of normalization of the data, with a nearly straight line indicating a fair normalization level. Gene expression omnibus datasets: 1, GSE3268; 2, GSE19804.

Results of GO analysis showed that the upregulated DEGs were significantly enriched in biological processes, including collagen metabolic processes, multicellular organismal macromolecule metabolic processes and nuclear division (Table I); the downregulated DEGs were significantly enriched in biological processes, including response to wounding, immune response, defense response and inflammatory response (Table I).

Table I GO and pathway analysis of the differentially expressed genes.

Table I

GO and pathway analysis of the differentially expressed genes.

Expression	Category	Term/gene and function	Count	P-value
Upregulated	KEGG_PATHWAY	hsa04110 - Cell cycle	12	6.94×10−7
	KEGG_PATHWAY	hsa04512 - ECM-receptor interaction	10	1.50×10−6
	KEGG_PATHWAY	hsa04510 - Focal adhesion	10	1.42×10−3
	KEGG_PATHWAY	hsa04115 - p53 signaling pathway	6	2.14×10−3
	KEGG_PATHWAY	hsa00240 - Pyrimidine metabolism	5	3.93×10−2
	GOTERM_BP_FAT	GO:0032963 - Collagen metabolic process	9	2.10×10−10
	GOTERM_BP_FAT	GO:0044259 - Multicellular organismal macromolecule metabolic process	9	5.19×10−10
	GOTERM_BP_FAT	GO:0000280 - Nuclear division	17	5.79×10−10
	GOTERM_BP_FAT	GO:0007067 - Mitosis	17	5.79×10−10
	GOTERM_BP_FAT	GO:0000278 - Mitotic cell cycle	21	7.04×10−10
	GOTERM_BP_FAT	GO:0000087 - M phase of mitotic cell cycle	17	7.55×10−10
	GOTERM_CC_FAT	GO:0005576 - Extracellular region	53	1.41×10−10
	GOTERM_CC_FAT	GO:0005578 - Proteinaceous extracellular matrix	19	7.80×10−9
	GOTERM_CC_FAT	GO:0031012 - Extracellular matrix	19	2.50×10−8
	GOTERM_CC_FAT	GO:0044421 - Extracellular region part	30	2.27×10−7
	GOTERM_CC_FAT	GO:0005819 - Spindle	12	4.55×10−7
	GOTERM_MF_FAT	GO:0004222 - Metalloendopeptidase activity	9	9.37×10−6
	GOTERM_MF_FAT	GO:0048407 - Platelet-derived growth factor binding	4	1.53×10−4
	GOTERM_MF_FAT	GO:0004175 - Endopeptidase activity	13	3.80×10−4
	GOTERM_MF_FAT	GO:0004857 - Enzyme inhibitor activity	11	3.81×10−4
Downregulated	KEGG_PATHWAY	hsa04060 - Cytokine-cytokine receptor interaction	20	6.99×10−5
	KEGG_PATHWAY	hsa04610 - Complement and coagulation cascades	8	2.47×10−3
	KEGG_PATHWAY	hsa04062 - Chemokine signaling pathway	13	4.53×10−3
	KEGG_PATHWAY	hsa04650 - Natural killer cell mediated cytotoxicity	10	9.69×10−3
	KEGG_PATHWAY	hsa04614 - Renin-angiotensin system	4	1.01×10−2
	GOTERM_BP_FAT	GO:0009611 - Response to wounding	48	2.23×10−17
	GOTERM_BP_FAT	GO:0006952 - Defense response	46	1.66×10−13
	GOTERM_BP_FAT	GO:0006954 - Inflammatory response	33	2.92×10−13
	GOTERM_BP_FAT	GO:0006955 - Immune response	43	4.20×10−10
	GOTERM_BP_FAT	GO:0048545 - Response to steroid hormone stimulus	21	3.81×10−9
	GOTERM_CC_FAT	GO:0005615 - Extracellular space	55	2.36×10−18
	GOTERM_CC_FAT	GO:0044421 - Extracellular region part	64	2.03×10−17
	GOTERM_CC_FAT	GO:0005576 - Extracellular region	93	3.37×10−15
	GOTERM_CC_FAT	GO:0005886 - Plasma membrane	131	2.25×10−12
	GOTERM_CC_FAT	GO:0005887 - Integral to plasma membrane	61	1.99×10−11
	GOTERM_MF_FAT	GO:0019838 - Growth factor binding	16	2.01×10−9
	GOTERM_MF_FAT	GO:0030246 - Carbohydrate binding	27	7.86×10−9
	GOTERM_MF_FAT	GO:0019955 - Cytokine binding	13	1.54×10−6
	GOTERM_MF_FAT	GO:0005509 - Calcium ion binding	39	1.04×10−5
	GOTERM_MF_FAT	GO:0030247 - Polysaccharide binding	14	1.11×10−5

[i] BP, biological process; CC, cellular component; MF, molecular function; Count, numbers of differentially expressed genes; ECM, extracellular matrix; GO, gene ontology; hsa, Homo sapiens; KEGG, Kyoto Encyclopedia of Genes and Genomes; FAT, functional annotation tool.

Pathway analysis showed that the upregulated DEGs were significantly enriched in cell cycle, extracellular matrix - receptor interaction and the p53 signaling pathway (Table I); the downregulated DEGs were significantly enriched in cytokine receptor interaction, complement and coagulation cascades as well as chemokine signaling pathways (Table I).

Construction of PPI network and screening of module

The PPI network was constructed based on the predicted interactions of the identified DEGs (Fig. 2). Genes of IL-6, FOSB, CDK1, MMP9 and ICAM1 were found to have a high degree of interaction in lung cancer. A sub-network containing 34 nodes and 547 edges was screened from the PPI network, such as PTTG1 (Fig. 3). The DEGs in the sub-net were significantly enriched in biological processes, such as the cell cycle, and pathway analysis showed that they were significantly enriched in cell cycle and oocyte meiosis (Table II).

Figure 2 Protein-protein interaction network of the DEGs. Blue nodes represent products of upregulated DEGs and pink nodes represent products of down-regulated DEGs. The size of each node is proportional to the degree of nodes. DEG, differentially expressed gene.

Figure 3 Sub-network screened from protein-protein interaction network. Nodes refer to the products of upregulated differentially expressed genes.

Table II GO and pathway analysis of genes in sub-network.

Table II

GO and pathway analysis of genes in sub-network.

Category	Term/gene and function	Count	P-value
KEGG_PATHWAY	hsa04110 - Cell cycle	10	1.09×10−11
KEGG_PATHWAY	hsa04114- Oocyte meiosis	6	1.09×10−5
KEGG_PATHWAY	hsa04914 - Progesterone-mediated oocyte maturation	4	1.83×10−3
KEGG_PATHWAY	hsa04115 - p53 signaling pathway	3	1.65×10−3
KEGG_PATHWAY	hsa00240 - Pyrimidine metabolism	3	3.10×10−2
GOTERM_BP_FAT	GO:0000278 - Mitotic cell cycle	19	7.13×10−21
GOTERM_BP_FAT	GO:0007049 - Cell cycle	22	1.65×10−19
GOTERM_BP_FAT	GO:0000280 - Nuclear division	16	2.14×10−19
GOTERM_BP_FAT	GO:0007067 - Mitosis	16	2.14×10−19
GOTERM_BP_FAT	GO:0000087 - M phase of mitotic cell cycle	16	2.82×10−19
GOTERM_CC_FAT	GO:0005819 - Spindle	12	9.20×10−15
GOTERM_CC_FAT	GO:0000777 - Condensed chromosome kinetochore	8	3.94×10−11
GOTERM_CC_FAT	GO:0015630 - Microtubule cytoskeleton	14	5.31×10−11
GOTERM_CC_FAT	GO:0000779 - Condensed chromosome, centromeric region	8	1.01×10−10
GOTERM_CC_FAT	GO:0000922 - Spindle pole	7	1.01×10−10
GOTERM_MF_FAT	GO:0005524 - Adenosine triphosphate binding	15	4.89×10−7
GOTERM_MF_FAT	GO:0032559 - Adenyl ribonucleotide binding	15	5.78×10−7
GOTERM_MF_FAT	GO:0030554 - Adenyl nucleotide binding	15	1.10×10−6
GOTERM_MF_FAT	GO:0001883 - Purine nucleoside binding	15	1.32×10−6
GOTERM_MF_FAT	GO:0001882 - Nucleoside binding	15	1.44×10−6

[i] BP, biological process; CC, cellular component; MF, molecular function; Count, numbers of DEGs; GO, gene ontology; hsa, Homo sapiens; KEGG, Kyoto Encyclopedia of Genes and Genomes; FAT, functional annotation tool.

TF-target gene regulatory network analysis

Associations between 44 TFs and their 47 target DEGs were collected from the TF regulatory network (Fig. 4). TFs of FOSB and LMO2, which exhibited a high degree of interaction, were selected from this network. Furthermore, the results also showed that MMP9 was the target of FOSB and PTTG1 was the target of LMO2.

Figure 4 Transcriptional regulatory network analysis. Blue nodes represent products of upregulated DEGs and pink nodes represent products of downregu-lated DEGs. Triangle arrowheads indicate transcription factors and circles indicate target genes. DEG, differentially expressed gene.

Discussion

Lung cancer is the leading cause of cancer-associated mortality; however, the underlying molecular mechanisms of its development and progression have remained to be fully elucidated (1). The present study used a bioinformatics approach to predict the potential therapeutic targets and explore the possible molecular mechanisms for lung cancer. A total of 466 DEGs between tumorous and normal tissues was identified, among which 310 genes were downregulated and 156 were upregulated. By constructing a PPI network and a TF regulatory network, key genes, including IL6, MMP9 and PTTG1, were identified.

IL-6 is a multifunctional cytokine that was characterized as a regulator of immune and inflammatory responses (27,28). It is involved in the regulation of cell proliferation, survival and metabolism, and IL-6 signaling has an important role in tumorigenesis (29). Chung et al (30) found that IL-6 activated PI3K, which promoted apoptosis in human prostate cancer cell lines. Furthermore, studies have shown that IL-6 inhibited the growth of numerous types of cancer, including lung (31), breast (32) and prostate cancer (33). In the present study, IL-6 was shown to be downregulated in squamous cell and non-small cell lung cancer, and GO analysis showed that IL-6 was significantly enriched in biological processes, including defense response, inflammatory response, immune response and regulation of cell proliferation, which was consistent with a previous study (29). Combined with the above studies, it is indicated that IL-6 may be a diagnostic biomarker and therapeutic target in lung cancer.

MMP9 has a key role in cell migration, proliferation, differentiation, angiogenesis, apoptosis and host defense (34). Dysregulatoin of MMPs has been implicated in numerous diseases, including chronic ulcers and cancer (35–37). Downregulation of MMPs has been shown to inhibit metastasis, while upregulation of MMPs led to enhanced cancer cell invasion (37). In the present study, MMP9 was overexpressed and regulated by FOSB in lung cancer tissues. Kim et al (38) found that FOSB was downregulated in pancreatic cancer and promoted tumor progression. Kataoka et al (39) found that FOSB gene expression in cancer stroma is a independent prognostic indicator for patients with epithelial ovarian cancer receiving standard therapy. Combined with the above studies, the present study indicated that MMP9 may have important roles in the progression of lung cancer, and that it may be utilized as a therapeutic target.

PTTG1 has tumorigenic activity and is highly expressed in various tumor types (40). Studies have shown that PTTG1 was overexpressed in esophageal cancer and associated with endocrine therapy resistance in breast cancer (41,42). Yoon et al (40) showed that the PTTG1 oncogene promoted tumor malignancy via epithelial-to-mesenchymal expansion of the cancer stem cell population. Hamid et al (43) found that PTTG1 promoted tumorigenesis in human embryonic kidney cells. A study by Li et al (44) indicated that PTTG1 promoted migration and invasion of human non-small cell lung cancer cells. Panguluri et al (45) showed that PTTG1 was an important target gene for ovarian cancer therapy. In the present study, PTTG1 was found to be overexpressed in lung cancer tissues and regulated by LMO2. LMO2 is an important regulator in determining cell fate and controlling cell growth and differentiation (46). Nakata et al (47) found that LMO2 was a novel predictive biomarker with the potential to enhance the accuracy of prognoses for pancreatic cancer. Yamada et al (48) showed that LMO2 is a key regulator of tumour angiogenesis. Combined with the above studies, the present study indicated that PTTG1 may have important roles in the progression of lung cancer and that it may represent a therapeutic target.

In conclusion, the bioinformatics analysis of the present study indicated that IL-6, MMP9 and PTTG1 may have key roles in the progression and development of lung cancer. They may be used as prognostic biomarkers as well as specific therapeutic targets for the treatment of lung cancer. However, molecular biology experiments are required to confirm these findings.

References