Next Article in Journal
Extracellular Vesicles in Non-Small-Cell Lung Cancer: Functional Role and Involvement in Resistance to Targeted Treatment and Immunotherapy
Previous Article in Journal
Structure and Optimization of Checkpoint Inhibitors
 
 
Journals
Cancers
Volume 12
Issue 1
10.3390/cancers12010037
cancers-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA)

by
Magdalena Niemira
1,*,
Francois Collin
2,
Anna Szalkowska
1,
Agnieszka Bielska
1,
Karolina Chwialkowska
2,
Joanna Reszec
3,
Jacek Niklinski
4,
Miroslaw Kwasniewski
2 and
Adam Kretowski
1,5
1
Clinical Research Centre, Medical University of Bialystok, 15-276 Bialystok, Poland
2
Centre for Bioinformatics and Data Analysis, Medical University of Bialystok, 15-276 Bialystok, Poland
3
Department of Medical Pathomorphology, Medical University of Bialystok, 15-276 Bialystok, Poland
4
Department of Clinical Molecular Biology, Medical University of Bialystok, 15-276 Bialystok, Poland
5
Department of Endocrinology, Diabetology and Internal Medicine, Medical University of Bialystok, 15-276 Bialystok, Poland
*
Author to whom correspondence should be addressed.
Cancers 2020, 12(1), 37; https://doi.org/10.3390/cancers12010037
Submission received: 12 November 2019 / Revised: 16 December 2019 / Accepted: 19 December 2019 / Published: 21 December 2019
Browse Figures
Versions Notes

Abstract

:
Non-small-cell lung cancer (NSCLC) represents a heterogeneous group of malignancies consisting essentially of adenocarcinoma (ADC) and squamous cell carcinoma (SCC). Although the diagnosis and treatment of ADC and SCC have been greatly improved in recent decades, there is still an urgent need to identify accurate transcriptome profile associated with the histological subtypes of NSCLC. The present study aims to identify the key dysregulated pathways and genes involved in the development of lung ADC and SCC and to relate them with the clinical traits. The transcriptional changes between tumour and normal lung tissues were investigated by RNA-seq. Gene ontology (GO), canonical pathways analysis with the prediction of upstream regulators, and weighted gene co-expression network analysis (WGCNA) to identify co-expressed modules and hub genes were used to explore the biological functions of the identified dysregulated genes. It was indicated that specific gene signatures differed significantly between ADC and SCC related to the distinct pathways. Of identified modules, four and two modules were the most related to clinical features in ADC and SCC, respectively. CTLA4, MZB1, NIP7, and BUB1B in ADC, as well as GNG11 and CCNB2 in SCC, are novel top hub genes in modules associated with tumour size, SUVmax, and recurrence-free survival. Our research provides a more effective understanding of the importance of biological pathways and the relationships between major genes in NSCLC in the perspective of searching for new molecular targets.

1. Introduction

Lung cancer is the second commonly diagnosed cancer in both men and women, as well as the leading cause of cancer-related mortality in the world. According to the American Cancer Society, approximately 228,150 new cases of lung cancer were diagnosed, and over 142,000 deaths were registered in the United States in 2019 alone [1]. Lung cancer is classified as either small-cell lung carcinoma (SCLC) or non-small-cell lung carcinoma (NSCLC). NSCLC is diagnosed in four out of five cases of lung cancer. The main subtypes of NSCLC are squamous cell carcinoma (SCC), adenocarcinoma (ADC), and large-cell carcinoma (LCC) [2]. Unfortunately, despite the introduction of new treatment options in recent years, the survival rates of patients with NSCLC remain unsatisfactory. Currently, the gold standard is the therapy of NSCLC based on histopathological features. However, progression, metastasis susceptibility, and potential response to treatment cannot be fully predicted based on histopathological observation. Furthermore, even if NSCLC patients have similar clinical characteristics and have tumours at a similar stage, they may experience different clinical outcomes [3]. Recent studies displayed that the histologic types of lung cancer differ in their molecular features, such as mutational signatures, DNA copy number, and gene expression [4,5,6,7]. Due to this fact, the molecular genetic classification of cancers is attracting widespread interest nowadays. The next-generation sequencing (NGS) plays a key role in cancer profiling. RNA-seq avoids the limitation of microarray analysis and allows enhanced detection of low-abundance transcripts. However, the attempt to explain the complex relationship between genes is still a challenge [8,9]. The differential expression analysis allows considering each gene individually and does not account for the connection between genes. Since it is known that biological processes are based on the numerous interactions between many cellular components, gene network analysis could provide precious information on cancer pathogenesis [10].
The Weighted Correlation Network Analysis (WGCNA) [11,12] corresponds to a data reduction method and unsupervised classification method. It simplifies the interpretation of thousands of gene responses to a dozen of synthetic groups (or modules) of genes. The net establishes connections between genes-genes are connected if their expression is correlated. Genes can be more or less intensively connected depending on the value of the correlation (the weights). The connectivity between genes is then interpreted into a distance, and the distance is used to group genes into modules. This is how a high number of genes can be reduced into a small number of clusters whose expression is quantified by the Eigengenes (first principal component within the module). It relies on the assumption that highly correlated genes within a module are involved in common biological processes. For example, Tian et al. [13] used WGCNA to identify two gene co-expression modules involved in the process of lung SCC metastasis and suggested that CFTR, SCTR, and FIGF genes could be used as a potential prognostic biomarker. Using WGCNA in the studies of progression and prognosis of breast cancer indicated the relevant correlation of one module of genes with patient survival time. Moreover, the eight candidate biomarkers of the pathogenesis of breast cancer have been identified [14]. Yin et al. [15] found through WGCNA that five high degree hub genes may play a key role in the hepatocellular carcinoma (HCC) progression. Additional studies using WGCNA have shown that SKA1, ERCC6L, and GTSE-1 may be potential diagnostic markers for renal cell carcinoma [16]. Furthermore, WGCNA has been also applied for RNA-seq and miRNA-seq data of colon adenocarcinoma (CAC). These studies defined two gene modules and one miRNA module correlated with the pathological stage of CAC, as well as the putative representative biomarkers related to prognosis were defined [17]. In addition, Nakamura et al. [18] showed in proteomics that data analysis using WGCNA was notably effective to group proteins into modules and to identify of hub proteins, which played key roles in small-cell lung carcinoma and large-cell neuroendocrine lung carcinoma.
The broad goal of the current study was to identify transcriptome profile associated with the histological subtypes of NSCLC in a well-characterised group of patients from MOBIT project using RNA-seq. The Molecular Biomarkers for Individualized Therapy (MOBIT) [19] is an ambitious study targeting the development of personalised diagnostics based on high-quality biobanking and integrated analysis of ‘omics’ and imaging PET/MRI (Positron Emission Tomography-Magnetic Resonance Imaging) analysis, including the phenotypic, clinical, and epidemiological traits. As part of the project MOBIT, biospecimens of lung tumour and normal tissues were collected under strictly monitored conditions. MOBIT’s unique approach to the acquisition of tumour samples considers the following parameters: time of anaesthesia start, time of surgery start, time of resection start, time of vessel ligation, and time of resection end. The above parameters affect the molecular reality of individual tumours. The tissue ischemia time affects gene and protein expression patterns within minutes following surgical tumour excision [20]. The know-how of managing the oncology biobank was obtained from a scientific partner-Indivumed GmbH [21,22]. Differential expression analysis was performed to detect differentially expressed genes (DEGs) in ADC and SCC tissues vs. normal tissues. In addition, we applied the weighted gene co-expression network analysis (WGCNA) to find highly correlated gene modules with clinical traits and to identify the key genes in selected modules with the highest connectivity in the module in question.

2. Results

2.1. Baseline Characteristic

During the study program from 2015–2016, 114 patients with NSCLC were enrolled in the MOBIT project and underwent surgical tumour resection [19]. The mean age at diagnosis was 65.4 years, and 80 patients (70.2%) were male. One hundred and seven patients were smokers, with the average number of pack years was 39.6, whereas seven patients never smoked. Fourty-five (39.5%) had adenocarcinoma, and 69 (60.5%) had squamous cell carcinoma. Stage I disease was identified in 39 patients, stage II in 52 patients, stage III in 17 patients, and stage IV in six patients (Table 1). During lung tumour excision procedures, one to three fragments of lung tissues were collected from each patient and classified by pathologists as normal (Tn) and tumour (Tp, peripheral and/or Tc, central).

2.2. Gene Expression Profile Data in NSCLC Subtypes

A total of 228 samples, 114 tumour and 114 normal tissues, were included in the differential expression analysis with the count depth of higher than 15 million reads. Analysis was done in the software R along with two packages, edgeR [23] and DESeq2 [24]. Due to the large size of the tested genes, raw p values were adjusted according to Benjamini and Hochberg’s method. The selection criteria were strengthened with a threshold of FDR ≤ 0.05 and |log2FC| ≥ 1. Volcano plots were plotted to illustrate the distribution of each gene according to the log2FC and FDR (Figure 1A). These stringent criteria generated a list of 7720 genes differentially expressed in SCC in comparison to normal tissue, with 3233 up-regulated and 4487 down-regulated genes; 5856 of genes were significantly differentially expressed in ADC in comparison to normal tissue, with 2452 up-regulated and 3404 down-regulated genes (Figure 1B,C). Interestingly, among the up-regulated and down-regulated genes, 4757 (53.9%) were overlapping between SCC and ADC, while 4062 (46.1%) were specific to SCC or ADC.
The top 20 genes with the greatest upregulation and downregulation in cancer tissues are summarised in Figure 2. The common upregulated genes in SCC and ADC included MAGEA9, KRT20, HOXD11, CST4, MAGEA1, MAGEA10, MAGEA4, FGF19, and CALM5. These genes are mainly associated with extracellular matrix and proliferation, promotion of tumorigenesis, and apoptotic processes [25,26,27]. Interestingly, as many as four top up-regulated genes in lung cancer tissue were homologous genes from the MAGE-A family. It should be noted that the remaining genes in ADC from the top 20 were also up-regulated in SCC but at a lower level. Nevertheless, we found that CSBF, PAGE2, CALM3, DSG3, CTAG2, and FAM83C were up-regulated only in SCC.
The common downregulated genes in SCC and ADC included SFTPC, CLDN18, CHIAP2, SLC6A4, CD300LG, RS1, AGER, and MYOC. Majority of these genes are mainly involved in differentiation, cell death, cell cycle progression, regulation of immune response, migration, and lung adenocarcinoma formation [28,29,30]. The most down-regulated gene in both subtypes is SFTPC, encoding surfactant protein c, a hydrophobic protein with the capacity of mitigating surface tension in the lung tissue [31]. We noticed that other top down-regulated genes in ADC were also down-regulated in SCC but at a lower level.

2.3. Gene Ontology and Canonical Pathway Analysis

To identify biologic pathways, networks, and functional categories of differentially expressed genes, we used a combination of bioinformatic tools: the Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., https://qiagenbioinformatics.com/products/ingenuity-pathway-analysis), the gene ontology (GO) enrichment analysis-biological processes (BP) by BiNGO tool [32] as plugin for Cytoscape version 3.7.2 (http://cytoscape.org/), and a novel tool to provide a comprehensive gene list annotation based on over 40 independent knowledge bases, Metascape (http://metascape.org) [33].
Among all the GO-terms (BP), only 97 (52.2%) overlapped between SCC and ADC. GO-enrichment analysis revealed that among the DEGs were those involved in cell cycle regulation, proliferation, metabolic processes, and DNA repair (Figure 3A). These processes are usually disrupted in cancer and are potentially involved in the initiation and/or promotion of tumorigenesis. When we looked at the GO-terms in each subtype of NSCLC, we noticed that, in SCC, the deregulated processes (31.7%) were focused around cornification (FDR = 1.47 × 10−12), epidermis development (FDR = 9.11 × 10−10), keratinisation (1.49 × 10−9), and epidermal cell differentiation (3.2 × 10−6) (Figure 3B). In contrast, GO-terms in ADC (16.1%) were related with the humoral immune response (1.21 × 10−17), production of molecular mediator of immune response (1.83 × 10−14), activation of immune response (1.32 × 10−13), and positive regulation of immune response (2.98 × 10−10) (Figure 3C). The most significant pathways identified by Metascape in SCC included cornification, epidermal cell differentiation, chemotaxis, matrisome-associated, core matrisome, developmental growth, and ECM regulators (Figure 3D). Humoral immunity, activation of immune response, leukocyte migration, chemotaxis, cell-cell adhesion, matrisome-associated, and mitotic nuclear division were the most significant pathways in ADC tissues (Figure 3E).
Ingenuity core analysis identified 70 and 54 altered canonical pathways in SCC and ADC, respectively (FDR ≤ 0.01). The top 10 pathways of each subtype with a highly significant correlation, in comparison to the normal tissue, are listed in Figure 4A. Among them, 42 (51.2%) canonical pathways were common to both subtypes. Eicosanoid signalling, agranulocyte and granulocyte adhesion and diapedesis, atherosclerosis signalling, and axonal guidance signalling were the most significant changed pathways in both subtypes of NSCLC.
IPA analysis for all differentially expressed genes in SCC vs. normal tissue revealed basal cell carcinoma signalling (FDR = 3.72 × 10−6, ratio = 33/70), p38 MAPK signalling pathway (FDR = 1.29 × 10−5, ratio = 32/54), p53 signalling pathway (FDR = 7.59 × 10−3, ratio = 34/62), and noradrenaline degradation (FDR = 7.05 × 10−4, ratio=16/30) as the top canonical pathways, not disclosed in ADC (Figure 4B). Overall, these pathways and GO-terms results indicate that epithelial squamous cells transformation and excessive proliferation play an important role in SCC pathogenesis.
The top dysregulated canonical pathways detected in ADC were MIF regulation of innate immunity (FDR = 2.88 × 10−5, ratio = 16/42), eNOS signalling (FDR = 2.54 × 10−3, ratio = 16/42), phospholipases (FDR = 2.95 × 10−3, ratio 20/57), and Wnt/β-catenin formation (FDR = 3.46 × 10−3, ratio 47/170). Interestingly, most of the biological processes and pathways activated in ADC were related to tumour immune response (Figure 4C).

2.4. Upstream Transcription Regulators

The upstream regulator analysis tool is a function in IPA which may, using analysing linkage to DEGs through coordinated expression, identify assumed upstream regulators that have been observed experimentally to affect gene expression [34]. We identified the upstream transcriptional regulators with FDR ≤ 0.01 and the absolute activation z-score > 2 or <–2. Among activated or inhibited regulators, only 90 (36.8%) were overlapping between SCC and ADC (Figure 5A). As shown in Figure 5B, the majority of regulators were transcriptional regulators, growth factors, kinases, and enzymes, etc., which are related to cell growth, cell cycle, invasion, and apoptosis. ERBB2, with the highest activation z-score, encodes a member of the epidermal growth factor (EGF) receptor family of receptor tyrosine kinases, which is involved in the development of human malignancies [35,36]. CDKN1A, with the lowest z-score, encodes cyclin-dependent kinase inhibitor 1, which plays a major part in cell cycle arrest in response to DNA damage [37].
However, as expected, we observed groups of regulators specific to each NSCLC subtypes. Within the IPA-predicted top 20 upstream regulators, almost half of them were different in SCC and ADC. For the DEGs in SCC, IPA identified 145 activated or inhibited upstream regulators, among which there are regulators such as RARA, KDM1A, HOXA10, EHF, EP400, GATA4, GATA6, and TP63 characteristic only for SCC. KDM1A, with the highest activation z-score, encodes lysine demethylase 1, which contributes to tumorigenesis in many types of cancer [38], while TP63 encodes p63 protein which is a master regulator in epithelial development and has shown significance in cancer progression [39].
For the DEGs in ADC, IPA identified 194 activated or inhibited upstream regulators, among which there are regulators such as FGF2, IL1B, SP1, HDAC1, VEGF, CG, and CTNNB1. FGF2, with the highest z-score in ADC encodes the fibroblast growth factor 2, which is associated with multiple cancer pathogenesis process including cell proliferation, differentiation and metastasis [40]. IL1B, the second upstream regulator with the high z-score in ADC, encodes IL-1β which is one of the inflammatory cytokines, plays a key role in carcinogenesis and tumour progression [41]. The target molecules list for each upstream regulator is provided in Table S1.

2.5. Construction of Weighted Gene Co-Expression Network

WGCNA was performed to construct co-expressed networks and identify co-expression modules. Hierarchical clustering of the samples was performed based on an euclidean distance computed on log10-transformed RNA-seq fractional counts (Figure 6A,B). Additionally, basic patient information was added below the resulting tree. The sample aggregation suggested by the resulting trees does not identify strong outliers, and the basic patient information does not suggest the inclusion of an atypical patient. The research of an adequate soft-threshold validating converged toward scale-free topology, pointing at values 10 and 7 for ADC and SCC, respectively (Figure 6C).
As shown in Figure 7, the first set of modules was obtained with Dynamic Tree Cut algorithm, then correlated modules (r > 0.75) were merged together (Merged dynamic); 32 and 15 modules were identified in ADC and SCC, respectively (Figure 7A,B).
The co-expression network was validated on LUAC (lung adenocarcinoma) and LUSC (lung squamous cell carcinoma) subset from the TCGA database. No outlier was identified, and current soft-threshold was validated. The network was therefore built with identical parameters (Figure S4). The gene intersection matrices showed associations between modules found in TCGA or MOBIT datasets. This was further evidenced by comparing to the results after random permutation of attributed gene clusters in TCGA, which exemplified the expected intersection matrix in the absence of association between TCGA- and MOBIT-based networks. Finally, some modules obtained in TCGA could be strongly associated with corresponding gene module in MOBIT, thus validating MOBIT modules (e.g., TCGA yellowgreen to MOBIT darkturquoise for ADC, or TCGA white to MOBIT blue in SCC) (Figure 8).

2.6. Association of Modules with Clinical Traits

For each module, the gene co-expression was summarised by the Eigengene (i.e., the first component expression of genes belonging to that module) and we calculated the correlations of each Eigengene with clinical traits, such as tumour width, tumour length, SUVmax, BMI (body mass index), pack-year, recurrence-free survival, disease-free survival, T, and smoking status. We noticed that there were four and two key modules most associated with one or more traits in ADC and SCC, respectively. As presented in Figure 9A, the red and lightcyan modules were strongly correlated with the tumour size in ADC, whereas the darkorange module had the highest negative correlation with the recurrence-free survival (RFS) in ADC. Furthermore, the yellow module, as well as the blue and turquoise modules had the highest correlation with SUVmax in ADC and SCC, respectively (Figure 9A,B). These results allowed selection of modules which were interesting for further analysis.

2.7. Functional Enrichment Analysis of Critical Modules

The GO-biological process (BP) and canonical pathway analysis by IPA were performed in four modules (red, lightcyan, darkorange, and yellow) in ADC and two modules (blue and turquoise) in SCC. As shown in Figure 10A, the GO analysis results revealed that in ADC, red and lightcyan modules’ genes in the BP group were mainly enriched in immune response and the regulation of this process. For the darkorange module, which had strong negative correlation with RFS, the genes were mainly enriched in RNA biogenesis and processing. The mitotic cell cycle and DNA conformation change processes were enriched in modules yellow and turquoise in ADC and SCC, respectively. The blue module genes in SCC in BP group were mainly enriched in cell adhesion, communication, and signalling processes.
Following analysis by IPA revealed that in the red and lightcyan modules of ADC, canonical pathways were mainly related to the checkpoint of immune cell activation. The pathways in the darkorange module included the canonical pathways connected with the translation process. The enriched pathways in the yellow module in ADC, as well as in turquoise module in SCC were involved in cell cycle regulation process, involving cyclins and kinases. In SCC, the blue module included the canonical pathways connected with the signalling through G protein-coupled receptors and axon guidance molecules (Figure 10B).

2.8. Hub Gene Identification in the Selected Module

In addition, we also constructed a protein-protein interactions (PPI) network of the expressed genes in the four selected modules in ADC (red, lightcyan, darkorange, yellow), in two selected modules in SCC (blue, turquoise) (Figures S5 and S6), and identified top 10 hub genes (Figure 11). The PPI network was generated using the STRING database in Cytoscape version 3.7.2 (http://cytoscape.org/), and 11 scoring methods including the newly developed algorithm Maximal Clique Centrality (MCC) were employed by use of cytoHubba plugin [42]. Hub genes with high betweenness centrality are important because they have the shortest path connectors through a network.
The top 10 hub genes in the red module, in ADC comprised CTLA4, PAX5, CCR7, SELL, CD19, CD28, TNF, IL10, FOXP3, and IL7R. For the red module, the highly connected hub gene was CTLA4, which is a member of the immunoglobulin superfamily and encodes a transmembrane protein which negatively regulates the T cell function [43]. In the lightcyan module, the top 10 high-degree genes were MZB1, IGJ, CD79A, IRF4, IGLL5, PRDM1, XBP1, CD27, TNRFSF17, and SEL1L. Most of the 10 hub genes in this module encode immunoregulatory factors. For instance, MZB1, with the highest degree of connectivity in the protein-protein interactions (PPI) network, promotes IgM assembly and secretion, as well as is a relevant component of cytoplasm involved in the regulation of apoptosis [44]. NIP7, WDR12, EBNA1BP2, RSL1D1, GNL3, GTPB4, MAK16, NIFK, PNO1 and WDR43 were identified in the darkorange module in ADC. NIP7 encodes the nucleolar pre-rRNA processing protein which plays an important role in rRNA biosynthesis [45]. The top 10 hub genes in the yellow module of ADC included BUB1B, CCNB2, NCAPG, KIF20A, KIF11, CCNB1, UBE2C, CCNA2, AURKB, and PBK. BUB1 encodes a serine/threonine-protein kinase that is a mitotic checkpoint protein [46].
We found the top 10 hub genes in the turquoise module in SCC including CCNB2, KIF20A, KIF23, KIF11, CCNB1, CCNA2, AURKB, MELK, CDK1, and ASPM. CCNB2 encodes cyclin B2 protein, with the highest degree of connectivity in the PPI network in module turquoise in SCC, and may exert its effect by acting as a regulator of G2-M transition during the cell cycle [47,48]. In the blue module, the top 10 high-degree genes were GNG11, ADYC9, ADYC6, ADYC5, ADYC4, C3, CHRM2, PPB2, PF4, and AGT. The hub gene GNG11 which encodes G protein subunit gamma 11 plays an important role in the transmembrane signalling system [49].

3. Discussion

In order to find the unique genes signature in each subtype of lung cancer, we did RNA-seq analysis in ADC and SCC samples compared with normal lung tissue in a well-defined cohort of NSCLC patients. Gene ontology (GO), canonical pathway analysis with the prediction of upstream regulators, and weighted co-expression gene network analysis (WGCNA) with the selection of hub genes were used to explain the most likely significant transcriptome profile alterations and to link them with the clinical traits.
We identified 7720 and 5856 differentially expressed genes in lung ADC and SCC, respectively. The majority of genes were downregulated. Compared with lung adenocarcinoma, there were more DEGs in the lung squamous carcinoma tissues. The top upregulated genes in both subtypes of NSCL play a major part in cancer progression, as evidenced by human cancer studies. Among the DEGs identified in both subtypes, MAGEA9 showed the greatest upregulation in tumour tissues. Moreover, it was also noticed that other genes from MAGEA family members were also significantly higher than normal tissues in both subtypes of NSCLC (MAGEA1, MAGEA10, MAGEA11, MAGEA4). Many studies have shown that melanoma-associated antigen (MAGE) genes, which are the best-characterised members of cancer/testis (CT) antigens family, have no expression in somatic adult tissues, but they differentially expressed in a variety of human cancers [50]. MAGEA9 is frequently overexpressed in bladder cancer, cutaneous T-cell lymphomas, oesophageal adenocarcinomas, as well as in renal cell carcinoma [51,52,53,54]. Currently, there are only a few reports about MAGEA9 in NSCLC. Zhang et al [55] showed that the high expression of MAGEA9 protein in NSCLC tumour cells was present in SCC and it was correlated with tumour size and lymph node metastasis. Our results indicate that the members of MAGEA gene family, homologous MAGEA9 in particular, may play key role in lung ADC and SCC development. If the results from further investigations from in vivo and in vitro studies support our findings, the combination of the appropriate strategy with targeting MAGEA9 might be expected as a high efficacy of the chemo-immunotherapeutic approach. On the other hand, the lowest expression in SCC and ADC tissues was of SFTPC gene. Recent experiments have established that the expression of SFTPC was downregulated in human lung adenocarcinoma tissues and it was correlated with poor overall survival of ADC patients [56]. Overexpression of CFBF, PAGE2, CALM3, DSG3, CTAG2, and FAM83C was seen only in squamous cell carcinoma. The overexpression of PAGE2 and DSG3 was previously described in NSCLC. Djureinovic et al [57] confirmed the overexpression of PAGE2 protein in NSCLC tissues. Moreover, results from microarray-based studies have led to the identification of upregulated expression levels of the DSG3 gene in squamous cell carcinoma.
To further identify the differentiating molecular mechanism of ADC and SCC, we focused on the biological pathways. The genes in both subtypes were mainly involved in biological processes associated with cell division and cell proliferation. However, IPA results have suggested that eicosanoid signalling, as well as agranulocyte and granulocyte adhesion and diapedesis signalling pathways, may be involved in the development of NSCLC. The top activated canonical pathway in lung ADC and SCC was eicosanoid signalling pathway (Figure S1). In general, eicosanoids are biologically active lipids which play key role in the various pathological processes, such as inflammation and cancer. For instance, prostaglandin E2 (PGE2) affects the development of cancer, primarily by suppression of immune response which leads to the induction of cellular proliferation. Interestingly, according to the previous study, there were significant differences in the extent of upregulation of mPGES-1 and its enzymatic product PGE2 was found to reveal lung tumours [58]. Agranulocyte and granulocyte adhesion and diapedesis, which is another pathway identified in the IPA, appears to be involved in the tumour immunosurveillance (Figures S2 and S3).
One of the mechanisms of avoiding immune surveillance by tumour is the suppression of the activity of adhesion molecules in the multistep adhesion cascade during the emigration of leukocytes from the bloodstream into the tissue [59]. In our study, we found that a variety of genes attending in the adhesion cascade (CD34, E-selectin, P-selectin, ICAM1, ICAM2, JAM2, PNAd, chemokine) were downregulated in lung SCC and ADC. The previous study evaluated the low expression of P-selectin protein correlated with diminished leukocyte tumour infiltration and poor prognosis [60]. Similarly, it was indicated that ICAM-1 and ICAM-2 suppression was positively correlated with decreased leukocyte extravasation and associated with cancer development [61]. In addition, a previous study with renal cell carcinoma by Hellwig et al [62] revealed that the expression of CD34 was downregulated compared to normal renal tissue. Thus, our results indicated that NSCLC is significantly associated with mechanisms to escape immunity by tumour cells.
Many of the lung SCC-associated altered biological processes were involved in keratinisation, epidermal cell and keratinocyte differentiation. SCC is a malignant epithelial tumour, arising in tissues that provide a barrier between the organism and its environment [63]. At an earlier stage, lung SCC is accompanied by epidermal keratinisation. Park et al [64] have provided ample evidence that keratinisation is a significantly poor prognosis factor in lung SCC. Tumour with keratinisation was also more frequently observed among patients with a history of smoking. In addition, we found in our study that the top upstream regulators in lung SCC were KDM1A and TP63 genes. KDM1A is overexpressed in numerous epithelial cancers and promotes proliferation, migration and invasion. Egolf et al [65] reported that its downstream targets were involved in epidermal differentiation and cornification processes. Kong et al [66] indicated that treatment with KDM1A inhibitor (2-PCPA) reduced tumour NSCLC cell proliferation and migration. On the other hand, TP63 has been shown as a key player in homeostasis and development of squamous epithelium. Moreover, numerous investigations confirmed the upregulation of expression of ΔNp63α, the predominant TP63 isoform in squamous epithelium, in squamous cell cancers, including head and neck, skin, as well as the lung. Cancer-related target molecules of TP63 identified by IPA (Table S1) were correlated with the apoptotic process, keratinisation, tissue development, and cell population proliferation. There is evidence that TP63 can promote cell survival and proliferation, and is involved in cancer formation and progression [39].
The top activated GO-terms and canonical pathways in ADC detected a unique lung ADC landscape wherein deregulation of certain pathways involved targets associated specifically with tumour immune response. Our results indicated that the MIF regulation of innate immunity and eNOS signalling were the most significant changed pathways in ADC. MIF as a proinflammatory cytokine was determined to be upregulated in all stages of neoplasia in most types of cancers and metastatic conditions. It has been revealed that MIF directly links inflammation and initiation of tumorigenesis [67,68]. Interestingly, the coordinated COX-MIF-p53 axis affects tumour suppression downregulation, as well as COX-2 and PGE2 upregulation, finally leading to enhanced tumour growth, proliferation, and progression [69]. This crosstalk can be seen from Figure S3. Nitric oxide (NO) is an important cell signalling molecule whose level was frequently elevated in many tumours, including that of the lung [70]. Collected evidence has demonstrated the tumour-promoting role of NO in various types of cancer, including of the lung, spanning from tumour initiation of cellular transformation to tumour progression, through the metastatic cascade and resistance to radio/chemotherapy. It was also further demonstrated that high levels of serum NO were associated with advanced-stage lung cancer and poor survival rate of patients. Among numerous upstream regulators in lung ADC, we identified those that were associated with tumorigenesis and inflammation (Table S2). One of them was FGF2 which affected downstream molecules like BIRC5, CDC25A, CXCL2, CXCL12, CDC25A, AXIN2, MMP1, MMP9, MMP13, and MIF. These molecules are related to a progression in various types of cancer. For instance, in small cell lung cancer, FGF2 increased the expression of anti-apoptotic proteins and triggered chemoresistance [71]. The activation of FGF receptors can activate multiple signal transduction pathways including PI3K, MAPK, and STAT pathway in prostate cancer, thereby causing tumour progression [72]. In addition, IL1B pro-inflammatory cytokine turned out to be a promising upstream regulator in lung ADC. It encodes the IL-1B protein which affects multiple aspects of the tumour microenvironment. IL-1B generated in the course of chronic inflammation supports tumour growth and metastasis predominantly by tumour-infiltrating macrophages [41,73].
Despite the fact that the traditional DEGs analysis has provided enormously relevant information, only the application of weighted co-expression network analysis (WGCNA) allowed for identifying the pattern of correlation among genes. WGCNA is notably useful for the identification of the modules of co-expressed genes that are correlated with clinical traits and consequently biological tumour behaviour. We identified four modules (red, lightcyan, darkorange, yellow) in ADC, and 2 modules (blue, turquoise) in SCC, the most related with tumour size (width, length), SUVmax parameter and recurrence-free survival (RFS). We found strong correlations of module red and lightcyan in ADC with tumour size. In addition, the results showed that the darkorange module in ADC was significantly associated with RFS. These results demonstrated that key genes within the above modules may serve as potential markers of tumour progression. We determined that the genes in the red and lightcyan module were enriched in the biological pathways involved in the regulation of immune response, complement activation, and altered T cells and B cells signalling. These key signalling pathways are closely related to the tumour microenvironment [74]. Similarly, Choi et al. [75] identified through WGCNA magenta and brown module of co-expressed genes in lung adenocarcinoma which were significantly correlated with tumour metabolism and immune microenvironment. The surge of interest in cancer immunotherapy is mainly focused on tumour immunosurveillance. Numerous studies have noted the tumour-infiltration of immune cells in the lung cancer microenvironment, but these cells play a dualistic role, they suppress and/or promote tumour progression and growth based on the production of growth factors, cytokines and chemokines into the tumour microenvironment [76].
To identify key genes in the five different colour modules, a PPI network was constructed and hub genes were identified with Cytoscape software and cytoHubba plugin, respectively. We found that MZB1 was the top high-degree gene in ADC lightcyan module. According to previous studies, MZB1, an immunoregulatory molecule, was expressed at significantly decreased levels in gastric cancer (GC) tissues. It has been reported that siRNA-mediated knockdown of MZB1 significantly increased proliferation, invasion and migration of GC cell lines [77]. The WGCNA used by Zhai et al [78] also demonstrated that MZB1 was the hub gene in the module strongly associated with immune response in lung adenocarcinoma. According to the network analysis of the red module, CTLA4 and PAX5 were identified as the highest degree hub genes. Until recently, CTLA-4 was known as the inhibitor of anti-tumour immunity, but the effectiveness of the anti-CTLA-4 therapies is limited to a small number of cancer types [79]. However, it has been recently shown that CTLA4 may play the tumour suppressor role across different cancer, including lung adenocarcinoma [80]. The protein encoded by PAX5 acts as a nuclear transcription factor that is involved in B cell development. However, it is shown recently that PAX5 was down-regulated significantly in lung adenocarcinoma but not squamous carcinoma [81,82]. Zhao et al [83] revealed that PAX5 functions as a tumour suppressor in NSCLC cells, both in in vitro and in vivo assay. Interestingly, it was shown that re-expression of PAX5 in lung cells inhibited β-catenin signalling pathway, which plays important role not only in proliferation, differentiation, and migration but also in regulating immune cell infiltration of the tumour microenvironment [84]. For the darkorange module of ADC, the most correlated with recurrence-free survival, NIP7 was the highest degree gene. Interestingly, NIP7 has been also identified as the gene among in a three-gene signature significantly associated with RFS in liposarcoma [85].
Next, we observed that the yellow module of ADC and the turquoise module of SCC were enriched in mitotic cell cycle and nuclear division biological processes and positively correlated with SUVmax. SUVmax is the most common semi-quantitative parameter determined by 18F-FDG PET/MRI. SUVmax of a tumour is a product of several fundamental factors including glucose metabolism and the type or a number of cells present in the tumour. Kosal et al [86] reported that SUVmax of NSCLC tumour was positively correlated with the number of mitosis, the largest tumour diameter, and pathological stage of the disease. In consistence with those studies, the present study found that SUVmax was also positively correlated with genes involved in mitosis and cell cycle regulation process. BUB1 was the highest connective gene in the PPI network in the yellow module of ADC, whereas CCNB2 in the turquoise module of SCC. Both are mitotic genes activated during mitosis. CCNB2 which encodes the cyclin B2 protein is additionally involved in checkpoint control and its overexpression was positively associated with the tumour size and lymph and node metastasis in NSCLC [47,87]. The presence of top 10 genes from the yellow module of ADC and the turquoise module of SCC has been also identified as hub genes in NSCLC by Li et al [88] and Xu et al [89] investigations. Li and co-workers determined through WGCNA the survival-related module containing genes enriched in the cell cycle. In addition, these modules were defined as hub genes such as MELK and UBE2C, which were also identified in the yellow and turquoise module in our studies. Xu et al also observed MELK, CCNB1, CCNA2, AURKA, and KIF11 as key genes associated with survival, which may be potential biomarkers or therapeutic targets for ADC. Furthermore, AURKB, KIF23, and CCNA2 from the turquoise module of SCC were determined as hub genes by Gao et al through WGCNA in modules of lung squamous cell carcinoma that were strongly associated with cell division and overall survival [90].
In the blue module of SCC, which was also correlated with SUVmax parameter, the highest-degree gene was the GNG11. Low expression of GNG11 was associated with worse overall survival for the female lung cancer patients who never smoked [91,92]. However, the specific biological function of GNG11 in lung cancer still remains unknown.

4. Materials and Methods

4.1. Study Cohort

For this study, 114 patients underwent tumour resection surgery, during which two to three tissue fragments were collected. Pathological specimens were reviewed to confirm diagnosis of ADC, SCC or non-malignant lung tissue. Other inclusion criteria consisted of the material containing at least 50% of tumour cells for the RNA extraction; completely resected tumour (free resection margins); and follow-up including monitoring events for cancer recurrence, and lung cancer-related deaths. All patients underwent a PET/MRI examination. All study participants provided written informed consent and received detailed information on the study and associated risk prior to the enrollment. The study was conducted in accordance with the Declaration of Helsinki as a statement of ethical principles for medical research involving human subject and was approved by the local ethics committee of the Medical University of Bialystok, Poland (approval number: R-I-002/36/2014).

4.2. RNA Sample Preparation and Sequencing

Total RNA was purified from 114 tumour and 114 normal tissues with the RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer’s instruction. RNA concentration, purity, and integrity were assessed by Qubit (Invitrogen, Carlsbad, CA, USA) and Tape Station 2200 (Agilent Technologies, Santa Clara, CA, USA). RNA-seq libraries were constructed from 1 μg of total RNA with RNA integrity number (RIN) >8, and using the Illumina TruSeq® Stranded Total RNA Library Prep Gold (#20020598, Illumina, San Diego, CA, USA). Indexed libraries were pooled, clustered with the use of the cBot, and sequenced on the Illumina HiSeq 4000 platform generating 150 bp paired-end reads (2 × 75 bp). The quality of the obtained reads was assessed using FastQC (Babraham Institute, Cambridge, United Kingdom) before the analysis, as well as after different processing steps. As all reads had high quality, they were only soft trimmed and filtered with BBduk (US Department of Energy Joint Genome Institute, Walnuk Creek, CA, USA). BBduk tool was also used to filter out fragments originated from rRNAs. Cleaned reads were mapped to the human reference genome GRCh37 using splice-aware aligner STAR [93]. Counts per gene were estimated with RNA-Seq by Expectation Maximization version 1.3 [94] and used for data visualisation and further differential gene expression analysis. The quality of mapping was assessed with Qualimap version 2.2.1 [95], flagstat module from Samtools [96], and RSeQC [97]. Mapped reads were also visually inspected with Integrative Genomics Viewer [98]. Samples size factors were estimated, and counts were normalised using Relative Log Expression (RLE). Differential expression analysis assuming negative binominal distribution was performed with R packages DESeq2 [19], while the package EdgeR [23] (Trimmed Mean of M-values normalisation) was also used for verification of estimation robustness. After model fitting, dispersion estimates were obtained, and general linearised model was applied. Likehood ratio test (LTR) was used to call differentially expressed (DEGs) under |log2FC| ≥ 1 and false discovery rate (FDR) corrected p-value ≤ 0.05. Volcano plots generated were based on DESeq2 estimation and drawn using the R package ggplot2 v.3.2.1 (Hadley Wickham, RStudio Inc., Boston, MA, USA).

4.3. Functional Enrichment Analysis and Identification of Upstream Regulators

To understand the functional significance of dysregulated genes, the gene ontology (GO) enrichment analysis-biological processes (BP) by BiNGO tool [32] as plugin for Cytoscape version 3.7.2 (http://cytoscape.org/), the Ingenuity Pathway Analysis (IPA) (QIAGEN Inc., https://qiagenbioinformatics.com/products/ingenuity-pathway-analysis), and Metascape tool (http://metascape.org) [33] were used. BiNGO plug-in determined overrepresentative GO categories in a set of genes, based on a two-sided hypergeometric test with a Benjamini-Hochberg adjustment. IPA was used to perform the core analysis to identify canonical pathways and to predict the upstream regulators. Metascape is a web-based portal that provides a comprehensive gene list annotation based on 40 independent databases for studying functional enrichment. Venn diagrams and heat maps were created with the BioVinci Software (BioTuring, San Diego, CA, USA).

4.4. Construction of Weighted Gene Co-Expression Networks and Identification of Modules Associated with Clinical Traits

Gene co-expression network analysis was specifically performed on tumour tissues using the R package WGCNA [11,12]. The expression matrix was restricted to only expressed genes, when the number of reads sequenced was higher than 10 in a minimum of 15 samples, normalised for sample depth (count per million read, CPM) and log-transformed, including a pseudo-count of four (Log10(CPM + 4)). Then, the optimal soft threshold for adjacency computation was graphically determined. The transformed expression matrix was inputted into the WGCNA package functions, modules, and corresponding Eigengenes were obtained. The cutreeDynamic function was used for tree pruning of the gene hierarchical clustering dendrograms resulting in co-expression modules; correlated modules (r > 0.75) were then merged. The dissimilarity of the module Eigengenes (ME) was calculated using the moduleEigengenes function in the R WGCNA package. Association between Eigengenes values with clinical traits was assessed by Pearson’s correlation.

4.5. Protein-Protein Interaction Network Construction for Selected Modules and Hub Genes Identification

In order to identify highly connected hub genes in the protein-protein network (PPI), the cytoHubba plugin based on Cytoscape was used. At first, the correlation of the gene to a module Eigengenes, which is known as the module membership (MM), was assessed by Pearson’s correlation. In this study, hub genes were selected for MM > 0.55 in the specific module.

5. Conclusions

Taken together, our RNA-seq data has generated a wealth of information for potentially new biomarkers, networks, and pathways dysregulated in lung ADC and SCC tissues which can be further explored to unravel the molecular mechanisms regulating the development of each subtype of NSCLC. Key biological processes, canonical pathways, and upstream regulators, which all play important roles in cancer progression, were predicted by bioinformatics tools. In addition, we applied WGCNA for exploring molecular networks associated with clinical traits such as tumour size, SUVmax, BMI, smoking status, recurrence-free survival, and disease-free survival. Overall, the results of this study emphasise the need to further explore the selected genes and pathways, particularly those related to the tumour microenvironment and immune escape mechanisms in ADC and SCC.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-6694/12/1/37/s1, Figure S1. Eicosanoind signalling pathway in lung SCC and ADC identified by IPA, Figure S2. Agranulocyte adhesion and diapedesis pathway in lung SCC and ADC identified by IPA, Figure S3. MIF regulation of innate immunity pathway in lung ADC identified by IPA, Figure S4. Cluster dendograms of the gene clusters of (A) LUAD and (B) LUSC subset from TCGA database, Figure S5. Protein-protein interaction (PPI) networks of genes in the red (A), lightcyan (B), darkorange (C), and yellow (D) modules in lung ADC. The networks were constructed using Cytoscape v. 3.7.2. software, Figure S6. Protein-protein interaction (PPI) networks of genes in the blue (D) and turquoise (E) modules in lung SCC. The networks were constructed using Cytoscape v. 3.7.2. software, Table S1. Upstream regulator analysis of DEGs in lung SCC predicted by IPA, Table S2. Upstream regulator analysis of DEGs in lung ADC predicted by IPA.

Author Contributions

Conceptualisation, M.N., M.K., and A.K.; methodology, M.N., F.C., M.K., and J.R.; formal analysis, M.N., F.C., and K.C.; investigation, M.N., A.S., A.B., and J.R.; writing-original draft preparation, M.N.; writing-review and editing, A.K. and M.K.; visualisation, M.N. and F.C.; supervision, J.N. and A.K. All authors have read and agreed to the published version of the manuscript.

Funding

The project funded by the National Centre for Research and Development in the framework of Programme “Prevention practices and treatment of civilization diseases”—STRATEGMED (contract no. STRATEGMED2/266484/2/NCBR/2015).

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Siegel, R.L.; Miller, K.D.; Jemal, A. Cancer statistics, 2019: {Cancer} {Statistics}, 2019. CA. Cancer J. Clin. 2019, 69, 7–34. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Herbst, R.S.; Morgensztern, D.; Boshoff, C. The biology and management of non-small cell lung cancer. Nature 2018, 553, 446–454. [Google Scholar] [CrossRef] [PubMed]
  3. Herbst, R.S.; Heymach, J.V.; Lippman, S.M. Lung cancer. N. Engl. J. Med. 2008, 359, 1367–1380. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Faruki, H.; Mayhew, G.M.; Serody, J.S.; Hayes, D.N.; Perou, C.M.; Lai-Goldman, M. Lung Adenocarcinoma and Squamous Cell Carcinoma Gene Expression Subtypes Demonstrate Significant Differences in Tumor Immune Landscape. J. Thorac. Oncol. 2017, 12, 943–953. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Gentles, A.J.; Bratman, S.V.; Lee, L.J.; Harris, J.P.; Feng, W.; Nair, R.V.; Shultz, D.B.; Nair, V.S.; Hoang, C.D.; West, R.B.; et al. Integrating Tumor and Stromal Gene Expression Signatures With Clinical Indices for Survival Stratification of Early-Stage Non-Small Cell Lung Cancer. J. Natl. Cancer Inst. 2015, 107, djv211. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Charkiewicz, R.; Niklinski, J.; Claesen, J.; Sulewska, A.; Kozlowski, M.; Michalska-Falkowska, A.; Reszec, J.; Moniuszko, M.; Naumnik, W.; Niklinska, W. Gene Expression Signature Differentiates Histology But Not Progression Status of Early-Stage NSCLC. Transl. Oncol. 2017, 10, 450–458. [Google Scholar] [CrossRef]
  7. Wilkerson, M.D.; Yin, X.; Walter, V.; Zhao, N.; Cabanski, C.R.; Hayward, M.C.; Miller, C.R.; Socinski, M.A.; Parsons, A.M.; Thorne, L.B.; et al. Differential pathogenesis of lung adenocarcinoma subtypes involving sequence mutations, copy number, chromosomal instability, and methylation. PLoS ONE 2012, 7, e36530. [Google Scholar] [CrossRef] [Green Version]
  8. Ozsolak, F.; Milos, P.M. RNA sequencing: Advances, challenges and opportunities. Nat. Rev. Genet. 2011, 12, 87–98. [Google Scholar] [CrossRef]
  9. Ramani, V.; Qiu, R.; Shendure, J. High-throughput determination of RNA structure by proximity ligation. Nat. Biotechnol. 2015, 33, 980–984. [Google Scholar] [CrossRef] [Green Version]
  10. Barabási, A.-L.; Gulbahce, N.; Loscalzo, J. Network medicine: A network-based approach to human disease. Nat. Rev. Genet. 2011, 12, 56–68. [Google Scholar] [CrossRef] [Green Version]
  11. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, e559. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Langfelder, P.; Horvath, S. Fast R Functions for Robust Correlations and Hierarchical Clustering. J. Stat. Softw. 2012, 46. [Google Scholar] [CrossRef] [Green Version]
  13. Tian, F.; Zhao, J.; Fan, X.; Kang, Z. Weighted gene co-expression network analysis in identification of metastasis-related genes of lung squamous cell carcinoma based on the Cancer Genome Atlas database. J. Thorac. Dis. 2017, 9, 42–53. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Wang, C.C.N.; Li, C.Y.; Cai, J.-H.; Sheu, P.C.-Y.; Tsai, J.J.P.; Wu, M.-Y.; Li, C.-J.; Hou, M.-F. Identification of Prognostic Candidate Genes in Breast Cancer by Integrated Bioinformatic Analysis. J. Clin. Med. 2019, 8, e1160. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Yin, L.; Cai, Z.; Zhu, B.; Xu, C. Identification of Key Pathways and Genes in the Dynamic Progression of HCC Based on WGCNA. Genes (Basel) 2018, 9, e92. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Yin, X.; Wang, J.; Zhang, J. Identification of biomarkers of chromophobe renal cell carcinoma by weighted gene co-expression network analysis. Cancer Cell Int. 2018, 18, e206. [Google Scholar] [CrossRef]
  17. Zhou, X.; Huang, X.; Liang, S.; Tang, S.; Wu, S.; Huang, T.; Mo, Z.; Wang, Q. Identifying miRNA and gene modules of colon cancer associated with pathological stage by weighted gene co-expression network analysis. Onco Targets. Ther. 2018, 11, 2815–2830. [Google Scholar] [CrossRef] [Green Version]
  18. Nakamura, H.; Fujii, K.; Gupta, V.; Hata, H.; Koizumu, H.; Hoshikawa, M.; Naruki, S.; Miyata, Y.; Takahashi, I.; Miyazawa, T.; et al. Identification of key modules and hub genes for small-cell lung carcinoma and large-cell neuroendocrine lung carcinoma by weighted gene co-expression network analysis of clinical tissue-proteomes. PLoS ONE 2019, 14, e0217105. [Google Scholar] [CrossRef]
  19. Niklinski, J.; Kretowski, A.; Moniuszko, M.; Reszec, J.; Michalska-Falkowska, A.; Niemira, M.; Ciborowski, M.; Charkiewicz, R.; Jurgilewicz, D.; Kozlowski, M.; et al. Systematic biobanking, novel imaging techniques, and advanced molecular analysis for precise tumor diagnosis and therapy: The Polish MOBIT project. Adv. Med. Sci. 2017, 62, 405–413. [Google Scholar] [CrossRef]
  20. Marko-Varga, G. BioBanking as the central tool for translational medicine CTM issue 2013. Clin. Transl. Med. 2013, 2, e4. [Google Scholar] [CrossRef] [Green Version]
  21. Neumeister, V.M.; Juhl, H. Tumor Pre-Analytics in Molecular Pathology: Impact on Protein Expression and Analysis. Curr. Pathobiol. Rep. 2018, 6, 265–274. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. David, K.A.; Unger, F.T.; Uhlig, P.; Juhl, H.; Moore, H.M.; Compton, C.; Nashan, B.; Dörner, A.; de Weerth, A.; Zornig, C. Surgical procedures and postsurgical tissue processing significantly affect expression of genes and EGFR-pathway proteins in colorectal cancer tissue. Oncotarget 2014, 5, 11017–11028. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  23. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010, 26, 139–140. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, e550. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Weon, J.L.; Potts, P.R. The MAGE protein family and cancer. Curr. Opin. Cell Biol. 2015, 37, 1–8. [Google Scholar] [CrossRef] [Green Version]
  26. Zhang, X.; Kong, M.; Zhang, Z.; Xu, S.; Yan, F.; Wei, L.; Zhou, J. FGF19 genetic amplification as a potential therapeutic target in lung squamous cell carcinomas. Thorac. Cancer 2017, 8, 655–665. [Google Scholar] [CrossRef]
  27. Zhang, Y.Q.; Zhang, J.J.; Song, H.J.; Li, D.W. Overexpression of CST4 promotes gastric cancer aggressiveness by activating the ELFN2 signaling pathway. Am. J. Cancer Res. 2017, 7, 2290–2304. [Google Scholar]
  28. Yin, Y.; Betsuyaku, T.; Garbow, J.R.; Miao, J.; Govindan, R.; Ornitz, D.M. Rapid induction of lung adenocarcinoma by fibroblast growth factor 9 signaling through FGF receptor 3. Cancer Res. 2013, 73, 5730–5741. [Google Scholar] [CrossRef] [Green Version]
  29. Shang, J.; Wang, Z.; Chen, W.; Yang, Z.; Zheng, L.; Wang, S.; Li, S. Pseudogene CHIAP2 inhibits proliferation and invasion of lung adenocarcinoma cells by means of the WNT pathway. J. Cell. Physiol. 2019, 234, 13735–13746. [Google Scholar] [CrossRef]
  30. Zhai, S.; Xu, H.; Jiang, X.; Sun, S.; Qin, L.; Wei, D.; Hu, L. Expression Depression of CD300LG-γ in Human Pulmonary Carcinoma. Monoclon. Antib. Immunodiagn. Immunother. 2016, 35, 94–99. [Google Scholar] [CrossRef]
  31. Schürch, D.; Ospina, O.L.; Cruz, A.; Pérez-Gil, J. Combined and Independent Action of Proteins SP-B and SP-C in the Surface Behavior and Mechanical Stability of Pulmonary Surfactant Films. Biophys. J. 2010, 99, 3290–3299. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Maere, S.; Heymans, K.; Kuiper, M. BiNGO: A Cytoscape plugin to assess overrepresentation of Gene Ontology categories in Biological Networks. Bioinformatics 2005, 21, 3448–3449. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Zhou, Y.; Zhou, B.; Pache, L.; Chang, M.; Khodabakhshi, A.H.; Tanaseichuk, O.; Benner, C.; Chanda, S.K. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun. 2019, 10, e1523. [Google Scholar] [CrossRef] [PubMed]
  34. Krämer, A.; Green, J.; Pollard, J.; Tugendreich, S. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 2014, 30, 523–530. [Google Scholar] [CrossRef]
  35. Chuang, J.C.; Stehr, H.; Liang, Y.; Das, M.; Huang, J.; Diehn, M.; Wakelee, H.A.; Neal, J.W. ERBB2-Mutated Metastatic Non-Small Cell Lung Cancer: Response and Resistance to Targeted Therapies. J. Thorac. Oncol. 2017, 12, 833–842. [Google Scholar] [CrossRef] [Green Version]
  36. Liu, S.; Li, S.; Hai, J.; Wang, X.; Chen, T.; Quinn, M.M.; Gao, P.; Zhang, Y.; Ji, H.; Cross, D.A.E.; et al. Targeting HER2 Aberrations in Non–Small Cell Lung Cancer with Osimertinib. Clin. Cancer Res. 2018, 24, 2594–2604. [Google Scholar] [CrossRef] [Green Version]
  37. LaBaer, J.; Garrett, M.D.; Stevenson, L.F.; Slingerland, J.M.; Sandhu, C.; Chou, H.S.; Fattaey, A.; Harlow, E. New functional activities for the p21 family of CDK inhibitors. Genes Dev. 1997, 11, 847–862. [Google Scholar] [CrossRef] [Green Version]
  38. Hayami, S.; Kelly, J.D.; Cho, H.-S.; Yoshimatsu, M.; Unoki, M.; Tsunoda, T.; Field, H.I.; Neal, D.E.; Yamaue, H.; Ponder, B.A.J.; et al. Overexpression of LSD1 contributes to human carcinogenesis through chromatin regulation in various cancers. Int. J. Cancer 2011, 128, 574–586. [Google Scholar] [CrossRef]
  39. Bergholz, J.; Xiao, Z.-X. Role of p63 in Development, Tumorigenesis and Cancer Progression. Cancer Microenviron. 2012, 5, 311–322. [Google Scholar] [CrossRef] [Green Version]
  40. Grose, R.; Dickson, C. Fibroblast growth factor signaling in tumorigenesis. Cytokine Growth Factor Rev. 2005, 16, 179–186. [Google Scholar] [CrossRef]
  41. Bent, R.; Moll, L.; Grabbe, S.; Bros, M. Interleukin-1 Beta-A Friend or Foe in Malignancies? Int. J. Mol. Sci. 2018, 19, e2155. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Chin, C.-H.; Chen, S.-H.; Wu, H.-H.; Ho, C.-W.; Ko, M.-T.; Lin, C.-Y. cytoHubba: Identifying hub objects and sub-networks from complex interactome. BMC Syst. Biol. 2014, 8, S11. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Buchbinder, E.I.; Desai, A. CTLA-4 and PD-1 Pathways: Similarities, Differences, and Implications of Their Inhibition. Am. J. Clin. Oncol. 2016, 39, 98–106. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Flach, H.; Rosenbaum, M.; Duchniewicz, M.; Kim, S.; Zhang, S.L.; Cahalan, M.D.; Mittler, G.; Grosschedl, R. Mzb1 Protein Regulates Calcium Homeostasis, Antibody Secretion, and Integrin Activation in Innate-like B Cells. Immunity 2010, 33, 723–735. [Google Scholar] [CrossRef] [Green Version]
  45. Stålberg, P.; Grimfjärd, P.; Santesson, M.; Zhou, Y.; Lindberg, D.; Gobl, A.; Öberg, K.; Westin, G.; Rastad, J.; Wang, S.; et al. Transfection of the Multiple Endocrine Neoplasia Type 1 Gene to a Human Endocrine Pancreatic Tumor Cell Line Inhibits Cell Growth and Affects Expression of JunD, δ-Like Protein 1/Preadipocyte Factor-1, Proliferating Cell Nuclear Antigen, and QM/Jif-1. J. Clin. Endocrinol. Metab. 2004, 89, 2326–2337. [Google Scholar] [CrossRef] [Green Version]
  46. Bolanos-Garcia, V.M.; Blundell, T.L. BUB1 and BUBR1: Multifaceted kinases of the cell cycle. Trends Biochem. Sci. 2011, 36, 141–150. [Google Scholar] [CrossRef] [Green Version]
  47. Giotti, B.; Chen, S.-H.; Barnett, M.W.; Regan, T.; Ly, T.; Wiemann, S.; Hume, D.A.; Freeman, T.C. Assembly of a parts list of the human mitotic cell cycle machinery. J. Mol. Cell Biol. 2019, 11, 703–718. [Google Scholar] [CrossRef]
  48. Ni, Z.; Wang, X.; Zhang, T.; Li, L.; Li, J. Comprehensive analysis of differential expression profiles reveals potential biomarkers associated with the cell cycle and regulated by p53 in human small cell lung cancer. Exp. Ther. Med. 2018, 15, 3273–3282. [Google Scholar] [CrossRef]
  49. Hossain, M.N.; Sakemura, R.; Fujii, M.; Ayusawa, D. G-protein gamma subunit GNG11 strongly regulates cellular senescence. Biochem. Biophys. Res. Commun. 2006, 351, 645–650. [Google Scholar] [CrossRef]
  50. Meek, D.W.; Marcar, L. MAGE-A antigens as targets in tumour therapy. Cancer Lett. 2012, 324, 126–132. [Google Scholar] [CrossRef]
  51. Bergeron, A.; Picard, V.; LaRue, H.; Harel, F.; Hovington, H.; Lacombe, L.; Fradet, Y. High frequency of MAGE-A4 and MAGE-A9 expression in high-risk bladder cancer. Int. J. Cancer 2009, 125, 1365–1371. [Google Scholar] [CrossRef] [PubMed]
  52. Eichmüller, S.; Usener, D.; Thiel, D.; Schadendorf, D. Tumor-specific antigens in cutaneous T-cell lymphoma: Expression and sero-reactivity. Int. J. Cancer 2003, 104, 482–487. [Google Scholar] [CrossRef] [PubMed]
  53. Lin, J.; Lin, L.; Thomas, D.G.; Greenson, J.K.; Giordano, T.J.; Robinson, G.S.; Barve, R.A.; Weishaar, F.A.; Taylor, J.M.G.; Orringer, M.B.; et al. Melanoma-Associated Antigens in Esophageal Adenocarcinoma: Identification of Novel MAGE-A10 Splice Variants. Clin. Cancer Res. 2004, 10, 5708–5716. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  54. Hatiboglu, G.; Pritsch, M.; Macher-Goeppinger, S.; Zöller, M.; Huber, J.; Haferkamp, A.; Pahernik, S.; Wagener, N.; Hohenfellner, M. Prognostic value of melanoma-associated antigen A9 in renal cell carcinoma. Scand. J. Urol. 2013, 47, 311–322. [Google Scholar] [CrossRef]
  55. Zhang, S.; Zhai, X.; Wang, G.; Feng, J.; Zhu, H.; Xu, L.; Mao, G.; Huang, J. High expression of MAGE-A9 in tumor and stromal cells of non-small cell lung cancer was correlated with patient poor survival. Int. J. Clin. Exp. Pathol. 2015, 8, 541–550. [Google Scholar]
  56. Li, B.; Meng, Y.-Q.; Li, Z.; Yin, C.; Lin, J.-P.; Zhu, D.-J.; Zhang, S.-B. MiR-629-3p-induced downregulation of SFTPC promotes cell proliferation and predicts poor survival in lung adenocarcinoma. Artif. Cells Nanomed. Biotechnol. 2019, 47, 3286–3296. [Google Scholar] [CrossRef] [Green Version]
  57. Djureinovic, D.; Hallström, B.M.; Horie, M.; Mattsson, J.S.M.; La Fleur, L.; Fagerberg, L.; Brunnström, H.; Lindskog, C.; Madjar, K.; Rahnenführer, J.; et al. Profiling cancer testis antigens in non-small-cell lung cancer. JCI Insight 2016, 1, e86837. [Google Scholar] [CrossRef] [Green Version]
  58. Wu, Y.-C.; Su, L.-J.; Wang, H.-W.; Jeff Lin, C.-F.; Hsu, W.-H.; Chou, T.-Y.; Huang, C.-Y.F.; Lu, C.-L.; Hsueh, C.-T. Co-Overexpression of Cyclooxygenase-2 and Microsomal Prostaglandin E Synthase-1 Adversely Affects the Postoperative Survival in Non-small Cell Lung Cancer. J. Thorac. Oncol. 2010, 5, 1167–1174. [Google Scholar] [CrossRef] [Green Version]
  59. Strell, C.; Entschladen, F. Extravasation of leukocytes in comparison to tumor cells. Cell Commun. Signal. 2008, 6, e10. [Google Scholar] [CrossRef] [Green Version]
  60. Peeters, C.F.J.M.; Ruers, T.J.M.; Westphal, J.R.; de Waal, R.M.W. Progressive loss of endothelial P-selectin expression with increasing malignancy in colorectal cancer. Lab. Investig. 2005, 85, 248–256. [Google Scholar] [CrossRef] [Green Version]
  61. Steege, J.C.A.B.-T.; Baeten, C.I.M.; Thijssen, V.L.J.L.; Satijn, S.A.; Verhoeven, I.C.L.; Hillen, H.F.P.; Wagstaff, J.; Griffioen, A.W. Angiogenic Profile of Breast Carcinoma Determines Leukocyte Infiltration. Clin. Cancer Res. 2004, 10, 7171–7178. [Google Scholar] [CrossRef] [Green Version]
  62. Hellwig, S.M.; Damen, C.A.; van Adrichem, N.P.; Blijham, G.H.; Groenewegen, G.; Griffioen, A.W. Endothelial CD34 is suppressed in human malignancies: Role of angiogenic factors. Cancer Lett. 1997, 120, 203–211. [Google Scholar] [CrossRef]
  63. Languino, L.R.; Singh, A.; Prisco, M.; Inman, G.J.; Luginbuhl, A.; Curry, J.M.; South, A.P. Exosome-mediated transfer from the tumor microenvironment increases TGFβ signaling in squamous cell carcinoma. Am. J. Transl. Res. 2016, 8, 2432–2437. [Google Scholar]
  64. Park, H.J.; Cha, Y.-J.; Kim, S.H.; Kim, A.; Kim, E.Y.; Chang, Y.S. Keratinization of Lung Squamous Cell Carcinoma Is Associated with Poor Clinical Outcome. Tuberc. Respir. Dis. (Seoul) 2017, 80, 179–186. [Google Scholar] [CrossRef]
  65. Egolf, S.; Aubert, Y.; Doepner, M.; Anderson, A.; Maldonado-Lopez, A.; Pacella, G.; Lee, J.; Ko, E.K.; Zou, J.; Lan, Y.; et al. LSD1 Inhibition Promotes Epithelial Differentiation through Derepression of Fate-Determining Transcription Factors. Cell Rep. 2019, 28, 1981–1992.e7. [Google Scholar] [CrossRef] [Green Version]
  66. Kong, L.; Zhang, P.; Li, W.; Yang, Y.; Tian, Y.; Wang, X.; Chen, S.; Yang, Y.; Huang, T.; Zhao, T.; et al. KDM1A promotes tumor cell invasion by silencing TIMP3 in non-small cell lung cancer cells. Oncotarget 2016, 7, 27959–27974. [Google Scholar] [CrossRef] [Green Version]
  67. Grieb, G.; Merk, M.; Bernhagen, J.; Bucala, R. Macrophage migration inhibitory factor (MIF): A promising biomarker. Drug News Perspect. 2010, 23, e257. [Google Scholar] [CrossRef] [Green Version]
  68. Hu, C.-T.; Guo, L.-L.; Feng, N.; Zhang, L.; Zhou, N.; Ma, L.-L.; Shen, L.; Tong, G.; Yan, Q.-W.; Zhu, S.; et al. MIF, secreted by human hepatic sinusoidal endothelial cells, promotes chemotaxis and outgrowth of colorectal cancer in liver prometastasis. Oncotarget 2015, 6, 22410–22423. [Google Scholar] [CrossRef] [Green Version]
  69. Dave, M.; Islam, A.B.M.M.K.; Jensen, R.V.; Rostagno, A.; Ghiso, J.; Amin, A.R. Proteomic Analysis Shows Constitutive Secretion of MIF and p53-associated Activity of COX-2−/− Lung Fibroblasts. Genomics. Proteom. Bioinform. 2017, 15, 339–351. [Google Scholar] [CrossRef]
  70. Luanpitpong, S.; Chanvorachote, P. Nitric Oxide and Aggressive Behavior of Lung Cancer Cells. Anticancer Res. 2015, 35, 4585–4592. [Google Scholar]
  71. Pardo, O.E.; Wellbrock, C.; Khanzada, U.K.; Aubert, M.; Arozarena, I.; Davidson, S.; Bowen, F.; Parker, P.J.; Filonenko, V.V.; Gout, I.T.; et al. FGF-2 protects small cell lung cancer cells from apoptosis through a complex involving PKCɛ, B-Raf and S6K2. EMBO J. 2006, 25, 3078–3088. [Google Scholar] [CrossRef]
  72. Kwabi-Addo, B.; Ozen, M.; Ittmann, M. The role of fibroblast growth factors and their receptors in prostate cancer. Endocr. Relat. Cancer 2004, 11, 709–724. [Google Scholar] [CrossRef] [Green Version]
  73. Mantovani, A.; Barajon, I.; Garlanda, C. IL-1 and IL-1 regulatory pathways in cancer progression and therapy. Immunol. Rev. 2018, 281, 57–61. [Google Scholar] [CrossRef]
  74. Zheng, X.; Hu, Y.; Yao, C. The paradoxical role of tumor-infiltrating immune cells in lung cancer. Intractable Rare Dis. Res. 2017, 6, 234–241. [Google Scholar] [CrossRef] [Green Version]
  75. Choi, H.; Na, K.J. Integrative analysis of imaging and transcriptomic data of the immune landscape associated with tumor metabolism in lung adenocarcinoma: Clinical and prognostic implications. Theranostics 2018, 8, 1956–1965. [Google Scholar] [CrossRef] [Green Version]
  76. Ohue, Y.; Kurose, K.; Nozawa, R.; Isobe, M.; Nishio, Y.; Tanaka, T.; Doki, Y.; Hori, T.; Fukuoka, J.; Oka, M.; et al. Survival of Lung Adenocarcinoma Patients Predicted from Expression of PD-L1, Galectin-9, and XAGE1 (GAGED2a) on Tumor Cells and Tumor-Infiltrating T Cells. Cancer Immunol. Res. 2016, 4, 1049–1060. [Google Scholar] [CrossRef] [Green Version]
  77. Kanda, M.; Tanaka, C.; Kobayashi, D.; Tanaka, H.; Shimizu, D.; Shibata, M.; Takami, H.; Hayashi, M.; Iwata, N.; Niwa, Y.; et al. Epigenetic suppression of the immunoregulator MZB1 is associated with the malignant phenotype of gastric cancer. Int. J. Cancer 2016, 139, 2290–2298. [Google Scholar] [CrossRef] [Green Version]
  78. Zhai, Y.; Chen, Y.; Li, Q.; Zhang, L. Exploration of the hub genes and miRNAs in lung adenocarcinoma. Oncol. Lett. 2019, 18, 1713–1722. [Google Scholar] [CrossRef] [Green Version]
  79. Wei, S.C.; Duffy, C.R.; Allison, J.P. Fundamental Mechanisms of Immune Checkpoint Blockade Therapy. Cancer Discov. 2018, 8, 1069–1086. [Google Scholar] [CrossRef] [Green Version]
  80. Xiong, D.; Wang, Y.; You, M. Tumor intrinsic immunity related proteins may be novel tumor suppressors in some types of cancer. Sci. Rep. 2019, 9, e10918. [Google Scholar] [CrossRef]
  81. Kanteti, R.; Nallasura, V.; Loganathan, S.; Tretiakova, M.; Kroll, T.; Krishnaswamy, S.; Faoro, L.; Cagle, P.; Husain, A.N.; Vokes, E.E.; et al. PAX5 is expressed in small-cell lung cancer and positively regulates c-Met transcription. Lab. Investig. 2009, 89, 301–314. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  82. O’Brien, P.; Morin, P.; Ouellette, R.J.; Robichaud, G.A. The Pax-5 Gene: A Pluripotent Regulator of B-cell Differentiation and Cancer Disease. Cancer Res. 2011, 71, 7345–7350. [Google Scholar] [CrossRef] [Green Version]
  83. Zhao, L.; Li, S.; Gan, L.; Li, C.; Qiu, Z.; Feng, Y.; Li, J.; Li, L.; Li, C.; Peng, W.; et al. Paired box 5 is a frequently methylated lung cancer tumour suppressor gene interfering β-catenin signalling and GADD45G expression. J. Cell. Mol. Med. 2016, 20, 842–854. [Google Scholar] [CrossRef]
  84. Pai, S.G.; Carneiro, B.A.; Mota, J.M.; Costa, R.; Leite, C.A.; Barroso-Sousa, R.; Kaplan, J.B.; Chae, Y.K.; Giles, F.J. Wnt/beta-catenin pathway: Modulating anticancer immune response. J. Hematol. Oncol. 2017, 10, 101. [Google Scholar] [CrossRef] [Green Version]
  85. Liu, J.; Li, R.; Liao, X.; Jiang, W. Comprehensive Bioinformatic Analysis Genes Associated to the Prognosis of Liposarcoma. Med. Sci. Monit. 2018, 24, 7329–7339. [Google Scholar] [CrossRef]
  86. Koksal, D.; Demirag, F.; Bayiz, H.; Ozmen, O.; Tatci, E.; Berktas, B.; Aydoğdu, K.; Yekeler, E. The correlation of SUVmax with pathological characteristics of primary tumor and the value of Tumor/ Lymph node SUVmax ratio for predicting metastasis to lymph nodes in resected NSCLC patients. J. Cardiothorac. Surg. 2013, 8, e63. [Google Scholar] [CrossRef] [Green Version]
  87. Qian, X.; Song, X.; He, Y.; Yang, Z.; Sun, T.; Wang, J.; Zhu, G.; Xing, W.; You, C. CCNB2 overexpression is a poor prognostic biomarker in Chinese NSCLC patients. Biomed. Pharmacother. 2015, 74, 222–227. [Google Scholar] [CrossRef]
  88. Wang, T.; Chen, B.; Yang, Y.; Chen, H.; Wang, Y.; Cviko, A.; Quade, B.J.; Sun, D.; Yang, A.; McKeon, F.D.; et al. Histologic and immunophenotypic classification of cervical carcinomas by expression of the p53 homologue p63: A study of 250 cases. Hum. Pathol. 2001, 32, 479–486. [Google Scholar] [CrossRef]
  89. Shi, Y.; Li, Y.; Yan, C.; Su, H.; Ying, K. Identification of key genes and evaluation of clinical outcomes in lung squamous cell carcinoma using integrated bioinformatics analysis. Oncol. Lett. 2019, 18, 5859–5870. [Google Scholar] [CrossRef] [Green Version]
  90. Gao, M.; Kong, W.; Huang, Z.; Xie, Z. Identification of key genes related to the mechanism and prognosis of lung squamous cell carcinoma using bioinformatics analysis. arXiv 2019, 1911, 05621. [Google Scholar]
  91. Hsu, Y.-L.; Hung, J.-Y.; Lee, Y.-L.; Chen, F.-W.; Chang, K.-F.; Chang, W.-A.; Tsai, Y.-M.; Chong, I.-W.; Kuo, P.-L. Identification of novel gene expression signature in lung adenocarcinoma by using next-generation sequencing data and bioinformatics analysis. Oncotarget 2017, 8, 104831–104854. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  92. Shi, K.; Li, N.; Yang, M.; Li, W. Identification of Key Genes and Pathways in Female Lung Cancer Patients Who Never Smoked by a Bioinformatics Analysis. J. Cancer 2019, 10, 51–60. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  93. Dobin, A.; Davis, C.A.; Schlesinger, F.; Drenkow, J.; Zaleski, C.; Jha, S.; Batut, P.; Chaisson, M.; Gingeras, T.R. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29, 15–21. [Google Scholar] [CrossRef] [PubMed]
  94. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, e323. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  95. Okonechnikov, K.; Conesa, A.; García-Alcalde, F. Qualimap 2: Advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics 2015, 32, btv566. [Google Scholar] [CrossRef]
  96. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R. 1000 Genome Project Data Processing Subgroup The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009, 25, 2078–2079. [Google Scholar] [CrossRef] [Green Version]
  97. Wang, L.; Wang, S.; Li, W. RSeQC: Quality control of RNA-seq experiments. Bioinformatics 2012, 28, 2184–2185. [Google Scholar] [CrossRef] [Green Version]
  98. Robinson, J.T.; Thorvaldsdóttir, H.; Winckler, W.; Guttman, M.; Lander, E.S.; Getz, G.; Mesirov, J.P. Integrative genomics viewer. Nat. Biotechnol. 2011, 29, 24–26. [Google Scholar] [CrossRef] [Green Version]
Cancers 12 00037 g001 550
Figure 1. Differential gene expression in squamous cell carcinoma (SCC) and adenocarcinoma (ADC). (A) Volcano plot of differentially expressed genes (DEGs) in each subtype of non-small-cell lung carcinoma (NSCLC). The horizontal line at false discovery rate (FDR) = 0.05; vertical line at |log2FC| = 1; (B) Pie charts represent the numbers of genes found to be upregulated or downregulated. (C) Venn diagram shows overlap between differentially expressed genes in SCC and ADC. The numbers of DEGs include both up- and down-regulated genes in tumour tissue compared to normal tissue.
Figure 1. Differential gene expression in squamous cell carcinoma (SCC) and adenocarcinoma (ADC). (A) Volcano plot of differentially expressed genes (DEGs) in each subtype of non-small-cell lung carcinoma (NSCLC). The horizontal line at false discovery rate (FDR) = 0.05; vertical line at |log2FC| = 1; (B) Pie charts represent the numbers of genes found to be upregulated or downregulated. (C) Venn diagram shows overlap between differentially expressed genes in SCC and ADC. The numbers of DEGs include both up- and down-regulated genes in tumour tissue compared to normal tissue.
Cancers 12 00037 g001
Cancers 12 00037 g002 550
Figure 2. Top 20 genes up- and downregulated in SCC and ADC tissue compared to normal tissue (FDR ≤ 0.01). In the heat map colour indicated log fold change, shades of blue represent log2FC < −1, and shades of red represent log2FC > 1 and missing values are in white. FC, fold change.
Figure 2. Top 20 genes up- and downregulated in SCC and ADC tissue compared to normal tissue (FDR ≤ 0.01). In the heat map colour indicated log fold change, shades of blue represent log2FC < −1, and shades of red represent log2FC > 1 and missing values are in white. FC, fold change.
Cancers 12 00037 g002
Cancers 12 00037 g003 550
Figure 3. Gene ontology (GO)-terms analysis. (A) Venn diagram showing the number of GO-biological process (BP) differentially expressed (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and adenocarcinoma (ADC) and the number of overlapping GO-terms between this both subtypes of NSCLC. (B) Top 10 GO-terms identified only in SCC vs. normal tissue. (C) Top 10 GO-terms identified only in ADC vs. normal tissue. (D) The network of enriched terms coloured by cluster identified only in SCC vs. normal tissue using Metascape tool. (E) The network of enriched terms coloured by cluster identified only in ADC vs. normal tissue using Metascape tool.
Figure 3. Gene ontology (GO)-terms analysis. (A) Venn diagram showing the number of GO-biological process (BP) differentially expressed (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and adenocarcinoma (ADC) and the number of overlapping GO-terms between this both subtypes of NSCLC. (B) Top 10 GO-terms identified only in SCC vs. normal tissue. (C) Top 10 GO-terms identified only in ADC vs. normal tissue. (D) The network of enriched terms coloured by cluster identified only in SCC vs. normal tissue using Metascape tool. (E) The network of enriched terms coloured by cluster identified only in ADC vs. normal tissue using Metascape tool.
Cancers 12 00037 g003
Cancers 12 00037 g004 550
Figure 4. The most representative altered canonical pathways associated with SCC and ADC identified using Ingenuity Pathway Analysis (IPA). (A) Venn diagram showing the number of canonical pathways (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and adenocarcinoma (ADC), and the overlapping of those pathways between these both subtypes of NSCLC. (B) Top 10 canonical pathways identified only in SCC vs. normal tissue. (C) Top 10 canonical pathways identified only in ADC vs. normal tissue. ‘Ratio’ refers to the number of molecules from the dataset that map to the pathway listed, divided by the total number of molecules that define the canonical pathway from within the IPA knowledge base (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/).
Figure 4. The most representative altered canonical pathways associated with SCC and ADC identified using Ingenuity Pathway Analysis (IPA). (A) Venn diagram showing the number of canonical pathways (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and adenocarcinoma (ADC), and the overlapping of those pathways between these both subtypes of NSCLC. (B) Top 10 canonical pathways identified only in SCC vs. normal tissue. (C) Top 10 canonical pathways identified only in ADC vs. normal tissue. ‘Ratio’ refers to the number of molecules from the dataset that map to the pathway listed, divided by the total number of molecules that define the canonical pathway from within the IPA knowledge base (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/).
Cancers 12 00037 g004
Cancers 12 00037 g005 550
Figure 5. IPA-predicted top 20 upstream regulators of gene expression in SCC and ADC. (A) Venn diagram showing the number of upstream regulators (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and the number of overlapping regulators between these both subtypes of NSCLC. (B) The heat map shows the predicted activation and inhibition states of categorised molecules in IPA as transcription regulators, growth factors, cytokines, kinases and others. Colours in heat map indicate z-scores.
Figure 5. IPA-predicted top 20 upstream regulators of gene expression in SCC and ADC. (A) Venn diagram showing the number of upstream regulators (FDR ≤ 0.01) in squamous cell lung carcinoma (SCC) and the number of overlapping regulators between these both subtypes of NSCLC. (B) The heat map shows the predicted activation and inhibition states of categorised molecules in IPA as transcription regulators, growth factors, cytokines, kinases and others. Colours in heat map indicate z-scores.
Cancers 12 00037 g005
Cancers 12 00037 g006 550
Figure 6. Cluster analysis of samples of (A) ADC and (B) SCC to detect outliers (white-to-red linear gradient colour associated with corresponding clinical variable, grey when missing data). (C) Determination of soft-thresholding power in weighted gene co-expression network analysis (WGCNA).
Figure 6. Cluster analysis of samples of (A) ADC and (B) SCC to detect outliers (white-to-red linear gradient colour associated with corresponding clinical variable, grey when missing data). (C) Determination of soft-thresholding power in weighted gene co-expression network analysis (WGCNA).
Cancers 12 00037 g006
Cancers 12 00037 g007 550
Figure 7. Hierarchical clustering dendrograms of identified co-expressed genes in modules in (A) ADC and (B) SCC lung cancer. Each coloured row represents a colour-coded module which contains a group of highly connected genes. A total of 32 and 15 modules was identified in ADC and SCC, respectively.
Figure 7. Hierarchical clustering dendrograms of identified co-expressed genes in modules in (A) ADC and (B) SCC lung cancer. Each coloured row represents a colour-coded module which contains a group of highly connected genes. A total of 32 and 15 modules was identified in ADC and SCC, respectively.
Cancers 12 00037 g007
Cancers 12 00037 g008 550
Figure 8. The percentage of genes allocated to a cluster based on TCGA (The Cancer Genome Atlas) data (rows) found in cluster build on MOBIT data (columns).
Figure 8. The percentage of genes allocated to a cluster based on TCGA (The Cancer Genome Atlas) data (rows) found in cluster build on MOBIT data (columns).
Cancers 12 00037 g008
Cancers 12 00037 g009 550
Figure 9. Heatmaps of the correlation between Eigengene and clinical traits of (A) SCC and (B) ADC. Each row corresponds to a module eigengene, and each column corresponds to a clinical characteristic. Each cell contains the corresponding correlation and p-value.
Figure 9. Heatmaps of the correlation between Eigengene and clinical traits of (A) SCC and (B) ADC. Each row corresponds to a module eigengene, and each column corresponds to a clinical characteristic. Each cell contains the corresponding correlation and p-value.
Cancers 12 00037 g009
Cancers 12 00037 g010 550
Figure 10. Top five enrichment results of (A) GO-Biological Process analysis by BINGO and (B) canonical pathway analysis by IPA for red, lightcyan, darkorange, and yellow modules in ADC, as well as for blue and turquoise modules in SCC.
Figure 10. Top five enrichment results of (A) GO-Biological Process analysis by BINGO and (B) canonical pathway analysis by IPA for red, lightcyan, darkorange, and yellow modules in ADC, as well as for blue and turquoise modules in SCC.
Cancers 12 00037 g010
Cancers 12 00037 g011 550
Figure 11. The networks of the top 10 hub genes in modules. Nodes represent genes and red colour indicates the highest rank.
Figure 11. The networks of the top 10 hub genes in modules. Nodes represent genes and red colour indicates the highest rank.
Cancers 12 00037 g011
Table
Table 1. Patient baseline characteristics.
Table 1. Patient baseline characteristics.
CharacteristicNSCLC patients (n = 114)
Age, year (mean ± SD)65.4 ± 7.9
Male, n (%)80 (70.2)
Smoking, PY (mean ± SD)39.6 ± 21.0
BMI (mean ± SD)26.9 ± 4.7
SUV (max) (mean ± SD)8.8 ± 4.6
Histology
Adenocarcinoma (AC), n (%)45 (39.5)
Squamous Cell Carcinoma (SCC), n (%)69 (60.5)
TNM stage, n (%)
I39 (34.2)
II52 (45.6)
III17 (14.9)
IV6 (5.3)
Recurrence (%)15.2
Deaths (%)13.3
SD, standard deviation; PY, pack year; TNM, tumour node metastasis; BMI, body mass index; SUV, standardized uptake value; NSCLC, non-small cell lung cancer

Share and Cite

MDPI and ACS Style

Niemira, M.; Collin, F.; Szalkowska, A.; Bielska, A.; Chwialkowska, K.; Reszec, J.; Niklinski, J.; Kwasniewski, M.; Kretowski, A. Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers 2020, 12, 37. https://doi.org/10.3390/cancers12010037

AMA Style

Niemira M, Collin F, Szalkowska A, Bielska A, Chwialkowska K, Reszec J, Niklinski J, Kwasniewski M, Kretowski A. Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers. 2020; 12(1):37. https://doi.org/10.3390/cancers12010037

Chicago/Turabian Style

Niemira, Magdalena, Francois Collin, Anna Szalkowska, Agnieszka Bielska, Karolina Chwialkowska, Joanna Reszec, Jacek Niklinski, Miroslaw Kwasniewski, and Adam Kretowski. 2020. "Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA)" Cancers 12, no. 1: 37. https://doi.org/10.3390/cancers12010037

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop