Transcriptomic and functional network features of lung squamous cell carcinoma through integrative analysis of GEO and TCGA data

Lung squamous cell carcinoma (LUSC) is associated with poor clinical prognosis and lacks available targeted therapy. Novel molecules are urgently required for the diagnosis and prognosis of LUSC. Here, we conducted our data mining analysis for LUSC by integrating the differentially expressed genes acquired from Gene Expression Omnibus (GEO) database by comparing tumor tissues versus normal tissues (GSE8569, GSE21933, GSE33479, GSE33532, GSE40275, GSE62113, GSE74706) into The Cancer Genome Atlas (TCGA) database which includes 502 tumors and 49 adjacent non-tumor lung tissues. We identified intersections of 129 genes (91 up-regulated and 38 down-regulated) between GEO data and TCGA data. Based on these genes, we conducted our downstream analysis including functional enrichment analysis, protein-protein interaction, competing endogenous RNA (ceRNA) network and survival analysis. This study may provide more insight into the transcriptomic and functional features of LUSC through integrative analysis of GEO and TCGA data and suggests therapeutic targets and biomarkers for LUSC.


GSE8569
Journal downloaded from the repository. Principal component analysis (PCA) was done for the datasets for dimensionality reduction and quality control. If the quality of a particular sample is not good enough, it would be excluded for subsequent analysis. Details of each microarray study, including sample descriptions are provided in Table 1.
Our workflow for bioinformatics analysis of publicly available datasets from both GEO and TCGA databases is illustrated in Fig. 1.

Differential expression analysis.
Heterogeneity and potential variables are commonly recognized as major sources of bias and variability in high-throughput experiments. Since the datasets we recruited for our multi-datasets analysis were based on different platforms and samples were handled on different days, in different groups or by different people. Besides, datasets GSE40275 and GSE61223 only have 5 and 2 tumor samples respectively and using few samples can affect the performance of statistical analysis and provides unreliable results. Therefore, we first integrated all samples of seven datasets to significantly improve the number of samples (61 normal samples vs. 88 tumor samples) so as to avoid generating less reliable results followed by batch normalization in the R computing environment using sva package 9 . The unnormalized raw data was summarized as the form of the matrix and can be acquired in Supplementary Table 1. Next, we performed the differential analysis (|Log 2 FC| > 2, adjusted p-value < 0.05) by comparing tumor tissues to normal tissues in the R computing environment using limma package 10 . Integration of the differentially expressed genes in TCGA database. The Cancer Genome Atlas (TCGA), a project supported by the National Cancer Institute (NCI) and National Human Genome Research Institute (NHGRI), has generated comprehensive, multi-dimensional maps of the key genomic changes in various types of cancers. In order to obtain a consensus of differentially expressed genes, gene expression quantification data and clinical information of LUSC patients in TCGA database were downloaded using TCGAbiolinks 11 . All data were normalized and processed with TCGAbiolinks pipeline. The TCGAbiolinks principle of differential analysis is to first convert the count matrix into an edgeR object 12 , then each gene gets assigned the same dispersion estimate, then performs pair-wise tests for differential expression between two groups, and finally takes the output using the False Discovery Rate (FDR) correction, and returns the top differentially expressed genes 11 . The parameters set for differential expression analysis were FDR < 0.05 with |Log 2 FC| > 2. Subsequently, we combined the differentially expressed genes acquired from GEO and TCGA databases to get the convergence gene signatures.  Construction of ceRNA network. To find out whether these 129 genes exist competing endogenous regulating network mediated by long non-coding RNAs (lncRNAs) and micro RNAs (miRNAs). A competing endogenous RNA (ceRNA) network was built using GDCRNATools 16 . The major criteria of building ceRNA network in GDCRNATools are: (1) The lncRNA and mRNA must share a significant number of miRNAs. (2) Expression of lncRNA and mRNA should be positively related. (3) miRNAs should play similar roles in regulating the expression of lncRNA and mRNA. We followed the pipeline of GDCRNATools to first identify differentially expressed lncRNAs (DElncRNAs) and differentially expressed miRNAs (DEmiRNAs) using the edgeR 12 method (FDR < 0.05 with |Log 2 FC| > 2). Next, we used the function of GDCRNATools to construct the network, total read counts for 5p and 3p strands of DEmiRNAs were obtained from isoform quantification files, miRcode was used to collect predicted and experimentally validated lncRNA targets 17 . StarBase v2.0 was used to predict miRNA-mRNA interactions 18 . Visualization of the ceRNA was performed by Cytoscape 19 .  Survival analysis. To see whether these 129 genes and DElncRNAs were related to prognostic significance, survival analysis was performed in the R environment using TCGAbiolinks 11 . We used clinical information to plot the survival curves for 1/3 of patients with higher expression of a specific gene versus the 1/3 of patients with lower expression of this gene (p < 0.05).

Principal component analysis verifying independence of each group.
To distinguish the significant difference between normal and tumor samples of GEO data, PCA was performed to reduce the dimensionality and evaluate the independence of each group. The results showed that normal samples vs. tumor samples in the datasets (GSE8569, GSE21933, GSE33532, GSE40275, GSE62113, GSE74706) displayed a significant difference except for dataset GSE33479, whose two tumor samples GSM828337 and GSM828345 were close to normal samples, so we removed these two samples for the subsequent analysis (Fig. 2B). The contribution rate for each principal component is on the vertical axis ( Fig. 2A). The cumulative contribution rates of the PC1 and PC2 of each of the seven datasets are 27.64%, 39.50%, 28.94%, 65.74%, 61.05%, 57.85% and 45.44% respectively. The horizontal axis stands for the number of principal components required to reach a cumulative proportion of 100%. It was obvious that the first two components were enough to separate the two groups, indicating each group is independent of each other (Fig. 2B).

Convergence of gene expression signatures across different studies of LUSC. Data integration
is becoming increasingly necessary to investigate the complex genetic mechanisms by applying appropriate statistical method 20 . As the outputs of individual experiments can be rather noisy, it is essential to look for findings that are supported by several pieces of evidence to increase the signal and lessen the fraction of false positive findings. We used batch correction to reduce variability and then used limma package 10 (|Log 2 FC| > 2, adjusted P value < 0.05) to identify differentially expressed genes. Table 1 demonstrates the number of differentially expressed genes identified from each of the seven datasets and TCGA data. Volcano plots in Fig. 3A showed the number of differentially expressed genes identified from each of the seven datasets and the number of differentially expressed genes after batch correction. We found 94 up-regulated genes and 39 down-regulated genes after batch normalization (Fig. 3A). For TCGA data, we found a total of 2242 differentially expressed genes with 1477 of them up-regulated and 765 genes down-regulated. Here, we demonstrate the names of genes with |Log 2 FC| > 8 (Fig. 3B). Venn diagram demonstrates the intersections of genes between GEO data and TCGA data, and 129 co-differentially expressed genes (91 up-regulated and 38 down-regulated) were found (Fig. 3C). Chromosome mapping of consensus genes revealed chromosome distribution, with chromosomes 1 containing the greatest number of dysregulated genes in LUSC (Fig. 3D). Interestingly, while two genes on the X chromosome showed dysregulation in LUSC (FHL1 and FIGF), not a single Y chromosome gene was affected.
In Fig. 4A,B we displayed the expression changes of these genes in GEO and TCGA data, respectively. More information including the fold change and FDR of these 129 genes can be found in Supplementary Table 2. These 129 genes were further subjected to functional annotation and protein to protein interaction analysis to determine the biological significance of this cross-study convergence in the pathogenesis of LUSC.
The complete results of GO and KEGG analyses can be found in Supplementary Table 3. Functional enrichment analysis indicated that expression changes of these genes could lead to increased activities of proliferation of cells, cell proliferation of tumor cell lines, invasion of cells, cell survival, migration of cells and cell movement in LUSC and decreased activities of organism death, cell movement of leukocytes, apoptosis of tumor cell lines, cell movement of blood cells, leukocyte migration, migration of blood cells and necrosis. All these functions are critically important in tumor cell survival, invasion and immune escape (Fig. 5C). Specific data of functional enrichment analysis can be found in Supplementary Table 4. Figure 5D showed the protein-protein interaction network. PPI enrichment analysis was done with the following databases: BioGrid16 21 , inWeb_IM17 22 and OmniPath18 23 . Molecular Complex Detection (MCODE) algorithm 24 was further applied to identify densely connected network components if there are more than two proteins in a network. We found that CCNB2, PLK1, KIF2C, CENPA, CENPF, BUB1, BUB1B, BIRC5, CENPE, ZWINT, AURKB, CHEK1, EXO1, RAD51, and RFC4 can interact with each other and this interaction was predominantly associated with protein serine/threonine kinase activity.
Survival analysis. Base on TCGA data and clinical information, we analyzed the survival curves for patients by comparing 1/3 of patients with higher expression of a certain gene to 1/3 of patients with lower expression. Of the 129 genes, we found that 60 genes were statistically related to the overall survival rate (p < 0.05). Here, we exhibited 20 examples of these genes (Fig. 7), the remaining can be found in Supplementary Figs 1 and 2. Expression changes of these 60 genes can be obtained in Supplementary Table 7. For these 60 genes, EZH2, ABCC5, and KIF23 were in the ceRNA network and could be modulated by corresponding lncRNAs and miR-NAs. EZH2, ABCC5, and KIF23 were up-regulated in LUSC and patients with low expression levels of these three genes had shorter survival times (Fig. 7, Supplementary Figs 1 and 2). We also found that LncRNAs KC6, PART1, SFTA1P, and SNHG1 were statistically related to the overall survival rate ( Supplementary Fig. 3, p < 0.05). Functional enrichment analysis indicated that the 60 overall survival related-genes were involved in the cell proliferation of tumor cell lines, perinatal death, invasion of cells, organism death, proliferation of cells, neonatal death and migration of cells ( Supplementary Fig. 4 and Table 8).

Discussion
LUSC has been regarded as the "neglected sibling" compared with lung adenocarcinoma due to lack of effective targeted treatment options. The mutations of epidermal growth factor receptor (EGFR) kinase, as well as fusions in the anaplastic lymphoma kinase (ALK), lead to a dramatic change in the treatment of patients with lung adenocarcinoma [25][26][27] . Unfortunately, EGFR mutations and ALK fusions are typically not present in LUSC 28 , and novel targeted agents for adenocarcinoma of the lung ineffective against LUSC. So, new classes of biomarkers with high efficiency, high specificity, and high sensitivity are required as novel molecules for diagnosis and prognosis of LUSC.
Integrating multiple individual data has been showed to improve detection power 29 . Integration of multiple arrays is considered a better approach of enhancing the reliability of results than individual array analysis. PCA is a sophisticated technique widely used for reducing the dimensions of multivariate problems and evaluating independence without losing much information 30 . In our present studies, PCA results showed that tumor groups were independent of normal groups in each of the seven datasets (GSE8569, GSE21933, GSE33479, GSE33532, GSE40275, GSE62113, GSE74706). We identified 129 (91 up-regulated and 38 down-regulated) intersections of genes between GEO data and TCGA data. Chromosome mapping of consensus genes showed chromosomes 1  containing the greatest number of dysregulated genes in LUSC. Previously studies confirmed that the skewed X chromosome inactivation was associated with early development of lung cancer in females. The X chromosomal inactivation assay may be used to screen for females predisposed to malignancies including lung cancer 31 . Our results indicated that the dysregulation of FHL1 and FIGF on X chromosome may be associated with LUSC in females. On the other hand, Mosaic loss of the Y chromosome has a striking association with aging and cigarette smoking 32 . In our present study, that no differentially expressed gene was found in Y chromosome may be related to loss of Y chromosome gene. We found that up-regulated genes were predominantly enriched in the activities of nuclear division, organelle fission, mitotic nuclear division, ATPase activity, microtubule binding and microtubule motor activity in LUSC. Meantime, down-regulated genes were enriched in humoral immune response, regulation of inflammatory response, regulation of cell growth, carboxylic acid binding, and response to transforming growth factor beta in LUSC. Previous studies showed that mitotic nuclear division is associated with cell proliferation, dysfunction of this process can lead to mitotic checkpoint failure and cause chromosome missegregation 33,34 . Microtubules function in the precise segregation of chromosomes during cell division, transport of cellular cargos, and positioning and movement of intracellular organelles 35 . Microtubule-binding drugs including the Vinca alkaloids and taxanes can suppress the dynamic instability of microtubules and induce apoptosis 36 . KEGG pathway enrichment analysis suggested significant enrichment in pathways including cell cycle and p53 signaling pathway. Our results indicated that the changes in biological processes, cellular components, molecular functions, and pathways may play critically important roles in the pathogenesis of LUSC. Protein-protein interaction network illustrated the overview of their functional connections. Module analysis of the PPI network suggested that protein serine/ threonine kinase activity might be involved in LUSC development. Above are critical cellular processes for maintenance of cell homeostasis, dysregulation of these processes tends to promote carcinogenesis 37,38 . Our findings highlighted the probable importance of the regulation of these key biological behaviors by aberrantly expression in LUSC which warranted further investigations to confirm.
Previous studies confirmed that Enhancer of zeste homolog 2 (EZH2), which is a histone methyltransferase, can regulate gene expression by catalyzing tri-methylation of histone H3 at Lys 27 (H3K27me3) 39 . Behrens, C. et al. found that over expression of EZH2 was associated with tumor progression in lung cancer 40 . However, interestingly, it has been reported that EZH2 can also act as a tumor suppressor gene 41 . In our study, EZH2 was over-expressed and its higher expression predicted longer survival time for LUSC patients, indicating its potential tumor suppressor role in LUSC. ABCC5 functions have been regarded as a mediator of breast cancer skeletal metastasis. ABCC5 may be a potential therapeutic target for breast cancer bone metastasis 42 . KIF23 (Kinesin family member 23) is an important regulator of cellular cytokinesis, and it has been considered a tumor gene is glioma 43 . But its relationship with LUSC is largely unknown at present. A growing number of studies have confirmed that the lncRNAs-miRNAs-mRNAs regulation network functions in tumor pathogenesis and progression 38,44,45 . In our present study, ceRNA network found that PTHLH, EZH2, CEP55, CCNA2, PFN2, ABCC5, ANLN, UCK2, DSG2, GREM1, MYBL2, PITX1, CHEK1, KIF23 could be modulated by lncRNAs through corresponding miRNAs. This regulation network could provide us more knowledge of the sophisticated regulation patterns in LUSC. Strikingly, we also identified that 60 genes were statistically related to the overall survival rate. These overall survival-related genes were involved in the invasion of cells, proliferation of cells, respiratory of system tumor, differentiation of cells, and apoptosis. Previous studies reported that PART1 was associated with poor prognosis and tumor recurrence in stage I-III non-small cell lung cancer 46 . SFTA1P were regarded as a tumor suppressor. Down-regulation of SFTA1P may be associated with decreased TP53 expression 47 . LncRNA SNHG1 promoted non-small cell lung cancer progression 48 . In our present study, we found that over-expression of KC6, PART1, and SNHG1 were associated with poor prognosis in LUSC. However, lower expression of SFTA1P was associated with poor prognosis in LUSC.
In summary, our study analyzed the array-based and sequence-based data of LUSC supported by GEO and TCGA databases. We discovered a number of co-differentially expressed genes and important pathways in LUSC. Based on these genes, we performed a series of analyses, which may contribute to the finding of molecular mechanisms underlying the initiation and development of LUSC.