BACKGROUND: Although with a good prognosis of papillary thyroid cancer (PTC), patients with PTC and also experiencing lymph node metastasis (LNM) had higher recurrence and mortality rates. Therefore it was essential to explore novel biomarkers or methods to predict and evaluate the situation in the stages of PTC. METHODS: In this study, mRNA sequence datasets from The Cancer Genome Atlas (TCGA) database and Gene Expression Omnibus (GEO) were utilized to obtain differentially expressed genes (DEGs) between PTC tumors and normal specimens and DEGs related to lymph node metastasis were identified using weighted gene co-expression network analysis (WGCNA) according to the clinical information. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were applied to query the biological functions and pathways. Furthermore, a protein-protein interaction (PPI) network was constructed using a STRING database and a prognosis model was established using the least absolute shrinkage and selection operator (LASSO) Cox regression analysis based on the LNM-related DEGs. Finally, six hub genes were identified and verified in vitro experiments. RESULTS: A novel six-gene signature model including COL8A2, MET, FN1, MPZL2, PDLIM4 and CLDN10 was established based on a total of 52 DEGs from the intersection of LNM-related modules identified by WGNCA from TCGA, THCA and GSE60542 to predict the situation of lymph node metastasis in PTC. Those six hub genes were all more highly expressed in PTC tumors and played potential biological functions on the development of PTC in in vitro experiments, which had potential values as diagnostic and therapeutic targets.