Antibodies used in this study included Vimentin (#5741), SLUG (#9585), β-catenin (#8480), and Anti-rabbit Alexa Fluor^® 594 Conjugate (#8889) from Cell Signaling Technology and MMP9 (10375-2-AP), SNAI1 (13099-1-AP), SNAI3 (21350-1-AP), and SETDB2 (14428-1-AP) from Proteintech. The siRNAs were obtained from Shanghai GenePharma.

The cell lines were purchased from American Type Culture Collection (ATCC). The cell lines BEAS-2B, A549, HCC827, H1299, H1975, EPLC-272H, H226, H157, and H2170 were cultured in a humidified incubator at 5% CO₂ and 37°C using RPMI1640 medium with 10% fetal bovine serum (FBS) and 1% Pen-Strep penicillin-streptomycin (Gibco). Transfections were performed with RNAiMAX Reagent (Invitrogen) according to the manufacturer’s instructions. Each experiment was performed at least three times. siRNA against SETDB2 (5’-GUGUACGCUGUCUAGAUGATT-3’), siRNA against SNAI3 (5’-GACGCAGAGAGAAAUCAAUTT-3’) and siRNA negative control (5’-UUCUCCGAACGUGUCACGUTT-3’) were obtained from Shanghai GenePharma.

The research route is shown in Figure 1 . The TCGA dataset referenced in the study is available in a public repository from UCSC Xena (https://xena.ucsc.edu/), including transcriptome expression RNA-seq and clinical data of 1076 NSCLCs. The expression profile and survival information of patients with NSCLC, including 130 samples, were downloaded from the UCSC database (https://ucscpublic.xenahubs.net). Transcription factor sets were defined using Trrust (https://www.grnpedia.org/trrust/), CISBP (http://cisbp.ccbr.utoronto.ca/index.php), and JASPAR (http://jaspar.genereg.net/). In total, 1942 (after de-redundancy) human TFs were downloaded from the database for subsequent analysis.

Data from 130 samples from UCSC were used to verify the efficacy of the TFs prognostic scores. After standardizing the maximum and minimum data values, a log-rank test was performed on the overall survival (OS) data, and a Kaplan–Meier curve was drawn. p < 0.05 was considered a statistically significant correlation. The survival receiver operating characteristic (ROC) software package was used to plot the ROC curve, observe the prognostic efficacy with 5-year survival as the threshold, and calculate the area under the curve (AUC) score. The ROC curve was used to describe the sensitivity and specificity of the survival prediction based on the TFs score.

Differentially expressed genes (DEGs) between normal and tumor samples were estimated using the limma R package. Using a Bayesian approach to computation, candidate DEGs were identified by p < 0.05 and |fold-change| > 1.5. After determining the intersection with the TF set, the “pheatmap” R package was used to draw a heatmap of partial differential TF gene expression in tumor and normal tissues. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID: https://david.ncifcrf.gov/).

The “WGCNA” R package (25) was used to perform co-expression network analysis on TCGA containing differentially expressed TFs. First, no abnormal samples were found during sample clustering. A picksoft threshold was used to select the optimal soft threshold β for subsequent network construction; the expression profile was analyzed using scale-free clustering and a dynamic shearing tree to identify co-expressed gene modules. The co-expression network matrix calculation was used on the Pearson correlation. In addition, the relationship between these gene modules and clinical phenotypes was positioned the module as a survival module; the genes of these modules were used as candidate prognostic-related TFs.

For the cancer samples of TCGA-NSCLC (997 samples after removing the no-OS data), the TF genes were analyzed using Cox single-factor survival analysis based on the survival package; the prognostic genes with single-factor significance (p < 0.05) were selected as candidate genes. The Kaplan–Meier method was used to establish survival curves, and the log-rank test was used to calculate the significance of the differences. Forest plot software was to plot a prognostic forest plot of the Cox single-factor results of related genes. Subsequently, prognosis-related TFs were included in the LASSO regression analysis using the “glmnet” R package to obtain the prognostic-related TFs, obtain the corresponding genes’ regression factor parameters, and build a TFs prognosis model based on this. Finally, Cox proportional hazards regression was used to analyze the independent predictive power of the TFs score.

Total cellular RNA was extracted from cells using an RNA-Quick Purification Kit (Yishan Biotechnology). cDNA was prepared using a Transcriptor First Strand cDNA Synthesis Kit (Roche, Basel, Switzerland). Reverse transcriptase polymerase chain reaction of the selected genes was performed by reverse transcriptase polymerase chain reaction using ABI QuantStudio 5 (Applied Biosystems). In addition, SYBR Green fluorescence was measured and quantified using the comparative Ct method (2^-ΔΔCt), with GAPDH expression as an internal control. The assay was performed in triplicate, and the primers used are listed in Table 1 .

Cells were inoculated into 12-well plates. Next, the cells were processed using the EdU cell proliferation kit (C10310 RiboBio). Finally, the assay was performed with a fluorescent microscope. The whole process was carried out according to the manufacturer’s instructions.

Cells were taken at logarithmic growth stage, spread in a 6-well plate, and waited for the cells to grow above 95% confluence for scratching. Cells were washed 3 times with PBS and photographed; at this point the photo was taken for 0 h. The cells were added to the medium and incubated, and the healing of the scratches was observed and photographed. The images were used to analyze the percentage of healing of the scratches. Images at time zero (0 h) and at 24 or 36 h (Δ h) were captured. The area of wound was quantified by Image J software. The percentage of wound closure: percentage of wound closure = [(A (0 h) – A (Δ h)/A (0 h)] × 100%, A (0 h) is the area of wound measured immediately after scratching, and A (Δ h) is the area of wound measured 24 or 36 h after scratching.

Transwell chamber filters were wrapped with Matrigel (BD Biosciences, Franklin Lakes, NJ). Lung cancer cells transfected with specific siRNAs were suspended in serum-free medium, and the cells were seeded into the upper chamber of the Transwell. The lower chamber was filled with medium containing 10% FBS. After 24 h of incubation, the cells were fixed with methanol and stained with crystal violet solution, after which the cells in the upper chamber were wiped off with a cotton swab and the remaining cells were photographed. Five high-power fields of view were taken for each small chamber selection.

Cells were cultured in confocal-specific dishes until 85% confluence was reached, when they were fixed in 2% paraformaldehyde fixative, after which the cells were washed 3 times with PBS. Cells were treated with 0.2% triton x-100 for 5 min at 25°C, washed 3 times with PBS, and closed with 0.8% bovine serum albumin for 1 h at 25°C. Cells were incubated overnight or for 1 h at 25°C by primary antibody and then they were washed 2 times (10 min each) at 25°C with 0.1% triton x-100. Cells were incubated for 1 h with the addition of secondary antibody, and 0.1% triton x-100 was added to wash twice for 10 min each time at 25°C. Finally, fluoroshield with 4′,6-diamidino-2-phenylindole (DAPI) was added to cover the bottom of the dish. The mean fluorescence intensity (MFI) of cells was calculated from Image J software.

Cell-culture dishes were placed on ice, and cells were washed with ice-cold PBS. Cells were added to ice-cold Radio Immunoprecipitation Assay (RIPA) lysate and lysed at 4°C for 10 min, vortexing every 5 min. The BCA was quantified, and protein lysate was added to the loading buffer for 10 minutes at 95°C. Electrophoresis was carried out at 80 V until bromophenol blue ran through the top layer of the gel, followed by a shift to 120 V until bromophenol blue ran out of the sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) gel. To transfer the proteins from the gel to the membrane, the polyvinylidene fluoride (PVDF) was activated with methanol for 1 min and rinsed with transfer buffer in a “sandwich” filter paper and sponge configuration, sponge–filter paper–gel-membrane filter paper–sponge. The PVDF membrane was placed in 5% skimmed milk for 1 h at 25°C. The PVDF membranes were washed with tris buffered saline with tween 20 (TBST), and then primary antibody (diluted in primary antibody diluent) was added and incubated overnight at 4°C. The membranes were washed 3 times in TBST for 5 minutes each time, secondary antibody prepared in 5% skimmed milk was added, and the membranes were incubated for 1 h at 25°C in a gentle mix. Membranes were washed 3 times in TBST for 5 minutes each time. Chemiluminescence images were obtained using a darkroom development technique. Western Blot signal was quantified using the Image J software.

GraphPad Prism (version 8.0) was used to analyze the data from our study. Results are expressed as the mean ± SD of at least three individuals. Student’s t-test was used to compare differences between the two groups. Survival curves were plotted using the Kaplan–Meier “survival” package in R (version 3.4.3). Log-rank test was used to assess statistical significance. Statistical significance was set at p < 0.05.

Data containing 1017 NSCLC and 59 normal samples from the TCGA dataset with clinical and RNA-seq datasets were used. Differential analysis was performed to find differentially expressed genes (DEGs) between normal and NSCLC samples, with the criteria of |fold change| >1.5 and p < 0.05; 8840 DEGs were identified. In total, 725 TFs ( Supplemental Table 1 ) were identified as NSCLC-related TFs ( Figure 2A ), of which the top 200 TFs were upregulated in NSCLC compared to normal tissues. The heatmap of the top 200 differentially expressed TFs is shown in Figure 2B . To further determine the function of these TFs, GO and KEGG pathway annotations were used to analyze the 725 differential TFs. These genes could be enriched in tumor malignancy-related GO biological processes, containing cell differentiation, proliferation, apoptosis, stemness, and angiogenesis ( Figure 2C ). Moreover, KEGG pathway enrichment analysis showed that 725 genes were involved in many pathways, including tumor transcriptional regulation interrelated to oncogenesis and development, the Hippo signaling pathway, TGF-β signaling pathway, cell cycle, immune-related Th17 cell differentiation, and Th1 and Th2 cell differentiation ( Figure 2D ).

We used WGCNA to build a co-expression network and modules of NSCLC-correlated differentially expressed TFs. A co-expression network was constructed with several modules. A scale-free network evaluation coefficient threshold > 0.9 was selected to make the co-expression network conform to the scale-free network standard. The 725 transcription factors were subjected to clustering analysis to classify them into six modules ( Figures 3A, B ) ( Table 2 ). The six modules were divided into two categories: green, brown, and yellow made up one category, while the other category included blue and turquoise ( Supplemental Table 2 ). Internal modules in the same category were positively correlated ( Figure 3C ). The correlation analysis of each module and clinical data was performed to find transcription factor modules that are highly correlated with clinical survival. There was a correlation between the gene significance of OS and module membership of the green module (cor = −0.071, p = 0.03). The OS time negatively correlated with “brown” (cor = −0.063, p = 0.05) and “yellow” (cor = −0.087, p = 0.006) ( Figure 3D ). The module membership of the green module was significantly related to the gene significance of the OS, and the module membership of the brown and yellow modules was significantly related to the gene significance of the OS time ( Figure 3E ). Therefore, 118 genes contained in these three modules (green, brown, and yellow) were selected as candidate survival-related TFs for further analysis. 118 genes were enriched in many immune-related pathways, including T-cell differentiation, lymphocyte differentiation, and T-cell cytokine production. Furthermore, the potential regulatory molecular functions of these TFs were found to involve their classical functions, including enhancer binding, histone deacetylation, nuclear receptor activity, and histone acetyltransferase ( Figures 3F, G ) ( Supplemental Table 3 ).

To search for TFs with key roles in the NSCLC prognosis, Cox single-factor regression analysis was performed on the above 118 key TFs; 9 TFs significantly (p < 0.05) related to survival were selected ( Figure 4A ). According to previous research, 21 genes highly related to NSCLC cell proliferation and OS of patients with NSCLC in the study were added, including 11 potential tumor suppressor genes (AFF3, AhR, AR, CBFA2T3, CHD4, KANK2, NR3C2, PTEN, PRDM16, RB1, and STK11) and 10 potential oncogenes (BARX1, DLX6, ELF3, EN1, ETV1, FOXBE1, IRX4, IRX5, and SALL1) (30). We performed LASSO regression analysis for these 30 TFs. Finally, five prognostic-related TFs (SETDB2, SNAI3, SCML4, ZNF540, and ETV1) were screened using LASSO regression ( Figure 4B ). The LASSO regression coefficients of each gene were shown in Figure 4C . Subsequently, clinical data were analyzed in relation to TFs scores, and the analysis revealed that the number of packs smoked per year (p = 0.003017), gender (p = 2.899365e-08), pathological stage TNM (T: p = 0.0004997501, N: p = 0. 006496752, M: p = 0.004497751), and targeted molecular therapy (p = 0.02204301) had significant differences in the distribution of clinical characteristics in the high and low TFs score groups. Conversely, there was no significant trend in other clinical characteristics ( Figure 4D ).

Furthermore, in the ROC analysis score for prognosis, the 3-year survival AUC value was 0.868, and that for the 5-year survival was 0.731 in the ROC curve ( Figure 4E ). To verify the prognostic model’s validity, the prognostic performance of the TFs score on another set of UCSC (https://ucscpublic.xenahubs.net) data containing 130 samples from NSCLC was analyzed. The TFs score was shown to have a predictive function for the survival of NSCLC patients (p < 0.037) ( Figure 4F ).

Finally, Cox multifactor regression was used to analyze the independent predictive performance of the model and explore the effect of other clinical factors on the prognosis in the TF score. We evaluated the prognostic value of the five TFs according to score, age, sex, number of packs smoked per year, N stage, radiotherapy, and targeted therapy. Univariate regression analysis showed that the TFs score was significantly associated with OS (HR = 1.2869, p = 0.0195); multivariate regression analysis showed that the TFs score was significantly correlated with OS (HR = 1.2681, p = 0.0348) ( Figures 5A, B ). The nomogram quantified the contribution of individual factors to clinical prognosis to verify the model’s validity, which had good predictive power ( Figure 5C ) ( Supplemental Table 4 ).

One important function of the immune system is to recognize and subsequently destroy tumors. In recent years, progress in tumor treatment has benefited from research involving immunotherapy (31–33). Therefore, TME plays a vital role in diagnosing and treating tumors. According to the TCGA-NSCLC dataset, we analyzed the differences in the immune microenvironment between the high- and low-score groups ( Figure 6A ). Using CIBERSORT calculations, significant differences were found in the infiltration levels of 22 immune cells in the high- and low-score groups. The immune cells with significant differences were memory B cells, naive CD4 T cells, resting memory CD4 T cells, regulatory T cells, resting natural killer (NK) cells, activated NK cells, monocytes, M0 macrophages, M1 macrophages, resting dendritic cells, activated dendritic cells, resting mast cells, activated mast cells, eosinophils, and neutrophils ( Figure 6B ) ( Supplemental Table 5 ).

To further explore the expression of the five TFs in clinical samples, LUAD and LUSC data from TCGA were used to compare their expression in normal and NSCLC samples. The results showed that the expression levels of SETDB2, SNAI3, SCML4, ZNF540, and ETV1 in LUAD ( Figure 7A ) and LUSC ( Figure 7B ) were significantly lower than those in normal tissues. Meanwhile, survival analysis of each TF showed that lower SETDB2, SNAI3, SCML4, and ZNF540 expression was associated with poor prognosis in lung cancer; higher expression of ETV1 was associated with poor prognosis in lung cancer via Affymetrix microarray data and via RNA-seq data ( Figures 7C, D ) (https://kmplot.com/analysis). It was worth noting that the expression level of ETV1 in normal tissues was higher than that in lung cancer tissues; however, patients with high expression had poor prognoses. This phenomenon suggested two roles of ETV1 in the occurrence and development of lung cancer. The TFs likely interacted with other proteins to form complexes that activated or inhibited the transcriptional regulation of genes and downstream genes to function (19). Therefore, the interaction protein network of the five TFs were analyzed using STRING ( Figure 7E ) ( Supplemental Table 6 ) (https://string-db.org). These data showed that SETDB2, SNAI3, SCML4, and ZNF540 have tumor suppressor functions in lung cancer.

To gain further support for the notion that SETDB2, SNAI3, SCML4, ZNF540, and ETV1 regulate the malignant phenotypes of lung cancer through the transcriptional repression of genes, experiments were performed. First, to further explore the function of SETDB2, SNAI3, SCML4, ZNF540, and ETV1 in lung cancer cell lines, the mRNA expression levels in the normal human bronchial epithelial cell line (BEAS-2B), lung adenocarcinoma cell lines (A549, HCC827, H1299, and H1975), and lung squamous cell lines (EPLC-272H, H226, H157, and H2170) were measured ( Figure 8A ). Consistent with the results of bioinformatics analysis, the expression levels of SETDB2 and SNAI3 were higher in BEAS-2B than in lung cancer cell lines ( Figure 8B ).

To determine how SETDB2 and SNAI3 regulate the cell proliferation of NSCLC, EdU assays were performed, which showed that knockdown of SETDB2 or SNAI3 in A549 cells had a strong promotion effect on the proliferation of lung cancer cells ( Figure 8C ). We then investigated the roles of SETDB2 and SNAI3 in cellular behavior of NSCLC using wound healing and transwell invasion assay in A549, HCC827, and H1299, which were transfected with control siRNA, SETDB2 siRNA, or SNAI3 siRNA. The wound healing assay showed that depletion of SETDB2 or SNAI3 promoted the migration of lung cancer cells ( Figure 8D ). Meanwhile, the increased invasion was induced by knocking-down SETDB2 or SNAI3 ( Figure 8E ).

We next investigated the possibility of SETDB2 and SNAI3 in NSCLC metastasis and stemness in vitro. The result of immunofluorescence assay indicated that metastasis-related markers (MMP9, Vimentin, β-catenin, SNAI1, and SLUG) were significantly upregulated with SETDB2 or SNAI3 knockdown in A549 cells ( Figure 8F ). The expression of metastasis-related markers was analyzed by RT-qPCR and Western blotting using A549 cells transfected with control siRNA, SETDB2 siRNA, or SNAI3 siRNA. SETDB2 or SNAI3 depletion led to an increase of metastasis-related markers ( Figures 8G, H ). In addition, the stem cell markers (SOX2, OCT4, and NANOG) were upregulated with the knockdown of SETDB2 or SNAI3 ( Figure 8I ). These results show that SETDB2 and SNAI3 inhibit proliferation, migration, invasive, metastasis, and cell stemness of lung cancer cells.