{'fraction_non_alphanumeric': 0.040190350546551075, 'fraction_numerical': 0.0354787033546928, 'mean_word_length': 5.11671469740634, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 9, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of potential pathogenic candidates or diagnostic biomarkers in papillary thyroid carcinoma using expression and methylation profiles

The mechanisms underlying the pathogenesis of papillary thyroid carcinoma (PTC) have not yet been elucidated. The aim of the current study was to identify potential pathogenic biomarkers in PTC by comprehensively analyzing gene expression and methylation profiles, and to increase the understanding of PTC pathogenesis. The gene expression profiles of the GSE97001 and GSE83520 datasets, the miRNA expression profiles of the GSE73182 dataset, and the DNA methylation profiles of the GSE86961 and GSE97466 datasets were downloaded from Gene Expression Omnibus database. The differentially expressed genes (DEGs) and the differentially expressed microRNAs (DEMs) were identified using the limma package in R, and the differentially methylated sites (DMSs) were identified using the β distribution and two-sample t-tests. The Database for Annotation, Visualization and Integrated Discovery, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome were subsequently used to perform functional and pathway enrichment analysis. The miRNA target genes were predicted using the online databases miRWalk. The protein-protein interactions (PPI) were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins. The regulatory network was constructed, and the gene expression and methylation levels of the key nodes were detected using reverse-transcription quantitative-polymerase chain reaction (PCR) and methylation-specific PCR. A total of 155 overlapping DEGs were identified between the GSE97001 and GSE83520 datasets, and 19 DEMs between PTC tissue and normal tissue samples were identified in the GSE73182 set. In the GSE86961 and GSE97466 datasets, 2,910 overlapping DMSs that were associated with 38 downregulated methylated genes were identified. The overlapping DEGs were enriched in 46 Gene Ontology terms and one KEGG pathway. A total of 60 PPI pairs were identified for the overlapping DEGs and 12 negative miRNA-gene pairs were identified for the DEMs. The expression levels of hsa-miR-199a-5p and decorin (DCN) were decreased in patients with PTC. C-X-C motif chemokine ligand 12 (CXCL12) was hypermethylated and had a decreased expression level in PTC tissues. LDL receptor related protein 4 (LRP4) and carbonic anhydrase 12 (CA12) were hypomethylated and had an increased expression level. The present study revealed that hsa-miR-199a-5p, DCN, CXCL12, LRP4 and CA12 may serve important roles in the pathogenesis of PTC.


Introduction
Thyroid cancer is one of the fastest growing malignancies and it was the fifth most common cancer worldwide in 2017 (1). The morbidity of thyroid cancer has increased significantly over the past 30 years, with a 300% increase in the USA and an increase from ~10/1,000,000 to 30-40/1,000,000 in China (2). Papillary thyroid carcinoma (PTC) is an endocrine malignancy and the most common type of thyroid cancer, accounting for ~80% of all thyroid cancer cases worldwide (3). Late diagnosis may be the primary cause of the increase in PTC morbidity (4). Therefore, early diagnosis and treatment are required to improve patient prognosis and increase survival rates (5). In recent years, there has been an increased understanding of the pathogenesis underlying PTC, particularly in microRNAs (miRNAs/miRs) and DNA methylation associated with the occurrence, progression, metastasis and prognosis of the disease. A previous study investigating 653 patients with PTC revealed that mutations in the promoter region of telomerase reverse transcriptase and the B-Raf proto-oncogene serine/threonine kinase (BRAF) V600E mutation contributed to the aggressiveness of PTC, particularly when the two mutations coexisted (6). Another study demonstrated that coiled-coil domain containing 67 (CCDC67) was downregulated in PTC, and revealed an inverse relationship between CCDC67 expression and PTC aggressiveness and BRAF T1799A mutation and that CCDC67 expression may inhibit PTC cells (7). miRNAs regulate gene expression by binding to the 3'untranslated regions of target mRNAs (8). Cheng et al (9) revealed that miR-150 inhibits cell proliferation, migration and invasion of PTC by targeting ρ-associated protein kinase 1 (ROCK1), and suggested that miR-150 may be an effective therapeutic target for the treatment of PTC. Liu et al (10) revealed that miR-363-3p inhibited the progression of PTC by targeting phosphoinositide-3-kinase catalytic subunit α. The promoter methylation status of genes associated with thyroid cancer including caspase 8, cyclin dependent kinase inhibitor 2A, death associated protein kinase 1 (DAPK1), estrogen receptor 1 (ESR1), solute carrier family 5 member 5, Ras association domain family member 1 (RASSF1) and TIMP metallopeptidase inhibitor 3, and the abnormal methylation of RASSF1, DAPK1 and ESR1 may be used as molecular markers for the early detection of thyroid cancer and the diagnosis of the follicular thyroid carcinoma (FTC)-Hurthle and FTC-Classic subtypes (11). The current study investigated differentially expressed genes (DEGs), differentially expressed miRNAs (DEMs) and differentially methylated sites (DMSs) in PTC samples. Functional enrichment analysis of DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome, and the miRNA-target interactions of DMSs and protein-protein interactions (PPI) of DEGs were identified separately using the miRWalk database and the Search Tool for the Retrieval of Interacting Genes/Proteins. The regulated network was conducted based on the aforementioned miRNA-target interactions and PPI to identify potential biomarkers of PTC. The potential biomarkers were further validated using clinical PTC samples using reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and methylation-specific PCR (MSP). The aim of the current study was to identify potential pathogenic biomarkers in PTC, and to increase the understanding of PTC pathogenesis.

Materials and methods
Expression and methylation profiles. The gene expression datasets GSE97001 (12) and GSE83520 (13), the miRNA expression profile dataset GSE73182 (14), and the DNA methylation profile datasets GSE86961 (15) and GSE97466 (16) were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo). A total of four PTC tissue and four normal thyroid tissue samples were included in the GSE97001 dataset and were analyzed using an Agilent-026652 Whole Human Genome Microarray 4x44K V2 platform. A total of 12 PTC tissue and 12 normal tissue samples in the GSE83520 dataset were analyzed using an Illumina HiSeq 2500 (Homo sapiens) platform, subsequently the data produced were annotated using the University of California, Santa Cruz (UCSC) Genome Browser (genome.ucsc.edu) and GENCODE (www.gencodegenes.org) respectively using UCSC.txt.gz and GENCODE.txt.gz files. The GSE73182 dataset contained 19 PTC tissue and five normal tissue samples, which were analyzed using an Agilent-035758 Human miRBase 16.0 plus 031181 platform. A total of 41 PTC tissue and 41 normal tissue samples in the GSE86961 dataset and 60 PTC tissue and 50 normal tissue samples in the GSE97466 dataset were analyzed using an Illumina HumanMethylation450 BeadChip (HumanMethylation450_15017482) platform. Data processing and differential analysis. The raw microarray data were converted into expression data and normalized using the preprocessCore function package (version 3.5; www. bioconductor.org/packages/release/bioc/html/preprocessCore. html). TopHat, an application within preprocessCore was used to map the RNA-Sequencing reads to the GRCh37/hg19 genome with the maximum mismatches of two. The sorted .bam files from TopHat were used to calculate the numbers of the reads which were mapped to the exons of a gene based on Cufflinks, another application within preprocessCore and normalized to reads per kilobase of the exon model per million mapped reads (RPKM), which is the representation of the expression values of the genes. Finally, through Cuffdiff in Cufflinks DEGs were obtained using the cut-offs |log2(RPKM Normal/RPKM PTC)|>1 and P<0.05. The average values of multiple probes that corresponded to the same gene were determined. The DEGs and DEMs between PTC tissue and normal samples were screened using the limma package (version 3.18.13; www.bioconductor.org/packages/2.13/bioc/html/limma.html). P<0.05 and |log fold-change (FC)|>2 were used to identify DEGs, and P<0.05 and |log FC|>1 were used to identify DEMs. The gplots package (version 3.0.1; cran.r-project.org/web/packages/gplots) was used to conduct the cluster analysis.
For original DNA methylation data, from the datasets GSE86961 and GSE97466, the β value of every methylated site was calculated and normalized with the β distribution test and compared using a two-sample Student's t-test. The DMSs were identified by comparing β values between cases and controls using |Δβ|>0.2 and P<0.05 as the cut-off criteria.
Functional and pathway enrichment analysis. DAVID (version 6.8; david.ncifcrf.gov) is a web-accessible program that integrates functional genomic annotations with intuitive graphical summaries (17). KEGG (www.genome.jp/kegg) is an online resource for interpreting high-level functions of biological systems from molecular level data, including complete and partial sequenced genome sequences, graphical cellular biochemical processes, and information on chemicals, enzymes, molecules, and enzymatic reactions (18). Reactome (www.reactome.org) is a free, open-source, curated and peer-reviewed pathway database, which provides intuitive bioinformatics tools for the visualization, interpretation and analysis of pathway knowledge to support basic research, genome analysis, modeling, systems biology and education (19). Gene Ontology (GO), KEGG and Reactome were used to perform enrichment analysis of the DEGs in DAVID.
Target prediction of DEMs. The targets of the DEMs were predicted using the miRWalk database (version 2.0; www. umm.uni-heidelberg.de/apps/zmf/mirwalk), a comprehensive atlas of predicted and validated miRNA-target interactions. Putative targets were predicted using the ten bioinformatics databases in the miRWalk database (http://www.umm. uni-heidelberg.de/apps/zmf/mirwalk/), in which the miRNA-target interactions proved by at least one of the ten bioinformatics databases were selected. Subsequently, the inverse miRNA-gene regulation pairs were screened out and the locations of the targets were selected.
Establishing the regulated network of DEGs and DEMs. The PPIs were analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (version 10.0; string-db.org), and the DEG PPI pairs with a score >500 were screened out. Based on inverse/opposite miRNA-gene pairs and PPI pairs (which was miRNA-gene pairs and PPI pairs with inverse/opposite regulation levels), the regulated network of the DEGs and DEMs was constructed and visualized using the Cytoscape software (version 3.5.1; www.cytoscape.org/download.php).

Investigation of the potential biomarkers in patients with PTC.
A total of 10 patients with PTC (6 females and 4 males) with a mean age of 40.6 years (range, 32-65 years) were admitted to Tianjin Medical University Cancer Institute and Hospital between July 2017 and December 2017. The study was approved by the Tianjin Medical University Cancer Institute and Hospital Ethics Committee (no. bc2018071), and all patients provided written informed consent. PTC tissues as well as paired paracancerous (1.0 cm from the tumor tissue) thyroid specimens were collected from the patients. Total RNA from the tissue samples was extracted using TRIzol (Invitrogen; Thermo Fisher Scientific, Inc.). RT-qPCR and MSP were performed to detect the mRNA level and methylation status of hsa-miR-199a-5p, decorin (DCN), C-X-C motif chemokine ligand 12 (CXCL12), LDL receptor related protein 4 (LRP4) and carbonic anhydrase 12 (CA12). The primers were designed and synthesized by Takara Biotechnology Co., Ltd. and are presented in Table I. RNA was reverse transcribed into cDNA using the PrimeScript ® 1st Stand cDNA Synthesis kit (Takara Biotechnology Co., Ltd.) using the following conditions: 30˚C for 10 min, 42˚C for 60 min and 95˚C for 5 min. qPCR was subsequently performed using the SYBR ® Premix Ex Taq™ kit (Takara Biotechnology Co., Ltd.) according to manufacturer's protocol, and the QuantStudio™ 5 Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). DNA methylation modification was performed using an EZ-DNA Methylation-Gold kit (Zymo Research Corp.), according to the manufacturer's protocol. Briefly, 20 µl DNA sample was added in a PCR tube to elute the DNA using M-Elution Buffer. Afterwards, bisulfite conversion was performed as follows: ~2 µg of DNA was diluted with double distilled water to 500 µl, and incubated at 42˚C in a water bath for 30 min, 30 µl 10 mM para-dihydroxybenzene and 520 µl 3.6 M bisulfite (pH 5.0) were subsequently added to the DNA solution, and finally incubated for 16 h at 50˚C in the dark. qPCR optimization was performed for the conversed DNA, using the following thermocycling conditions: 95˚C for 5 min; followed by 40 cycles at 95˚C for  Table II, according to the logFC value. Among the DEGs-overlap, DCN and CXCL12 had a decreased expression level in PTC tissues compared with that in normal tissues, and LRP4 and CA12 had an increased expression level in PTC tissues compared with that in normal tissues. A total of 19 DEMs (10 upregulated and 9 downregulated) were identified in the GSE73182 dataset (Table III). Among them, hsa-miR-199b-5p had a lower expression level in PTC tissues compared with that in normal tissues, and was the fourth most significant DEM according to its P-value (P=1.439x10 -5 ). hsa-miR-199a-5p had a lower expression level in PTC tissues compared with that in normal tissues. hsa-miR-199a-5p is highly homologous to hsa-miR-199b-5p, and hsa-miR-199a-5p was selected for further experimentation in the current study.
A total of 3,264 (230 upregulated and 3,034 downregulated) and 2,993 (240 upregulated and 2,753 downregulated) DMSs were screened out in the GSE86961 and GSE97466 datasets, respectively. There were 2,910 overlapping DMSs (212 upregulated and 2,698 downregulated) between the two datasets. The overlapping DMSs were termed DMSs-overlap. The locations of the DMSs-overlap are presented in Fig. 1. The DMSs-overlap were primarily located in the gene body (49%) and intragenic gene region (21%). The DMSs-overlap were located in 84 differentially methylation genes (DMGs), in which 38 DMGs expression was negatively regulated by the DMSs-overlap (Fig. 2). Phospholipase A2 receptor 1, solute carrier family 26 member 4 and thyroid peroxidase were downregulated and hypermethylated in PTC tissues compared with that in normal tissues, and the others were upregulated and hypomethylated.
GO and KEGG pathway enrichment. The DEGs-overlap were enriched in 46 GO terms and one KEGG pathway. The top 20 most significant GO terms are presented in Fig. 3. The biological process 'negative regulation of catabolic process' and the cellular component 'plasma membrane' were the top two most significant terms. The KEGG enriched pathway was 'cell adhesion molecules' and the associated DEGs were claudin 16, CLDN1, cadherin (CDH)2, CDH3 and syndecan 4.

miRNA-gene pairs, PPI pairs and the regulated network.
A total of 15 miRNA-gene pairs were identified for the DEMs, 12 of which were negative regulation pairs. A total of 60 PPI pairs were obtained for the DEGs-overlap. The regulated network based on the PPI pairs and the aforementioned negative miRNA-gene pairs is presented in Fig. 4. In the network, hsa-miR-199a-5p, DCN and CXCL12 were nodes with the highest scores, and CXCL12, CA12, DCN, KIT proto-oncogene, receptor tyrosine kinase and apolipoprotein E were the target genes of hsa-miR-199a-5p, predicted using the miRWalk database (version 2.0). LRP4 and CA12 overlapped between the regulated network and the 38 negative DMGs presented in Fig. 2. hsa-miR-199a-5p, DCN, CXCL12, LRP4 and CA12 were subsequently chosen for further experimentation in patients with PTC.
Expression and methylation levels of the potential biomarkers. The expression and methylation levels of the aforementioned potential biomarkers are presented in Fig. 5 and Table IV. The expression levels of hsa-miR-199a-5p and DCN were significantly decreased in PTC tissues compared     with that in normal paracancerous tissues (P<0.001); however, no statistically significant difference in the methylation levels were observed (P>0.05). CXCL12 was significantly downregulated and hypermethylated in PTC tissues compared with that in normal paracancerous tissues (P<0.001), and LRP4 and CA12 were significantly upregulated and hypomethylated in PTC tissues compared with that in normal paracancerous tissues (P<0.001).

Discussion
The present study investigated DNA methylation and gene expression profiles and identified hsa-miR-199a-5p, DCN, CXCL12, LRP4 and CA12 as potential biomarkers in PTC.
The expression levels of hsa-miR-199a-5p and DCN were decreased in patients with PTC. CXCL12 was downregulated and hypermethylated in PTC tissues, while LRP4 and CA12 were upregulated and hypomethylated. hsa-miR-199a-5p has been extensively studied in several types of cancer such as lung cancer (21) and breast cancer (22).
However, there are few reports on the role of hsa-miR-199b-5p in PTC. A recent study revealed that hsa-miR-199a-5p expression was significantly downregulated in PTC tissues and cells, and it inhibited the progression of PTC by downregulating snail family transcriptional repressor 1 (SNAI1) (23). A previous study demonstrated that hsa-miR-199a-5p was downregulated in tumor tissues obtained from patients with follicular thyroid carcinoma (FTC), and exhibited a negative association with its target, connective tissue growth factor (24). Therefore, hsa-miR-199a-5p may serve as a novel biomarker for the clinical management of thyroid carcinoma. Similar to previously published studies, the differential expression analysis and RT-qPCR performed in the current study revealed that the expression of hsa-miR-199a-5p was decreased in patients with PTC. The present study revealed that CXCL12, CA12, DCN, KIT proto-oncogene, receptor tyrosine kinase and apolipoprotein E are potential target genes of hsa-miR-199a-5p, which may aid to elucidate the pathways involved in the pathogenesis of PTC. Chung et al (25) reported that CXCL12 has a higher sensitivity and specificity as a novel diagnostic marker for PTC compared with cytokeratin 19, Hector Battifora mesothelial epitope-1 and galectin-3. Zhang et al (26) conducted a reduced representation bisulfite sequencing analysis and identified CXCL12 as an important gene in the co-function network, which suggested that CXCL12 may facilitate PTC progression by methylation-mediated epigenetic regulation of gene expression. The present study identified CXCL12 as a key pathogenic gene in the regulated network, and RT-qPCR and MSP results demonstrated that CXCL12 was downregulated and hypermethylated in PTC. Few studies have investigated the functions of DCN and LRP4 in PTC. Arnaldi et al (27) revealed that DCN was downregulated in a high percentage of patients with FTC. Jarzab et al (28) demonstrated that LRP4 was downregulated in PTC tissues, and was variably associated with the immune response of PTC (Spearman's correlation coefficient, 0.792). While data regarding the potential functions of CA12 in PTC  are limited, previous studies have demonstrated that CA12 is involved in the progression and prognosis of patients with breast and cervical cancer (29,30). The present study revealed that DCN was downregulated in PTC tissues, while LRP4 and CA12 were upregulated and hypomethylated. DCN, LRP4 and CA12 may therefore be closely associated with the pathogenesis of PTC through the methylation-mediated epigenetic regulation of gene expression. Ma et al (23) revealed that hsa-miR-199a-5p inhibited the progression of PTC by down-regulating SNAI1, as demonstrated by the RT-qPCR, cell migration, invasion and epithelial-mesenchymal transition results. According to the present bioinformatics analysis and experimental results, it can be hypothesized that hsa-miR-199a-5p may downregulate LRP4 and CA12, and upregulate DCN and CXCL12. However, the clear regulatory mechanisms need to be further verified in subsequent studies. The present study involved analysis of mRNA, miRNAs and gene methylation, and suggests that hsa-miR-199a-5p may be involved in PTC through the methylation-mediated epigenetic regulation of the aforementioned genes.
In conclusion, the present study revealed that hsa-miR-199a-5p, DCN, CXCL12, LRP4 and CA12 may be associated with the pathogenesis of PTC and may serve as potential diagnostic markers for the disease. However, future clinical studies are required to verify the results obtained in the current study.