For The Cancer Genome Atlas (TCGA) database, we downloaded the original transcriptome data from TCGA–GDC (TCGA–LUAD and –LUSC projects; STAR-Counts form). Before conducting the analysis, we merged and organized the downloaded raw transcriptional data into an expression matrix. A human genome reference document is used for ENSG number annotation. Meanwhile, we performed mean taking on duplicate genes, and genes with an average expression of less than 0.05 were deleted. The clinical formation was also obtained from the TCGA–GDC (bcr-xml form). For the Gene Expression Omnibus database, the data from GSE21656 and GSE108214 were selected, which provided the next-sequence data from cisplatin-resistant and wild-type lung cancer cells (Sun et al., 2012; Sarin et al., 2018). The probe annotation of GSE21656 was conducted using GPL6244, and GSE108214, using GPL17077. The baseline information on HNSCC patients from TCGA database is shown in Table 1.

The limma package was applied for differentially expressed gene (DEG) analysis with specific thresholds (Ritchie et al., 2015). By integrating patient expression profiles and prognostic data, univariate Cox regression analysis was utilized to identify the genes remarkably correlated with patient survival with a p < 0.05. Pathway enrichment was explored using the gene set enrichment analysis (GSEA) algorithm (Subramanian et al., 2005). The Gene Ontology (GO), Hallmark, and Kyoto Encyclopedia of Genes and Genomes (KEGG) gene sets were selected as the reference set. The immune microenvironment was quantified using the CIBERSORT algorithm based on the input expression matrix (Chen et al., 2018). Quantification of immune function was completed using the single-sample GSEA (ssGSEA) algorithm (Hänzelmann et al., 2013). Quantification of the stromal score, immune score, and estimate score was conducted using the ESTIMATE package (Yu et al., 2021a). The response of patients to immunotherapy was quantified using the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (Fu et al., 2020). The TIDE algorithm gives each patient a TIDE score. Lung cancer patients with TIDE scores more than zero were defined as responders to immunotherapy, and those with scores less than zero were the opposite. Quantification of patients on target drugs was realized using the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). The machine learning algorithm LASSO regression was utilized to reduce the data dimension (McEligot et al., 2020). The prognosis prediction model was identified using the multivariate Cox regression analysis. For a better clinical application, a nomogram that merges the risk score and clinical features was constructed. The prediction value of one certain continuous variable to survival was performed using the receiver operating characteristic (ROC) curve.

Single-cell analysis is directly completed through an online interactive website—the TISCH project (Sun et al., 2021). GSE148071 was conducted by Wu et al. (2021). In this project, they collected 42 NSCLC samples and characterized the entire tumor ecosystem through single-cell RNA sequencing.

The IHC image of PCDH7 in lung cancer and para-carcinoma tissue was directly downloaded from The Human Protein Atlas (HPA) project (Colwill and Gräslund, 2011). Subcellular localization of PCDH7 in the HPA project was obtained in the HeLa cells.

The used BEAS-2B, A549, H1299, H522, H460 cells, and A549-Res cell lines (A549 cell line that is resistant to cisplatin) were routinely stored in laboratories. The A549-Res (resistant to cisplatin) cell line was purchased from Shanghai MEIXUAN Biological Science and Technology Co, Ltd. All cells are cultured under normal conditions. The cell culture medium used is 1640-RPMI. According to cell growth, subculture was conducted every 3–4 days.

First, total RNA was extracted using the TRIzol reagent and then transcribed into cDNA for further analysis (Pan et al., 2020). The PCR system is a 20-μL system, and PCR detection is performed using SYBR Green. The primer sequence is as follows: PCDH7: F, 5′-GAC​TCT​GGG​CGT​CTC​TGA​AG-3’; R, 5′-CTC​AAC​TCC​GAC​TCT​GCT​CA; GAPDH: F, 5′-CTG​GGC​TAC​ACT​GAG​CAC​C-3’; R, 5′-AAG​TGG​TCG​TTG​AGG​GCA​ATG-3’.

The transfection of control and PCDH7 knockdown plasmids was conducted with Lipofectamine 2000 based on standard procedures (Pan et al., 2020). The target sequence is as follows: sh#1: 5′-GGA​GGC​TTC​TAA​GCC​AAA​T-3′, sh#2: 5′-GGA​CCA​TTT​ACT​CCA​CAA​T-3′, sh#3: 5′-CCA​CCG​TGG​TCC​TTA​ACA​T-3’.

The proliferation ability of A549 and H1299 cells was evaluated using the CCK8, colony formation, and EdU assays, according to standard procedures (Pan et al., 2020).

The invasion and migration ability of A549 and H1299 cells were evaluated using the transwell and wound-healing assays, according to standard procedures (Pan et al., 2020).

The detection of cisplatin was performed according to the procedures of Heinze et al. (2021).

For tumor inoculation, each nude mouse was subcutaneously injected with 10 × 10⁵ cells. The mice were then monitored for 4 weeks for tumor development. At the end of the experiment, the mice were euthanized, and tumors were excised, weighed, and photographed.

For all the analyses based on public data, R software version 4.0.4 was used. Moreover, SPSS and GraphPad Prism 8 software were also used for data statistics of the experimental data. Generally, the comparison with a p-value <0.05 was regarded as significant. Different testing methods were chosen based on statistical requirements for variables that meet different data distributions.

Through careful search, we found two datasets from the GEO database. The GSE21656 provided the transcriptional profile derived from cisplatin-resistant and wild-type H460 lung cancer cells. GSE108214 provided the transcriptional profile derived from cisplatin-resistant and wild-type A549 lung cancer cells. The data preprocessing process is shown in Supplementary Figures S1A, B. The DEG analysis identified 70 upregulated and 106 downregulated genes in cisplatin-resistant cells of the GSE21656 cohort (Figure 2A, H460). For the GSE108214 cohort, 1,430 upregulated and 1,280 downregulated genes were identified in cisplatin-resistant A549 cells (Figure 2B). Through intersection processing, we found that 16 genes showed consistent downregulation in drug-resistant cells in the GSE21656 and GSE108214 cohorts: NTS, TMPRSS15, TMEM27, MCAM, IGFBP3, S100A16, GALC, CDH11, DCLK1, MYO5C, CPVL, SEMA5A, ANO3, AQP3, IFITM2, and PCDH7; in total, 11 genes showed consistent upregulation in drug-resistant cells in the GSE21656 and GSE108214 cohorts: PKIA, CDH2, CALB2, CDK14, VAV3, KCNK1, COL12A1, ANO5, SNAP25, CP, and TPM2 (Figure 2C). Then, we compared the expression level of these common genes in NSCLC and para-carcinoma tissue. We noticed most of these genes had a significant difference in the expression level between tumor and normal tissues, revealing their underlying role in cancer development (Figure 2D). The results of univariate Cox regression showed that the genes PCDH7, TPM2, S100A16, CDH2, ANO3, CALB2, COL12A1, PKIA, and MCAM were remarkably correlated with NSCLC patient survival (Figure 2E, PCDH7: HR = 1.138; TPM2: HR = 1.148; S100A16: HR = 1.133; CDH2: HR = 1.109; ANO3: HR = 1.197; CALB2: HR = 1.063; COL12A1: HR = 1.065; PKIA, HR = 1.085; MCAM, HR = 1.115). Among these genes, PCDH7 has the most significant p-value and was consequently selected for the following analysis. In the TCGA–LUAD cohort (lung adenocarcinoma), patients with high expression of PCDH7 appear to have a poor prognosis (Figures 2F–H, overall survival, p = 0.003; disease-free survival, p = 0.002; progression-free survival, p = 0.003). However, this effect does not seem significant in the TCGA–LUSC cohort (lung squamous cell carcinoma) (Figures 2I, J). Furthermore, an association between PCDH7 and worse clinical features was found in a clinical correlation analysis (Figures 3A–D).

Subsequently, we tried to identify the biological role of PCDH7 in NSCLC. The results of qRT-PCR indicated that PCDH7 is overexpressed in lung cancer cells compared to normal lung cells (Supplementary Figure S2A). Moreover, the IHC image from the HPA database indicated a higher PCDH7 protein level in NSCLC than the control tissue (Supplementary Figures S3A, B). The inhibition efficiency of three sh-PCDH7 cells was quantified using qRT-PCR (Supplementary Figures S2B, C). In both A549 and H1299 cells, sh-PCDH7#2 showed the best performance and was consequently selected for further experiments. The transwell assay indicated that knockdown of PCDH7 could hamper the invasion and migration ability of NSCLC cells (Figure 3E). The result of the wound-healing assay obtained the same conclusion (Figure 3F). For cell proliferation, the results of the CCK8 assay indicated that inhibition of PCDH7 could reduce the cell proliferation ability of NSCLC cells (Figures 4A, B). In addition, the number and size of cell colonies in cells with PCDH7 knockdown were smaller than those in control cells (Figure 4C). The EdU assay indicated that inhibition of PCDH7 could remarkably inhibit the DNA replication capability of NSCLC cells (Figures 4D, E). In vivo experiments showed that PCDH7 knockdown cells formed a lighter tumor than control cells (Figures 4F, G). Then, we detected the IC₅₀ concentration of cisplatin in sh-PCDH7 and control cells. The result indicated that cells with sh-PCDH7 had a lower IC₅₀ than control cells, indicating that PCDH7 is associated with cisplatin resistance (Figures 5A, B). Moreover, A549 cells with cisplatin resistance had a higher PCDH7 level (Figure 5C).

Based on the public single-cell data from TISCH projects (Wu et al., 2021), we explored the single-cell expression pattern of PCDH7 in NSCLC. The result showed that PCDH7 was mainly expressed in malignant cells, fibroblasts, and CD8⁺ T cells (Figures 5D, E). The GSEA analysis based on the Hallmark set indicated that in the patients with high PCDH7 expression, the top five upregulated terms were epithelial–mesenchymal transition (EMT), apical junction, UV-response, angiogenesis, and TNF-α signaling (Figure 6A). For the GSEA analysis based on the GO set, the top three upregulated terms were all related to the spliceosome-related complex (Figure 6B) and the top three downregulated terms were related to the immunoglobulin complex (Figure 6C). For the GSEA analysis based on the KEGG set, the top three upregulated terms were small-cell lung cancer, focal adhesion, and ECM–receptor interaction, while the top three downregulated terms were ribosome, maturity-onset diabetes of the young, and linolenic acid metabolism (Figures 6D, E). The subcellular localization of PCDH7 in HeLa cells showed that it mainly localized in the plasma membrane (Supplementary Figures S3C, D).

The heatmap of the level of immune cells quantified by the CIBERSORT algorithm is shown in Figure 7A. Correlation analysis indicated that PCDH7 was positively correlated with resting NK cells, M0 macrophages, and neutrophils yet negatively correlated with CD8⁺ T cells and activated NK cells (Figure 7B). We also noticed a significant difference in several immune-related genes in patients with high and low PCDH7 expression levels (Figure 7C). For the immune function terms quantified by the ssGSEA algorithm, PCDH7 was positively correlated with para-inflammation, CCR, MHC_class_I, and APC_co_inhibition (Figure 7D). Furthermore, we found that PCDH7 was significantly correlated with key immune checkpoints PDCD1LG2 and CD274 (Figures 7E–H). A positive correlation was found between PCDH7 and stromal score, as well as ESTIMATE scores (Figures 8A–C). Moreover, we found that PCDH7 was positively correlated with the TIDE score (Figure 8D, R = 0.330, p < 0.001). Correspondingly, 44.6% of patients with low PCDH7 levels tended to respond to immunotherapy, but this percentage reduced to 26.9% in patients with high PCDH7 expression (Figure 8E). Meanwhile, a higher immune exclusion was found in patients with high PCDH7 expression (Figure 8F). The results of drug sensitivity indicated that patients with high PCDH7 expression might be more sensitive to bortezomib, docetaxel, and gemcitabine (Figure 8G).

Then, we tried to construct a prognosis model derived from PCDH7-related genes. The top 100 genes positively and negatively correlated with PCDH7 are shown in Figures 9A, B. Then, univariate Cox regression analysis was performed to identify the prognosis-related genes (Figure 9C; Supplementary File S1). LASSO regression analysis was conducted to reduce data dimensions (Figures 9D, E). Then, three PCDH7-related molecules were identified for a prognosis model: Risk score = GPX8 * 0.107 + BCAR3 * 0.184 + TNS4 * 0.05 (Figure 9F). The KM curve in the training cohort demonstrated a shorter survival rate of patients with a high risk score than those with a low risk score (Figure 10A, HR = 3.96, p < 0.001). The AUC values of 1-, 3-, and 5-year ROC curves were 0.750, 0.745, and 0.688, respectively, indicating a good prediction ability of our prognosis model (Figures 10B–D). This effect was also found in the validation cohort (Figures 10E–H, HR = 3.36, p < 0.001; AUC values of 1-, 3-, and 5-year ROC curves were 0.724, 0.733, and 0.673, respectively). A nomogram plot was established to get a better clinical application ability by integrating the risk score and clinical features (Figure 10I). For the 1-, 3-, and 5-year survival, a satisfactory fit was observed between the survival predicted by the nomogram and the actual survival (Figure 10J). Moreover, we noticed that the risk score is an independent marker for patient prognosis, which increases its potential for clinical applications (Figures 10K, L).