We downloaded 484 LUSC and 510 LUAD cohorts from the cBioPortal for Cancer Genomics (http://www.cbioportal.org/) (17). We then used R to standardize the dataset and determine the comparability of the data. By calculating the median absolute deviation, we selected the top 20% of genes shared by the two cancers. Then, according to the degree of sample aggregation estimated by the WGCNA package, 136 LUSCs and 184 LUADs were retained (Supplementary Table S1) (18). The two datasets were comparable when the data correlation between LUSCs and LUSCs was 0.46 (P < 0.01) (Supplementary Figure S1).

We set genes as nodes and relationships between genes as edges to build a co-expression network using the WGCNA package (18). The construction steps of the network mainly included correlation matrix calculation, soft threshold selection, adjacency matrix calculation, heterogeneous matrix calculation, dynamic branch cutting, and module merging (19). We then identified the network modules.

We identified the consensus module between LUSC and LUAD. Module overlap degree referred to the common gene number between modules (20). We then generated Z-scores (20) using the permutation test in the WGCNA package to assess the preservation of LUAD genes in the LUSC module. The ranges 5 < Z < 10 and Z > 10 were considered moderate and highly preserved, respectively.

Module eigengene (ME), the first principal component of modules, represents the feature expression mode of modules (21). Eigengene connectivity (KME) represents the Pearson correlation between genes (including genes not originally assigned to modules) and modules in the network. We calculated KME and P values for all genes in LUSC and LUAD, then ranked the KME values from largest to smallest and selected the top 100 genes for each module (P < 0.05). We performed a Gene Ontology (GO) enrichment analysis to select the top five enrichment terms for each module via the Clusterprofiler package (22).

Based on Z-scores, we selected candidate hub genes with KME-module correlations greater than 0.9 in high-conservation modules. We then performed Spearman's rank correlation analysis between the expression of candidate hub genes and the overall survival time of patients. Finally, the online software Kaplan–Meier Plotter (23) was used to analyze the effect of gene expression on patient survival. Cancer samples were divided into two groups based on high and low expression levels of genes. Genes with significantly different survival curves between the two groups were thought to be closely related to the survival of cancer patients (19). We used Cox univariate analysis in the SURVIVAL and SURVMISER software packages to verify whether hub gene expression is a major prognostic factor.

Based on the GEO database (https://www.ncbi.nlm.nih.gov/gds/), we obtained the LUSC and LUAD gene expression dataset (GSE134381), including 20 LUSC cancer samples and 20 normal samples, and 17 LUAD cancer samples and 17 normal samples. We evaluated the expression differences of hub genes in cancer and normal samples. We searched for chemical substances and human diseases that have regulatory relationships with hub genes in the Comparative Toxicogenomics Database (CTD) (24) and counted the number of corresponding files.

We identified 11 LUSC and 6 LUAD gene co-expression modules, respectively (Figure 1 and Supplementary Figures S1, S2). For visualization, we named modules with colors. Modules are clusters of densely interconnected genes that may be involved in a similar function. For example, the LUSC module “greenyellow” contains some genes, such as AKR1C1, AKR1C2, and AKR1C3, enriched in regulation of the extent of cell growth pathways. In addition, the LUAD module “yellow” consists of genes, such as CCL19, CCR7, CCL21, and LY9, that are enriched in T-lymphocyte-related pathways involved in T-cell activation, regulation of lymphocyte activation, and regulation of T-cell activation.

Modules highlighted similar expression patterns of genes in LUSC and LUAD, respectively. We then assessed how well the modules in one cancer subtype are preserved in another cancer subtype. As a qualitative evaluation, we imposed modules from LUSC onto the network for dataset LUAD and vice versa (Figures 1A,B). Some modules from one cancer subtype are still assembled in another cancer, well preserved (Figure 2; Supplementary Table S2). There is a high degree of gene overlap within the module between LUSC and LUAD. For example, the LUSC module “turquoise” is contained in the LUAD module “turquoise,” and the LUSC module “brown” is corresponding to the LUAD module “blue.” To quantify this result, we used the Z-score to summarize how well the modules hold in each other (Table 1). In general, the higher the value of the Z-score, the more preserved the module is between cancer subtypes. The Z-scores of the six modules are greater than 10, indicating that the LUSC network is significantly similar to LUAD (Table 1).

We explored the top five significant rich GO terms for each module. The functional overlap degree of modules between LUSC and LUAD corresponds to the degree of preservation of the module among them (Figure 3 and Supplementary Table S3). The higher the degree of preservation between the LUSC module and the LUAD module, the higher the functional overlap between them. For example, the LUSC module “turquoise” and the LUAD module “turquoise” have the highest Z-score, sharing 80% (4/5) of GO terms (Figure 3A). Furthermore, LUSC module “brown” and LUAD module “blue” are almost completely overlapping (Figure 2), which share common biological processes, such as extracellular matrix organization, extracellular structure organization, collagen fibril organization, ossification, and cartilage development (Figure 3B). In addition, low-preservation modules also have shared GO terms. For example, both the LUSC module “magenta” and the LUAD module “blue” are involved in extracellular structure in organization processes (Figure 3B).

The LUSC and LUAD networks share 46 hub genes with connectivity greater than 0.9. We then correlated these genes with patient survival time. The hub genes of COL5A2, TTLL3, SPEF1, TMEM190, CCDC65, CCDC33, and GLT8D2 have the highest Spearman rank correlations. Furthermore, COL5A2 is closely related to both LUSC and LUAD patients’ survival time (Figure 4), implying it may be a risk factor. The COL5A2 gene may lead to both LUSC and LUAD patients' poor prognosis [LUSC HR = 1.35 (1.02–1.78), logrank P = 0.035; LUAD HR = 1.5 (1.07–2.11), logrank P = 0.019]. The higher the COL5A2 expression, the worse the patient survival rate (Figure 4).

To determine whether the COL5A2 gene is a major prognostic factor, we correlated the clinical data of LUSC and LUAD with patient prognoses (Table 2). The survival analysis results showed that age and gender have no significant relationship with LUSC and LUAD patients' survival rates (Figure 5). Moreover, the COL5A2 gene expression in LUSC and LUAD samples is significantly higher than that in normal samples (P < 0.001, Figure 6). Therefore, COL5A2 may serve as a potential therapeutic target for NSCLC.

To investigate the role of COL5A2, we searched the CTD database. We found that COL5A2 may regulate chemicals and human diseases and it is related to tumor treatment drugs. In addition, many tumor-related substances or drugs can affect the expression of COL5A2 (Figure 7). For example, PIRINIXIC acid can increase COL5A2 expression by binding to PAPAYA protein and increasing its activity. In addition, tumor preventive drugs such as decitabine and fenretinide can downregulate COL5A2. Dexamethasone can reduce COL5A2 expression and suppresses the testosterone response.