We downloaded the miRNA and mRNA sequencing expression profiles and associated clinicopathological data of TCGA-LUSC from the GDC data portal at the National Cancer Institute (https://portal.gdc.cancer.gov/). We searched lung squamous cell carcinoma relevant gene microarray expression datasets from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) with the following keywords: “lung squamous cell carcinoma”. Filters were set to “series” and “Expression profiling by array”, “Non-coding RNA profiling by array” and “Homo sapiens”. The inclusive criteria of datasets is that the datasets have miRNA or mRNA expression values in lung squamous cell carcinoma tumors and normal tissues. Then we screened 5 eligible miRNA datasets and 10 eligible mRNA datasets. First, we screened differentially expressed miRNA (DE-miRNAs) or mRNA (DE-mRNAs) in each dataset. The DE-miRNAs and DE-mRNAs were obtained from microarray expression profiles using the web analysis tool GEO2R, which is used to compare groups of samples by the GEOquery and limma R packages from the Bioconductor project in the GEO database (http://www.ncbi.nlm.nih.gov/geo/geo2r/). The cut-off criteria were p-value < 0.05 and |log2 (fold change) | ≥ 1. Then, we summarized the differential miRNAs or mRNAs screened out from each dataset. We also collected differentially expressed miRNAs from miRCancer and Database of Differentially Expressed MiRNAs in human Cancers (dbDEMC) (Xie et al., 2013; Yang et al., 2017). An overview of the workflow steps is shown in Figure 1.

TarBase (Karagkouni et al., 2018) (http://www.microrna.gr/tarbase) is a reference database, specifically for index experiments support the miRNA targets. It integrates information on cell-type specific miRNA–gene regulation, while hundreds of thousands of miRNA-binding locations are reported. miRTarBase (Chou et al., 2018) is a comprehensive annotated and experimentally validated miRNA-target interaction database for miRNA-related researches. Tarbase and miRTarBase are used to construct the miRNA-mRNA regulatory networks. We analyzed the correlation between miRNA and mRNA in TCGA-LUSC, and screened the miRNA-mRNA regulatory networks according to the results. The online program DAVID (http://david.abcc.ncifcrf.gov/), a database for annotation, visualization and integrated discovery, is a comprehensive tool for researchers and scholars to understand the biological meaning behind multiple genes. We used DAVID, DIANA-MirPath (Vlachos et al., 2015), a miRNA pathway analysis web-server, and Hiplot (https://hiplot.com.cn), a comprehensive web platform for scientific data visualization, to perform Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes Genomes (KEGG) pathways analysis. Single sample gene set Enrichment analysis (ssGSEA) is an extension of the GSEA method that allows the definition of an enrichment score that represents the absolute degree of enrichment of the gene set in each sample within a given dataset (Hoadley et al., 2018; Xiao et al., 2020). The ssGSEA data of TCGA-LUSC were downloaded from UCSC Xena (https://xena.ucsc.edu/) to analyze the possible enrichment pathways of DE-miRNAs and DE-mRNAs.

FFPE specimens were obtained from LUSC patients undergoing radical surgery in the First Affiliated Hospital of Nanjing Medical University. All samples were collected with the written informed consent from the patient and the prior approval of the Institutional Review Boards of the First Affiliated Hospital of Nanjing Medical University (ID: 2016-SRFA-148) according to the Declaration of Helsinki. The clinical features of the 30 LUSC patients are shown in Table 1. Total RNA was extracted from FFPE specimens using the TIANGEN RNAprep Pure FFPE kit (Tiangen, Beijing, China). The acquired RNA from each sample was lysed in 40 µl RNase-free water and stored at −80°C until use. The concentration and purity of RNA were analyzed by the Nanodrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States ) (A260/A280 = 1.8–2, A260/A230 > 1.7).

According to the manufacturer’s instructions, the external validation was verified by qRT-PCR using PrimeScript RT reagent Kit (Takara) and SYBR Premix Ex Taq II (Takara) after adding a poly (A) tail to RNA by Poly (A) Polymerase Kit (Takara). The sequences of PCR primers are listed in Supplementary Table S1. This process was run on qTOWER³ 84 (Analytik Jena) at 95°C for 20 s, followed by 40 cycles of 10 s at 95°C, 20 s at 60°C. The expression levels of miRNAs and mRNAs in tissue were calculated using the 2^−ΔΔCt method (Livak and Schmittgen 2001). (U6 as endogenous reference miRNA for sample normalization; 18s rRNA as endogenous reference mRNA for sample normalization; ΔCt = Ct_miRNA/mRNA− Ct_normalizer; Ct: the threshold cycle).

We downloaded the infiltrating immune cell types data from the TCGA website and calculated using CIBERSORT (https://cibersort.stanford.edu/index.php/), a general calculation method used to quantify cell fractions from a large number of tissue gene expression profiles (GEP). Combining support vector regression with prior knowledge of the expression profile of purified white blood cell subsets, CIBERSORT was able to accurately estimate the immune component of tumor biopsy (Chen et al., 2018). The stromal and immune levels of TCGA-LUSC samples were evaluated using ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data) software. This method uses gene expression signatures to infer the fraction of stromal and immune cells in tumour samples (Yoshihara et al., 2013). We used the UCSC Xena platform (https://xena.ucsc.edu/) to obtain the methylation level of CpG sites in TCGA-LUSC samples, which is detected by the Illumina Infinium HumanMethylation450 BeadChips platform. The amount of somatic variants per megabase (MB) of the genome was measured by tumor mutation burden (TMB).

We used GraphPad Prism software, IBM SPSS Statistics v.26 software (IBM Corporation, Armonk, NY, United States ), and R language v3.6.3 (https://cran.r-project.org/) for data analysis. The t-test was performed to analyze DE-miRNAs and DE-mRNAs. The statistically significant criteria of DE-miRNAs and DE-mRNAs are |log2FC| >0.58 and p < 0.05. The area under the ROC curve (AUC) and decision curve analysis (DCA) based on logistic regression were used to evaluate the diagnostic efficacy of miRNA-mRNA networks. The Pearson correlation method was used to calculate the correlation between DE-mRNAs or DE-miRNAs and tumor-related phenotypes. Prognostic analysis was performed using the “survival” package.

We screened 5 miRNA and 10 mRNA expression datasets from GEO, and the information of 15 GEO datasets is shown in Table 2. The log2 fold change (LUSC vs. NC) was used to screen the DE-mRNAs and DE-miRNAs. As shown in Figure 2, the 7 DE-miRNAs and 270 DE-mRNAs were the intersections of the TCGA, GEO, dbDEMC, and miRcancer databases. The results of the preliminary screening are presented in Supplementary Table S2. We utilized DIANA-MirPath to predict the possible functions of the 7 DE-miRNAs, as well as the targets of DE-miRNAs enriched in Fatty acid biosynthesis, Cell cycle, Fatty acid metabolism etc,. (Supplementary Figure S1A). Then, we also analyzed the 270 DE-mRNAs using cluster profile on the Hiplot platform. KEGG pathway enrichment analysis revealed that the dysregulated DE-mRNAs were enriched in the cell cycle, DNA replication, p53 signaling pathway etc,. (Supplementary Figure S1C). The top 10 GO terms of the DE-miRNAs and DE-mRNAs were presented in Supplementary Figures S1B,D, including cellular component, extracellular exosome and extracellular space in the CC category, DNA metabolic process, DNA replication and cell division in the BP category, protein binding transcription factor activity, protein binding, and ATP binding in the MF category.

We further screened miRtarbase and Tarbase database to establish potential miRNA-mRNA regulatory networks, and screened out 5 miRNA-mRNA regulatory networks (Figure 3A). The complete prediction network of miR-205-5p is shown in Supplementary Table S3. Then, we filtered out 4 miRNA-mRNA networks with significant correlation (adjusted p-value < 0.05) in TCGA-LUSC listed in Figure 3B and full statistical results listed in Supplementary Table S4.

In order to prove the differential expression of the DE-miRNAs and DE-mRNAs, we further validated them in 30 pairs of matched tumors and adjacent normal tissues by qRT-PCR. As shown in Figure 4, the expression of miR-205-5p (p < 0.0001) and UBE2C (p = 0.0093) were upregulated in tumor tissues, while the expressions of miR-140-3p (p = 0.0293), PTPRM (p < 0.0001), GPD1L (p = 0.0002), and FOXF2 (p < 0.0001) were downregulated in tumor tissues. At the same time, miR-182-5p (p = 0.2286) and miR-210-3p (p = 0.0879) showed no significant difference in tumor tissues compared with normal tissues. Spearman correlation analysis of the interaction between DE-miRNAs and DE-mRNAs, showed that miR-205-5p was significantly correlated with PTPRM (p = 0.0186, r = −0.3031). IHC images from the HPA database showed that PTPRM expression in LUSC was lower than in normal control shown in Supplementary Figure S2.

MiR-205-5p and PTPRM were combined into a panel by logistic regression analysis, and the equation for calculating the probability of LUSC was: Logit(P) = 1.009 + 0.03 × miR-205-5p—4.767 × PTPRM. As is demonstrated in Figures 5A,B, the AUC of the panel was 0.994 (95% CI: 0.989–1.000, p < 0.0001) in TCGA-LUSC and 0.858 (95% CI: 0.746–0.951, p < 0.0001) in the external validation. The DCA showed that the miRNA-mRNA network had good diagnostic performance under arbitrary threshold probability (Figures 5C,D).

Based on the analysis of FIGO stages, there was no statistically significant difference in the expression of miR-205-5p and PTPRM in the early stage (I) and advanced stage (II + III + IV) (Supplementary Figures S3A,B). The expression of miR-205-5p in male patients was higher than that in female patients (Supplementary Figures S3C,D), and PTPRM was higher in age over 68 years (Supplementary Figures S3E,F). Univariate and multivariate cox regression analyses were used to estimate the Hazard ratio (HR) of different clinical features in the TCGA-LUSC. K-M survival analysis and cox regression analysis showed that miR-205-5p and PTPRM were not associated with prognosis (Supplementary Figure S4).

We downloaded the ssGSEA enrichment score of TCGA-LUSC data from UCSC Xena, and then analyzed the correlation between expression values of miR-205-5p and PTPRM and enrichment score. The results showed that miR-205-5p and PTPRM were highly correlated with immune-related pathways, such as yaci and bcma stimulation of b cell immune responses and Translocation of ZAP-70 to Immunological synapse (Figure 6A). Therefore, we further analyzed its correlation with immune cells to explore its role in tumor immunity. We conducted a differential analysis of the immune cell data in TCGA-LUSC and found that 13 types of immune cells differed between tumor and normal tissues listed in Supplementary Table S5. Then, the correlation analysis of differential immune cells with miR-205-5p and PTPRM were carried out in the following research. As shown in Figure 6B, the expression of miR-205-5p and PTPRM were correlated with dendritic cells activate. We used CIBERSORT to calculate the proportion of various immune cells in each TCGA-LUSC sample, ESTIMATE to estimate the stroma and immunity levels of the samples, and obtained the methylation levels of CpG sites in TCGA-LUSC samples from the UCSC Xena platform. In conclusion, miR-205-5p and PTPRM interact with DNA methylation, tumor immunity and inflammation in the tumor microenvironment. MiR-205-5p and PTPRM had no correlation with TMB (Figure 6C).