Abstract
Background:
With the advent of large-scale molecular profiling, an increasing number of oncogenic drivers contributing to precise medicine and reshaping classification of lung adenocarcinoma (LUAD) have been identified. However, only a minority of patients archived improved outcome under current standard therapies because of the dynamic mutational spectrum, which required expanding susceptible gene libraries. Accumulating evidence has witnessed that understanding gene regulatory networks as well as their changing processes was helpful in identifying core genes which acted as master regulators during carcinogenesis. The present study aimed at identifying key genes with differential correlations between normal and tumor status.
Methods:
Weighted gene co-expression network analysis (WGCNA) was employed to build a gene interaction network using the expression profile of LUAD from The Cancer Genome Atlas (TCGA). R package DiffCorr was implemented for the identification of differential correlations between tumor and adjacent normal tissues. STRING and Cytoscape were used for the construction and visualization of biological networks.
Results:
A total of 176 modules were detected in the network, among which yellow and medium orchid modules showed the most significant associations with LUAD. Then genes in these two modules were further chosen to evaluate their differential correlations. Finally, dozens of novel genes with opposite correlations including ATP13A4-AS1, HIGD1B, DAP3, and ISG20L2 were identified. Further biological and survival analyses highlighted their potential values in the diagnosis and treatment of LUAD. Moreover, real-time qPCR confirmed the expression patterns of ATP13A4-AS1, HIGD1B, DAP3, and ISG20L2 in LUAD tissues and cell lines.
Conclusion:
Our study provided new insights into the gene regulatory mechanisms during transition from normal to tumor, pioneering a network-based algorithm in the application of tumor etiology.
Introduction
Lung cancer has been one of the most leading causes responsible for cancer mortality globally for several decades (). According to the Surveillance, Epidemiology, and End Results (SEER) Program, 228,820 new cases and 135,720 deaths of lung and bronchus cancer patients have been estimated in the United States in 2020 (). In China, estimated new cases and deaths from lung cancer were 730,000 and 610,000 in 2015, respectively, accounting for over 30% of the world’s total (). Based on histological subtypes, approximately 85% are of non-small cell lung cancer (NSCLC), of which lung adenocarcinoma (LUAD) is the most prevalent and attributed to more than 50% of all lung cancer cases (). Despite that conventional therapies including surgery, radiotherapy, chemotherapy, immunotherapy, and targeted drugs have been applied to LUAD treatment, the 5-year survival rate for patients remains less than 15%, mainly due to lack of early detection and intervention (). In addition to environmental exposures and tobacco smoking, genetic susceptibility has been recognized as the most important risk factor associated with LUAD. Recently, a new conception called “oncogene addiction” referring to dependence of cancer cells on the activation of specific oncogenes has attracted increasing attention and underlined the significance of oncogenic drivers in LUAD (). At present, patients with LUAD harboring gene aberrations in EGFR (), KRAS (), ALK (), BRAF (), HER2 (), and MET () benefit the most from targeted therapies, and those genes have been considered as oncogenic drivers in the field of thoracic oncology. However, on account of the dynamic mutational spectrum comprised of a vast number of hidden drivers, only a minority of patients archived improved outcome under current standard therapies (). Thus, it is imperative to identify new risk genes for elucidating lung carcinogenetic mechanisms in order to guide researchers to develop new therapeutic strategies and physicians to tailor the treatment options.
Exponential advances in the high-throughput sequencing technologies and informatics have generated large-scale omics data which promotes a paradigm shift in the study of biomedical sciences and is available for deciphering molecular characteristics of oncogenesis that could be translated into clinical practice (; ). Among them, transcriptome data is quite informative because of the discovery of the significantly altered abundance of cellular components or pathways between disease and healthy tissues or two disease states (). In recent years, comparative analysis of transcriptome data has helped to discover prognostic markers and signatures by looking for differentially expressed genes in a variety of cancer types (; ). However, a fact that could not be ignored is that most of the biological activities require an orchestrated action of multiple genes, whose dysregulation could lead to the occurrence of cancer. Therefore, gene correlation approaches have been intensively used for transcriptional profiling, providing preliminary steps toward genetic interaction networks and offering clues about the function of unknown genes.
Complementary to gene correlation analysis, changes in correlation patterns under different conditions, referred to as “differential correlations,” are also attractive as for the reconstruction of genome architecture and identification of regulators or marker genes (; ). For example, one transcription factor defined as a master gene regulates a number of downstream targets. Also, these targets express in an ordered module where the regulatory mechanism is functional. However, in disease tissues where the regulatory mechanism is malfunctional, the gene expression module may be disordered or random. In this case, correlation changes can be detected by differential correlations rather than differential expression analysis. Recently, such alterations in network structures and variance in gene expression levels have been observed in cellular differentiation and cancer cell formation (; ), contributing to understanding the expression patterns between two states. However, differential correlations in LUAD are poorly determined, so it is imperative to unravel network dynamics that can be used to identify new candidate genes. Hitherto, an increasing number of systematic approaches have been developed for identifying such changes in network structures (; ). Application of these methods to LUAD and adjacent normal tissues may reveal variance in gene expression levels during carcinogenesis.
Weighted gene co-expression network analysis (WGCNA) is a systematic method widely utilized in oncology research that aims at finding co-expressed genes through calculating gene connectivity (; ). The soft thresholding of the Pearson correlation matrix is integrated into analysis for determining the connection strengths within gene pairs (; ). Based on Fisher’s z-test, R package DiffCorr is a useful tool in identifying differential correlations between disease states and providing a list of differentially correlated gene pairs (). In this study, we developed an in silico framework to identify key genes with differential correlations (Figure 1). First, we employed WGCNA to build a gene interaction network using the expression profile of LUAD from The Cancer Genome Atlas (TCGA). A total of 176 modules were detected in the network, among which yellow and medium orchid modules showed the most significant associations with LUAD. Then, we calculated differential correlations of genes in these two modules and identified significant differences between tumor and adjacent normal tissues using DiffCorr. Finally, dozens of new genes including ATP13A4-AS1, HIGD1B, DAP3, and ISG20L2, were identified and their expression patterns were confirmed by real-time qPCR. Further biological and survival analysis yielded their valuable effects on the progression of LUAD.
FIGURE 1
Materials and Methods
LUAD RNA-Sequencing Datasets
The RNA-sequencing data of LUAD was downloaded from the TCGA database1, including 468 tumor samples and 58 normal samples. As previously described, the gene expression levels were quantified as FPKM (fragments per kilobase per million mapped reads) using TopHat and HTSeq-count (; ). The TCGA sample information was listed in Supplementary Table 1. Within the 468 tumor samples, 3 samples did not have stage information, 165 samples were in stage I, 247 samples were in stage II, 39 samples were in stage III, and 14 samples were in stage IV. The sample distribution also stressed the importance of identification of new genes for early diagnosis.
Co-expression Network Analysis
Hub gene screening and co-expression of gene pair detection were performed by an R package WGCNA, which was designed for analyzing a weighted correlation network (). The premise of the network construction was that elements in the gene co-expression matrix were the weighted values of the correlation coefficient between gene pairs. The selection criterion of the weight was to make the connections among genes conform to the scale-free network distribution; that is, for the number of connections i, the probability p(i) was inverse to in. Practically, a weighted coefficient was selected to approximate the scale-free topology, which needs to satisfy the following condition: the log of the number of connected nodes log(i) was negatively correlated with the log of occurrence probability of nodes log(p(i)) (the correlation coefficient should be at least 0.8). At the same time, the average connection degree of genes in different modules should be quite high. Associated genes were clustered based on dissimilarity of the unsigned topological overlap matrix (TOM). Finally, network modules and the genes within them were identified.
Functional Enrichment Analysis of Module Genes
R package clusterProfiler was used to perform functional enrichment analysis on clustered genes in yellow and medium orchid modules. A hypergeometric distribution test was applied to detect enrichment terms, and P values were adjusted by false discovery rate (FDR) method with a cutoff FDR < 0.05 ().
Differential Correlation Analysis
R package DiffCorr was implemented for the identification and visualization of differential correlations in biological networks, and details were described in . Briefly, the DiffCorr package mainly contained three functions. First is calculation of differential correlations. The correlation coefficients for each of the two conditions (herein referred to tumor and normal), rA and rB, were transformed into ZA and ZB, respectively, based on Fisher’s transformation: . Differences between the two correlations could be tested using the equation , where nA and nB represent the sample size for each biomolecule pair in each condition. Then, the local false discovery rate (lfdr) derived from the fdrtool package was used for controlling true estimates and identifying significant differential correlations. Second is identifying eigen-molecules. Eigen-molecules or “eigengenes” in the network were calculated based on the first principal component of a data matrix of a module which was extracted from a hierarchical cluster analysis. In addition to pair-wise differential correlations between molecules, differential correlations between modules were tested by using these eigen-molecule modules. Third is scaling and clustering. Different pretreatment methods with downstream correlation analyses were integrated, including auto-scaling, range scaling, Pareto scaling, vast scaling, level scaling, and power transformation.
Visualization and Functional Analysis
Construction and analysis of networks were carried out using STRING (11.0)2 (). Visualization of networks was realized using Cytoscape (3.6.0) ().
The Survival Analysis of Key Genes With Differential Correlations
There were 468 LUAD patients with survival information. We performed survival analysis using the Cox proportional hazard regression model on these samples (). Each gene expression level was divided into two groups and repeated for 90 times based on the cutoff value from 5 to 95% of its expression level. A repeated log-rank test based on each cutoff value was processed, and a cutoff value with the lowest P-value was selected for subsequent univariate survival analysis. The Kaplan–Meier plot was used to describe the survival curves of these two groups of patients. The significance of the survival difference between these two patient groups was evaluated by the log-rank test P value. If the P value was less than 0.05, its survival was considered as significantly different. The R package survival was used to perform the survival analysis.
Cell Lines and Cell Culture
The LUAD cell lines (A549, SPCA1, PC9, H1299, and H1975) and the normal human bronchial epithelial cell line 16HBE were from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences, Shanghai Institute of Cell Biology (Shanghai, China). A549, H1975, and H1299 cells were grown in RPMI-1640 (Thermo, United States) containing 10% FBS, and 16HBE, SPCA1, and PC9 cells were grown in a DMEM medium (Thermo, United States) containing 10% fetal bovine serum (FBS). In addition, all cell media were supplemented with penicillin (100 U/mL) and streptomycin (100 U/mL) at 37°C in a 5% CO2 incubator.
LUAD Patients and Tissue Specimens
A total of 15 pairs of LUAD tumor and corresponding adjacent normal tissues were collected from the Third Affiliated Hospital of Soochow University between October 2019 and June 2020. Samples were snap-frozen and stored at −80°C until use in real-time qPCR (RT-qPCR) experiments. In addition, we conducted immunohistochemical staining of formalin-fixed paraffin-embedded LUAD patients and normal control samples. The Research Ethics Committee of the Third Affiliated Hospital of Soochow University approved this study, which was consistent with the Declaration of Helsinki. All patients provided written informed consent.
RNA Isolation, cDNA Synthesis, and RT-qPCR
We isolated total RNA with TRIzol reagent (Thermo, United States). First-strand cDNA was synthesized from 1 μg total RNA using a ReverTra Ace qPCR RT Kit (Toyobo, Osaka, Japan). The Fast SYBR Green Master Mix was used for RT-qPCR (Applied Biosystems Inc., CA, United States). The cycling conditions were 30 s of polymerase activation at 95°C followed by 40 cycles at 95°C for 5 s and 60°C for 30 s. GAPDH was used as an internal loading control. The relative level was calculated by the relative quantification 2-ΔΔCT method. All primer sequences were listed in Supplementary Table 2.
Statistical Analysis
Student’s t-test was applied to identify genes differentially expressed between normal and tumor samples. P values were adjusted by the Benjamini–Hochberg method (). Differentially expressed genes were defined as adjusted P value < 0.05. Fisher’s z-test was employed to evaluate differential correlations of gene pairs between normal and tumor samples. Moreover, lfdr < 0.05 was defined as significant differential correlations.
Results
Construction of the Co-expression Network
The weighted gene co-expression network was constructed from 58,387 coding and noncoding RNAs through the WGCNA approach. Here, the soft-thresholding power was set to be six to satisfy the scale-free topology of the network (Figure 2), in which R2 was used to check how well the network fit the scale freeness. When the soft-thresholding power was set to be six, the R2 was 0.91. Furthermore, we detected 176 modules in this network, whose relationship was shown in a cluster dendrogram (Figure 3A). The number of members in different modules varied widely. The members of each module were listed in Supplementary Table 3. Besides the gray module which comprised many unclassified members, the turquoise module contained a maximum of 4,594 genes, while a minimum of 30 genes were included in the dark sea green module.
FIGURE 2
FIGURE 3
Each module represented a group of genes with similar expression profiles across samples. Next, we quantified module-trait associations (Supplementary Figure 1), among which the yellow and medium orchid modules showed the most significant associations with LUAD. The corresponding correlation coefficients of yellow and medium orchid modules were 0.83 (P = 410–133) and −0.58 (P = 510–49), respectively. Clearly, Gene Significance (GS) and Module Membership (MM) analysis illustrated that genes highly significantly associated with LUAD were also the most important elements of modules associated with LUAD (Supplementary Figure 2). We also performed Gene Ontology (GO) enrichment analysis of these two modules. As presented in Figure 3B, genes in the yellow module were significantly enriched in the negative regulation of growth and vasculogenesis, while genes in the medium orchid module were enriched in the cellular metabolic process (geranyl diphosphate metabolic process, geranyl diphosphate biosynthetic process, farnesyl diphosphate biosynthetic process, and farnesyl diphosphate metabolic process) with an adjusted P value smaller than 0.05, conferring the importance of these biological functions on LUAD development.
Identification of Differential Correlations
The genes in the yellow and medium orchid modules were further chosen to evaluate their differential correlations. Using R package DiffCorr, these genes were grouped based on their expression patterns in each subtype (normal or tumor) using the cluster.molecule function. We used the one-correlation coefficient as a distance measure (the cutoff of the coefficient was 0.6) according to the cutree function. The get.eigen.molecule and get.eigen.molecule.graph functions were used for visualization of the module network (Figure 4). The comp.2.cc.fdr function provided the resulting pair-wise differential correlations in the yellow (Supplementary Table 4) and medium orchid modules (Supplementary Table 5).
FIGURE 4
The DiffCorr package also detected oppositely correlated pairs where, for example, two genes exhibited a positive correlation in normal samples and a negative correlation in tumor samples, or vice versa, a condition referred to as a “switching mechanism” (). These switched gene pairs, which especially show a differential expression at the same time, were worth noticing for their critical roles in understanding cellular functions in the development of LUAD. In total, we obtained 42 oppositely correlated gene pairs with a differential expression simultaneously from the yellow module and 62 from the medium orchid module (Supplementary Table 6). Their interaction networks are shown in Figure 5. The top 10 significant switching mechanisms of gene expression between normal and tumor samples from the yellow and medium orchid modules are shown in Table 1.
FIGURE 5
TABLE 1
| Molecule X | Molecule Y | r1 (normal) | r2 (tumor) | lfdra (difference) | Module color |
| HIGD1B | DUOX1 | –0.58 | 0.53 | 0 | Yellow |
| FAM162B | DUOX1 | –0.56 | 0.51 | 0 | Yellow |
| EPCAM | CAV1 | 0.65 | –0.40 | 0 | Yellow |
| FAM162B | DUOXA1 | –0.59 | 0.46 | 4.67E-12 | Yellow |
| PHACTR1 | HIGD1B | –0.58 | 0.45 | 1.47E-11 | Yellow |
| DUOX1 | COX4I2 | –0.55 | 0.47 | 3.16E-11 | Yellow |
| IL33 | DUOX1 | –0.53 | 0.48 | 8.28E-11 | Yellow |
| FXYD1 | ATP13A4-AS1 | –0.58 | 0.41 | 1.29E-10 | Yellow |
| MMP19 | CASS4 | –0.43 | 0.55 | 2.30E-10 | Yellow |
| ATP13A4-AS1 | AGER | –0.40 | 0.58 | 2.51E-10 | Yellow |
| VPS72 | PIP5K1A | –0.47 | 0.73 | 0 | Mediumorchid |
| VPS45 | PIP5K1A | –0.41 | 0.69 | 0 | Mediumorchid |
| UBE2Q1 | ILF2 | –0.60 | 0.58 | 0 | Mediumorchid |
| TARS2 | PIP5K1A | –0.54 | 0.66 | 0 | Mediumorchid |
| TARS2 | ENSA | –0.46 | 0.62 | 0 | Mediumorchid |
| PYGO2 | PIP5K1A | –0.59 | 0.62 | 0 | Mediumorchid |
| PSMD4 | PIP5K1A | –0.55 | 0.60 | 0 | Mediumorchid |
| PRUNE | PIP5K1A | –0.58 | 0.70 | 0 | Mediumorchid |
| PRUNE | CDC42SE1 | –0.55 | 0.59 | 0 | Mediumorchid |
| PIP5K1A | LYSMD1 | –0.74 | 0.71 | 0 | Mediumorchid |
Top 10 correlated gene pairs changed to the opposite direction from the yellow and medium orchid modules between normal and LUAD samples.
aLocal FDR.
Functional Analysis of Key Genes With Differential Correlations
Next, we focused on several key genes that acted as master regulators for their occupying most connections with other genes. In the yellow module, DUOX1, DUOXA1, HIGD1B, and ATP13A4-AS1 stood out in the network (Figure 5A). ATP13A4 antisense RNA 1 (ATP13A4-AS1) is an antisense long noncoding RNA (lncRNA) derived from ATP13A4 with an unknown function in cancer. Here, we can infer the roles of ATP13A4-AS1 from its interacted genes advanced glycosylation end-product specific receptor (AGER) and angiopoietin-like 7 (ANGPTL7) (Figure 5A). As a member of the immunoglobulin superfamily, cell surface receptor AGER is involved in multiple inflammatory responses whose abnormal expression has been reported to be closely associated with carcinogenesis (). Accumulating evidence has shown the downregulation of AGER in lung cancer, leading to enhanced proliferation, invasion and migration abilities, and decreased apoptosis of cells (; ). In addition, genetic polymorphisms of AGER have been reported to increase risks of lung cancer and breast cancer (; ). Angiopoietin-like 7 (ANGPTL7) belongs to a family of secreted angiopoietin-like proteins and plays vital roles in the modulation of hematopoietic stem cell maintenance, angiogenesis, and lipid metabolism (; ). Prior studies have revealed striking ANGPTL7 underexpression in various cancers such as colorectal cancer and breast cancer. Besides, upregulation of ANGPTL7 has been found in cancer cells after myeloid cell depletion, which affected liver metastasis by diminishing cell growth and vascular density, implying that ANGPTL7 could act as a mediator of metastatic progression and as a promising intervention target (). As AGER and ANGPTL7 both show strong links to cancer, it is plausible that their interacted gene ATP13A4-AS1 also engages in the progression of lung cancer by exerting similar impacts on cellular activities. Furthermore, the Kaplan–Meier plot exhibited that lower levels of ATP13A4-AS1 correlated with poor patient outcome (Figure 6A), corresponding to its lower expression in LUAD than normal samples confirmed by our RT-qPCR (Figure 7A). Hence, our study revealed, for the first time, the roles of ATP13A4-AS1 in carcinogenesis, providing experimental clues for further investigations.
FIGURE 6
FIGURE 7
Dual oxidase 1 (DUOX1), an oxidant-generating enzyme within the airway epithelium, has been previously reported to take part in innate airway host defense and epithelial homeostasis, the activation of which stimulated cell migration dependent on Src family tyrosine kinases and epidermal growth factor receptor (EGFR) signaling pathways (; ; ; ). DUOXA1, as a maturation factor of DUOX1, is transcriptionally and functionally linked to DUOX1 (). Recent evidence has indicated that DUOX1 and DUOXA1 were frequently silenced due to promoter hypermethylation in various epithelial cancers including lung cancer (; ; ). Consistent with pan-cancer analysis, DUOX1 and DUOXA1 were also downregulated in LUAD from our TCGA dataset, proving their potential tumor suppression roles.
HIG1 domain family member 1B (HIGD1B) is localized to the cell membrane whose expression is induced by hypoxia and glucose deprivation (). Moreover, an increased expression of HIGD1B has been observed in cervical cancer and plurihormonal pituitary adenomas (; ). In contrast with previous research, HIGD1B showed lower expression levels in LUAD compared with normal samples in TCGA datasets, which was also observed in our LUAD patients and cell lines (Figure 7B), probably due to different cancer characteristics, implying its dual function in carcinogenesis. Moreover, its interaction with DUOX1 and DUOXA1 emphasized similar contributions to inhibit the malignant progress of cancer. The survival analysis also demonstrated that the low expression was associated with high risk (Figure 6B). Thus, our study expanded the roles of HIGD1B in cancer biology which need further validations.
In the medium orchid module, PIP5K1A, CRTC2, DAP3, and ISG20L2 deserved much attention for their central positions in the network (Figure 5B). As a lipid kinase, the basic function of PIP5K1A is to phosphorylate PI4P to synthesize important signaling phospholipid PI(4,5)P2, which serves as the substrate for phosphoinositide 3-kinase (PI3K) for conversion into PI(3,4,5)P3 to promote cell proliferation and survival (). Intriguingly, PIP5K1A has been recently shown to interact directly with mutant KRAS and TP53, which were the most prevalent drivers in lung cancer, and facilitated downstream oncogenic signaling (; ).
CRTC2, belonging to the CREB-regulated transcription coactivator (CRTC) family, can bind the leucine zipper DNA-binding region of CREB which results in the enhancement of CREB transcriptional activity (; ). Recent studies have identified novel mutations and a high expression of CRTC2 in non-small cell lung cancer patients, leading to reinforced migration and invasion abilities of cancer cells (; ). Current evidence has pointed out oncogenic roles of CRTC2 in lung cancer, consistent with its elevated expression levels in our study.
As a molecule involved in controlling apoptosis and anoikis, death-associated protein 3 (DAP3) has been highly implicated in the context of carcinogenesis. However, controversy exists with regard to its roles in human cancer (). Increased expression levels of DAP3 have been observed in invasive glioblastoma as well as glioma cells with induced migratory phenotype (), while an inverse association between DAP3 expression and clinical outcome has been found in breast cancer, corresponding to the pro-apoptotic function of DAP3 (). A recent study has identified differentially expressed DAP3 for classification between LUAD and normal samples (), indicating its crucial roles in carcinogenesis. Consistently, our Kaplan–Meier analysis showed better survival in the low transcription group approaching significance (Figure 6C), supporting tumor promotion influences of DAP3 on LUAD. At the same time, our experiments also witnessed the upregulation of DAP3 in LUAD patients and cell lines (Figure 7C). Therefore, our study illustrated that, for the first time, DAP3 occupied the core position in the regulatory network of LUAD, explaining its mechanism of action to some extent.
Interferon-stimulated 20-kDa exonuclease-like 2 (ISG20L2) is a novel vertebrate nucleolar exoribonuclease and involved in ribosome biogenesis (; ). As one of target genes regulated by miR-139-3p, ISG20L2 has recently been related to hepatocellular carcinoma prognosis (), consistent with an immunogenomic landscape study which adds the roles of ISG20L2 in reflecting levels of infiltration in diverse immune cells (). Additionally, high expression of ISG20L2 has been found in several cancers and affected cancer occurrence in a variety of aspects (; ). Nevertheless, there is scant literature with regard to its role in LUAD. Here, we proposed for the first time the roles of ISG20L2 in promoting carcinogenesis and identified its interacted genes in the regulatory network which contributed to understanding further functional mechanisms. Accordingly, survival time was significantly higher in patients with low ISG20L2 expression, compared with the high ISG20L2 expression group (Figure 6D). Meanwhile, remarkable elevated levels of ISG20L2 were also verified by RT-qPCR in LUAD patients and cell lines (Figure 7D).
Discussion
Although multiple studies have continuously explored gene interaction networks which were constructed from a series of correlated gene pairs, little is known about the differential correlations between normal and tumor status which contributed to elucidating deep mechanisms of gene regulation and identifying master genes. Benefiting from the availability of high-throughput data and different kinds of algorithms, detecting correlation changes became feasible. In this study, we developed an in silico framework to identify key genes with differential correlations in LUAD. Both candidate genes and their targets from gene regulatory networks could be further used for experimental investigations of biological functions in order to guide the diagnosis and treatment of patients. As transcriptomic data only represents a single layer of genome and is not sufficient to depict the whole-cell atlas, multidimensional molecular data including proteomic and metabolomic data as well as genome-wide association studies are required to be integrated for elaborating transition mechanisms from normal to tumor. Meanwhile, a statistical evaluation should be more strongly emphasized in future studies of differential correlations.
Statements
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The studies involving human participants were reviewed and approved by The Third Affiliated Hospital of Soochow University. The patients/participants provided their written informed consent to participate in this study.
Author contributions
YoZ conceived the study. YiZ and JL conducted the database search and analysis. BX, TH, and XZ designed the algorithm and network construction. YL and HD conducted the experimental validations. YoZ and JJ wrote the manuscript. All authors contributed to biological analysis, interpretation of the results, read, and approved the final version of the manuscript.
Funding
This project was supported by grants from the National Natural Science Foundation of China (31701111, 31701151, and 81902386), Strategic Priority Research Program of Chinese Academy of Sciences (XDB38050200), National Key R&D Program of China (2018YFC0910403 and 2017YFC1201200), Shanghai Municipal Science and Technology Major Project (2017SHZDZX01), Young Talent Development Plan of Changzhou Health Commission (CZQM2020008), and Youth Talent Science and Technology Project of Changzhou Health Commission (QN202014).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.675438/full#supplementary-material
Supplementary Figure 1Module-trait associations. Each row corresponded to a module eigengene, column to a trait. Each cell contained the corresponding correlation and P-value. The table was color-coded by correlation according to the color legend.
Supplementary Figure 2A scatterplot of Gene Significance (GS) for Disease vs. Module Membership (MM) in the yellow (A) and mediumorchid (B) modules. There was a highly significant correlation between GS and MM in the module.
Supplementary Table 1The TCGA sample IDs and their clinical information.
Supplementary Table 2Primer sequences of GAPDH, ATP13A4-AS1, HIGD1B, DAP3, and ISG20L2.
Supplementary Table 3The members of each module.
Supplementary Table 4Pair-wise differential correlations in the yellow module.
Supplementary Table 5Pair-wise differential correlations in the mediumorchid module.
Supplementary Table 6Oppositely correlated gene pairs with differentially expression simultaneously from the yellow and mediumorchid modules.
References
1
AdhikariH.CounterC. M. (2018). Interrogating the protein interactomes of RAS isoforms identifies PIP5K1A as a KRAS-specific vulnerability.Nat. Commun.9:3646. 10.1038/s41467-018-05692-6
2
AndersS.PylP. T.HuberW. (2015). HTSeq–a Python framework to work with high-throughput sequencing data.Bioinformatics31166–169. 10.1093/bioinformatics/btu638
3
AndersenP. K.GillR. D. (1982). Cox’s regression model for counting processes: a large sample study.Ann. Statist.101100–1120.
4
AndoT.KatoR.HondaH. (2015). Differential variability and correlation of gene expression identifies key genes involved in neuronal differentiation.BMC Syst. Biol.9:82. 10.1186/s12918-015-0231-6
5
BockmayrM.KlauschenF.GyörffyB.DenkertC.BudcziesJ. (2013). New network topology approaches reveal differential correlation patterns in breast cancer.BMC Syst. Biol.7:78. 10.1186/1752-0509-7-78
6
BongarzoneS.SavickasV.LuziF.GeeA. D. (2017). Targeting the receptor for advanced glycation endproducts (RAGE): a medicinal chemistry perspective.J. Med. Chem.607213–7232. 10.1021/acs.jmedchem.7b00058
7
CalvayracO.PradinesA.PonsE.MazieresJ.GuibertN. (2017). Molecular biomarkers for lung adenocarcinoma.Eur. Respir. J.49:1601734. 10.1183/13993003.01734-2016
8
ChangX.ShiL.GaoF.RussinJ.ZengL.HeS.et al (2013). Genomic and transcriptome analysis revealing an oncogenic functional module in meningiomas.Neurosurg. Focus35:E3. 10.3171/2013.10.focus13326
9
ChenW.ZhengR.BaadeP. D.ZhangS.ZengH.BrayF.et al (2016). Cancer statistics in China, 2015.CA Cancer J. Clin.66115–132. 10.3322/caac.21338
10
ChengL. (2020). Omics data and artificial intelligence: new challenges for gene therapy.Curr. Gene Ther.20:1. 10.2174/156652322001200604150041
11
ChoiS.ChenM.CrynsV. L.AndersonR. A. (2019). A nuclear phosphoinositide kinase complex regulates p53.Nat. Cell Biol.21462–475. 10.1038/s41556-019-0297-2
12
CoutéY.KindbeiterK.BelinS.DieckmannR.DuretL.BezinL.et al (2008). ISG20L2, a novel vertebrate nucleolar exoribonuclease involved in ribosome biogenesis.Mol. Cell. Proteomics7546–559. 10.1074/mcp.M700510-MCP200
13
DanknerM.RoseA. A. N.RajkumarS.SiegelP. M.WatsonI. R. (2018). Classifying BRAF alterations in cancer: new rational therapeutic strategies for actionable mutations.Oncogene373183–3199. 10.1038/s41388-018-0171-x
14
DenkoN.SchindlerC.KoongA.LaderouteK.GreenC.GiacciaA. (2000). Epigenetic regulation of gene expression in cervical cancer cells by the tumor microenvironment.Clin. Cancer Res.6480–487.
15
DentinR.HedrickS.XieJ.YatesJ.IIIMontminyM. (2008). Hepatic glucose sensing via the CREB coactivator CRTC2.Science3191402–1405. 10.1126/science.1151363
16
DonchevaN. T.MorrisJ. H.GorodkinJ.JensenL. J. (2019). Cytoscape StringApp: network analysis and visualization of proteomics data.J. Proteome Res.18623–632. 10.1021/acs.jproteome.8b00702
17
DrilonA.CappuzzoF.OuS. I.CamidgeD. R. (2017). Targeting MET in lung cancer: will expectations finally be MET?J. Thorac. Oncol.1215–26. 10.1016/j.jtho.2016.10.014
18
FukushimaA. (2013). DiffCorr: an R package to analyze and visualize differential correlations in biological networks.Gene518209–214. 10.1016/j.gene.2012.11.028
19
GhazalpourA.DossS.ZhangB.WangS.PlaisierC.CastellanosR.et al (2006). Integrating genetic and network analysis to characterize genes related to mouse weight.PLoS Genet.2:e130. 10.1371/journal.pgen.0020130
20
GorissenS. H.HristovaM.HabibovicA.SipseyL. M.SpiessP. C.Janssen-HeiningerY. M.et al (2013). Dual oxidase-1 is required for airway epithelial cell migration and bronchiolar reepithelialization after injury.Am. J. Respir. Cell. Mol. Biol.48337–345. 10.1165/rcmb.2012-0393OC
21
GrasbergerH.RefetoffS. (2006). Identification of the maturation factor for dual oxidase. Evolution of an eukaryotic operon equivalent.J. Biol. Chem.28118269–18272. 10.1074/jbc.C600095200
22
HallidayP. R.BlakelyC. M.BivonaT. G. (2019). Emerging targeted therapies for the treatment of non-small cell lung cancer.Curr. Oncol. Rep.21:21. 10.1007/s11912-019-0770-x
23
HatoT.TabataM.OikeY. (2008). The role of angiopoietin-like proteins in angiogenesis and metabolism.Trends Cardiovasc. Med.186–14. 10.1016/j.tcm.2007.10.003
24
HerbstR. S.MorgenszternD.BoshoffC. (2018). The biology and management of non-small cell lung cancer.Nature553:446. 10.1038/nature25183
25
HochbergY.BenjaminiY. (1990). More powerful procedures for multiple significance testing.Stat. Med.9811–818. 10.1002/sim.4780090710
26
HristovaM.HabibovicA.VeithC.Janssen-HeiningerY. M.DixonA. E.GeisztM.et al (2016). Airway epithelial dual oxidase 1 mediates allergen-induced IL-33 secretion and activation of type 2 immune responses.J. Allergy Clin. Immunol.1371545–1556. 10.1016/j.jaci.2015.10.003
27
HsuM. K.WuI. C.ChengC. C.SuJ. L.HsiehC. H.LinY. S.et al (2015). Triple-layer dissection of the lung adenocarcinoma transcriptome: regulation at the gene, transcript, and exon levels.Oncotarget628755–28773. 10.18632/oncotarget.4810
28
IdekerT.KroganN. J. (2012). Differential network biology.Mol. Syst. Biol.8:565. 10.1038/msb.2011.99
29
JiangZ.ZhouX.LiR.MichalJ. J.ZhangS.DodsonM. V.et al (2015). Whole transcriptome analysis with sequencing: methods, challenges and potential solutions.Cell. Mol. Life Sci.723425–3439. 10.1007/s00018-015-1934-y
30
JiangZ.GuiS.ZhangY. (2012). Analysis of differential gene expression in plurihormonal pituitary adenomas using bead-based fiber-optic arrays.J. Neurooncol.108341–348. 10.1007/s11060-011-0792-1
31
KayanoM.TakigawaI.ShigaM.TsudaK.MamitsukaH. (2011). ROS-DET: robust detector of switching mechanisms in gene expression.Nucleic Acids Res.39:e74. 10.1093/nar/gkr130
32
KimD.PerteaG.TrapnellC.PimentelH.KelleyR.SalzbergS. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.Genome Biol.14:R36. 10.1186/gb-2013-14-4-r36
33
KimY. J.BaekD. S.LeeS.ParkD.KangH. N.ChoB. C.et al (2019). Dual-targeting of EGFR and Neuropilin-1 attenuates resistance to EGFR-targeted antibody therapy in KRAS-mutant non-small cell lung cancer.Cancer Lett.46623–34. 10.1016/j.canlet.2019.09.00
34
KooS. H.FlechnerL.QiL.ZhangX.ScreatonR. A.JeffriesS.et al (2005). The CREB coactivator TORC2 is a key regulator of fasting glucose metabolism.Nature4371109–1111. 10.1038/nature03967
35
KwonO. S.KimW.ChaH. J.LeeH. (2019). In silico drug repositioning: from large-scale transcriptome data to therapeutics.Arch. Pharm. Res.42879–889. 10.1007/s12272-019-01176-3
36
LangfelderP.HorvathS. (2008). WGCNA: an R package for weighted correlation network analysis.BMC Bioinformatics9:559. 10.1186/1471-2105-9-559
37
LiY.JinS.LeiL.PanZ.ZouX. (2015). Deciphering deterioration mechanisms of complex diseases based on the construction of dynamic networks and systems analysis.Sci. Rep.5:9283. 10.1038/srep09283
38
LimS. Y.Gordon-WeeksA.AllenD.KersemansV.BeechJ.SmartS.et al (2015). Cd11b(+) myeloid cells support hepatic metastasis through down-regulation of angiopoietin-like 7 in cancer cells.Hepatology62521–533. 10.1002/hep.27838
39
LingQ.ShiW.HuangC.ZhengJ.ChengQ.YuK.et al (2014). Epigenetic silencing of dual oxidase 1 by promoter hypermethylation in human hepatocellular carcinoma.Am. J. Cancer Res.4508–517.
40
LittleA. C.ShamD.HristovaM.DanyalK.HeppnerD. E.BauerR. A.et al (2016). DUOX1 silencing in lung cancer promotes EMT, cancer stem cell characteristics and invasive properties.Oncogenesis5:e261. 10.1038/oncsis.2016.61
41
LiuW.OuyangS.ZhouZ.WangM.WangT.QiY.et al (2019). Identification of genes associated with cancer progression and prognosis in lung adenocarcinoma: Analyses based on microarray from Oncomine and The Cancer Genome Atlas databases.Mol. Genet. Genomic Med.7:e00528. 10.1002/mgg3.528
42
LuxenS.BelinskyS. A.KnausU. G. (2008). Silencing of DUOX NADPH oxidases by promoter hypermethylation in lung cancer.Cancer Res.681037–1045. 10.1158/0008-5472.can-07-5782
43
ManzoniC.KiaD. A.VandrovcovaJ.HardyJ.WoodN. W.LewisP. A.et al (2018). Genome, transcriptome and proteome: the rise of omics data and their integration in biomedical sciences.Brief. Bioinform.19286–302. 10.1093/bib/bbw114
44
MarianiL.BeaudryC.McDonoughW. S.HoelzingerD. B.KaczmarekE.PonceF.et al (2001). Death-associated protein 3 (Dap-3) is overexpressed in invasive glioblastoma cells in vivo and in glioma cell lines with induced motility phenotype in vitro.Clin. Cancer Res.72480–2489.
45
MazieresJ.PetersS.LepageB.CortotA. B.BarlesiF.Beau-FallerM.et al (2013). Lung cancer that harbors an HER2 mutation: epidemiologic characteristics and therapeutic perspectives.J. Clin. Oncol.311997–2003. 10.1200/JCO.2012.45.6095
46
OberndorferF.MullauerL. (2018). Molecular pathology of lung cancer: current status and perspectives.Curr. Opin. Oncol.3069–76. 10.1097/CCO.0000000000000429
47
RodónL.SvenssonR. U.WiaterE.ChunM. G. H.TsaiW. W.EichnerL. J.et al (2019). The CREB coactivator CRTC2 promotes oncogenesis in LKB1-mutant non-small cell lung cancer.Sci. Adv.5:eaaw6455. 10.1126/sciadv.aaw6455
48
ShamD.WesleyU. V.HristovaM.van der VlietA. (2013). ATP-mediated transactivation of the epidermal growth factor receptor in airway epithelial cells involves DUOX1-dependent oxidation of Src and ADAM17.PLoS One8:e54391. 10.1371/journal.pone.0054391
49
ShawA. T.SolomonB. J.BesseB.BauerT. M.LinC. C.SooR. A.et al (2019). ALK resistance mutations and efficacy of lorlatinib in advanced anaplastic lymphoma kinase-positive non-small-cell lung cancer.J. Clin. Oncol.371370–1379. 10.1200/JCO.18.02236
50
ShiS.LuoW.ZhangR.WangC.ZhengY.SongY.et al (2018). CRTC2 promotes non-small cell lung cancer A549 migration and invasion in vitro.Thorac. Cancer9136–141. 10.1111/1759-7714.12550
51
SiegelR. L.MillerK. D.JemalA. (2020). Cancer statistics, 2020.CA A Cancer J. Clin.707–30. 10.3322/caac.21590
52
SimabucoF. M.MorelloL. G.AragãoA. Z.Paes LemeA. F.ZanchinN. I. (2012). Proteomic characterization of the human FTSJ3 preribosomal complexes.J. Proteome Res.113112–3126. 10.1021/pr201106n
53
SzklarczykD.GableA. L.LyonD.JungeA.WyderS.Huerta-CepasJ.et al (2019). STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets.Nucleic Acids Res.47D607–D613. 10.1093/nar/gky1131
54
ThomasR.SrivastavaS.KatreddyR. R.SobieskiJ.WeihuaZ. (2019). Kinase-Inactivated EGFR is required for the survival of wild-type EGFR-Expressing cancer cells treated with tyrosine kinase inhibitors.Int. J. Mol. Sci.20:2515. 10.3390/ijms20102515
55
WangQ.ZhuW.XiaoG.DingM.ChangJ.LiaoH. (2020). Effect of AGER on the biological behavior of non-small cell lung cancer H1299 cells.Mol. Med. Rep.22810–818. 10.3892/mmr.2020.11176
56
WazirU.OrakzaiM. M.KhanzadaZ. S.JiangW. G.SharmaA. K.KasemA.et al (2015). The role of death-associated protein 3 in apoptosis, anoikis and human cancer.Cancer Cell Int.15:39. 10.1186/s12935-015-0187-z
57
WazirU.JiangW. G.SharmaA. K.MokbelK. (2012). The mRNA expression of DAP3 in human breast cancer: correlation with clinicopathological parameters.Anticancer Res.32671–674.
58
XiaoH.WangB.XiongH. X.GuanJ. F.WangJ.TanT.et al (2021). A novel prognostic index of hepatocellular carcinoma based on immunogenomic landscape analysis.J Cell Physiol2362572–2591. 10.1002/jcp.30015
59
XuT.RuanH.GaoS.LiuJ.LiuY.SongZ.et al (2020). ISG20 serves as a potential biomarker and drives tumor progression in clear cell renal cell carcinoma.Aging (Albany NY)121808–1827. 10.18632/aging.102714
60
YangL.LiY.WeiZ.ChangX. (2018). Coexpression network analysis identifies transcriptional modules associated with genomic alterations in neuroblastoma.Biochim. Biophys. Acta Mol. Basis Dis.18642341–2348. 10.1016/j.bbadis.2017.12.020
61
YinN. C.LangX. P.WangX. D.LiuW. (2015). AGER genetic polymorphisms increase risks of breast and lung cancers.Genet. Mol. Res.1417776–17787. 10.4238/2015.December.22.2
62
YooS. K.FreisingerC. M.LeBertD. C.HuttenlocherA. (2012). Early redox, Src family kinase, and calcium signaling integrate wound responses and tissue regeneration in zebrafish.J. Cell. Biol.199225–234. 10.1083/jcb.201203154
63
YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS16284–287. 10.1089/omi.2011.0118
64
YuX.ZengT.WangX.LiG.ChenL. (2015). Unravelling personalized dysfunctional gene network of complex diseases based on differential network model.J. Transl. Med.13:189. 10.1186/s12967-015-0546-5
65
ZhangB.HorvathS. (2005). A general framework for weighted gene co-expression network analysis.Stat. Appl. Genet. Mol. Biol.4:Article17. 10.2202/1544-6115.1128
66
ZhangC. C.KabaM.GeG.XieK.TongW.HugC.et al (2006). Angiopoietin-like proteins stimulate ex vivo expansion of hematopoietic stem cells.Nat. Med.12240–245. 10.1038/nm1342
67
ZhangH.WangS.HuangT. (2020). Identification of chronic hypersensitivity pneumonitis biomarkers with machine learning and differential co-expression analysis.Curr. Gene Ther.10.2174/1566523220666201208093325[Epub ahead of print].
68
ZhangW.FanJ.ChenQ.LeiC.QiaoB.LiuQ. (2018). SPP1 and AGER as potential prognostic biomarkers for lung adenocarcinoma.Oncol. Lett.157028–7036. 10.3892/ol.2018.8235
69
ZhuY.ZhouC.HeQ. (2019). High miR-139-3p expression predicts a better prognosis for hepatocellular carcinoma: a pooled analysis.J. Int. Med. Res.47383–390. 10.1177/0300060518802727
Summary
Keywords
WGCNA, differential correlation, switching mechanism, gene regulation, lung adenocarcinoma
Citation
Zhou Y, Xu B, Zhou Y, Liu J, Zheng X, Liu Y, Deng H, Liu M, Ren X, Xia J, Kong X, Huang T and Jiang J (2021) Identification of Key Genes With Differential Correlations in Lung Adenocarcinoma. Front. Cell Dev. Biol. 9:675438. doi: 10.3389/fcell.2021.675438
Received
03 March 2021
Accepted
24 March 2021
Published
05 May 2021
Volume
9 - 2021
Edited by
Liang Cheng, Harbin Medical University, China
Reviewed by
Xiao Chang, Children’s Hospital of Philadelphia, United States; Yi Zhang, Fudan University, China
Updates
Copyright
© 2021 Zhou, Xu, Zhou, Liu, Zheng, Liu, Deng, Liu, Ren, Xia, Kong, Huang and Jiang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jingting Jiang, jiangjingting@suda.edu.cnTao Huang, tohuangtao@126.com
This article was submitted to Molecular Medicine, a section of the journal Frontiers in Cell and Developmental Biology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.