Identification of TNIK as a novel potential drug target in thyroid cancer based on protein druggability prediction.

ABSTRACT: Thyroid cancer is a common endocrine malignancy; however, surgery remains its primary treatment option. A novel targeted drug for the development and application of targeted therapy in thyroid cancer treatment remain underexplored.We obtained RNA sequence data of thyroid cancer from The Cancer Genome Atlas database and identified differentially expressed genes (DEGs). Then, we constructed co-expression network with DEGs and combined it with differentially methylation analysis to screen the key genes in thyroid cancer. PockDrug-Server, an online tool, was applied to predict the druggability of the key genes. Finally, we constructed protein-protein interaction (PPI) network to observe potential targeted drugs for thyroid cancer.We identified 3 genes correlated with altered DNA methylation level and oncogenesis of thyroid cancer. According to the druggable analysis and PPI network, we predicted TRAF2 and NCK-interacting protein kinase (TNIK) sever as the drug targeted for thyroid cancer. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis indicated that genes in protein-protein interaction network of TNIK enriched in mitogen-activated protein kinase signaling pathway. For drug repositioning, we identified a targeted drug of genes in PPI network.Our study provides a bioinformatics method for screening drug targets and provides a theoretical basis for thyroid cancer targeted therapy.

Introduction

Thyroid cancer is a common endocrine malignant tumor, which has gained widespread concern because of its fastest-growing. In the past few years, the morbidity and mortality of thyroid cancer have increased at a rate of more than 3% in the United States.[ In China, the morbidity has increased from ∼3.21/1,000,000 to 9.61/1,000,000 in the past decade.[ Despite great advancements have achieved in thyroidectomy, radioiodine therapy, and thyroid-stimulating hormone inhibition therapy, they still have not changed the increasing trend of morbidity and mortality per year. Recently, many patients can benefit from targeted drugs with the development of personalized therapy[; however, few targeted drugs have been approved for the treatment of thyroid cancer. Therefore, it is essential to explore novel drug targets for thyroid cancer treatment.

Epigenetic changes are regarded as significant contributors to tumor progression.[ DNA methylation is 1 of the most important epigenetic modifications and leads to the development and progression of thyroid cancer by altering the genic status. Ras association domain family 1 isoform A (RASSF1A) methylation is the potential molecular marker to characterize the histopathology of papillary thyroid cancer (PTC).[ While DNA methylation level of genes could be used to differentiate non-malignant tumors from thyroid cancer.[ Additionally, DNA methylation also serves as a cancer therapeutic target. Hegi et al determined that O6-methylguanine-DNA methyltransferase (MGMT) methylation is sensitive to temozolomide in patients with glioblastomas.[ Fujii et al demonstrated that patients with human epidermal growth factor 2 (HER2)-positive breast cancer with hydroxysteroid dehydrogenase type 4 (HSD17B4) methylation are sensitive to trastuzumab combined with chemotherapy.[ However, there are few studies on the relationship between DNA methylation and thyroid cancerous therapy.

Drug research and design are very complex and require a substantial amount of time and capital. Thus, shortening the time is required and reducing the cost for drug development is imperative. Drug repositioning can be used to confirm the new use of approved drugs and greatly reduces the time and cost for drug development, using the construction of protein-protein interaction (PPI) network.[ PPI network allows researchers to determine potential drugs for diseases through their interactions with known drug targets or proteins with indirect effects.[ In recent years, an increasing number of studies have focused on drug repositioning. Wang et al used weighted gene co-expression network analysis (WGCNA) and constructed a PPI network to detect 3 drug targets and 15 candidate drugs for melanoma treatment.[ Islam et al identified 238 gene signatures as therapeutic targets using a PPI network and 37 novel drugs as potential anticancer drugs for low-grade glioma.[

In the current study, we aimed to predict the drug target to improve target therapy of thyroid cancer. We analyzed the differentially expressed genes (DEGs) and obtained expression data from The Cancer Genome Atlas (TCGA) database. Then WGCNA and differentially methylation analysis were used to identify the key genes associated with thyroid cancer. Then, a drug target with druggable protein pockets was identified. Based on the drug repositioning, we constructed PPI network and predict potential drugs by approved drug targets in PPI network. Finally, we annotated the genes in the PPI network.

Materials and methods

Data and sources

The RNA sequencing data and DNA methylation data of thyroid cancer were downloaded from TCGA database (https://portal.gdc.cancer.gov). The RNA sequencing data included 510 thyroid cancerous samples and 58 normal samples, DNA methylation data included 507 thyroid cancerous samples and 56 normal samples. DEGs were calculated by “edgeR” package in R language.[ The cut-off criteria were set as |log2 Fold Change| > 1 and False discovery rate (FDR) <0.05.

Ethical approval or patient consent was not required because the data for the present research were obtained from a public database, and the data were available without personal identifiers.

Weighted gene co-expression network analysis

The “WGCNA” R package was used to construct a co-expression network to identify the hub genes in thyroid cancer.[ First, the DEGs and the Pearson correlation coefficient were used to confirm the most correlated genes and exclude weakly correlated genes. We calculated the soft thresholding power (β) using network topology analysis and converted the adjacency to the topological overlap matrix.[ Second, a gene dendrogram was constructed using hierarchical clustering. We used a dynamic tree-cut algorithm to separate the branches of the gene dendrogram from modules of co-expressed genes into different colors[; deepsplit = 2 and minimal module size = 30. Third, we estimated the similarity of the modules and merged the genes with high co-expression. We selected 0.25 as the threshold for the dissimilarity between the modules and merged highly co-expressed modules. Then, the module-traits relationship was constructed by measuring the relevance between the module eigengenes and thyroid cancer.[ The correlations indicated the gene significance and module membership with P value <0.05.

Identify differentially methylated gene

The differentially methylated genes (DMGs) between thyroid cancerous and normal samples were analyzed with t test (FDR < 0.05 and the absolute MTBeta-MNBeta >.3; MTBeta, MNBeta: β means of tumor samples and normal samples). Methylation levels were calculated as follows: we used the β value to estimate the methylation level of a given CpG probe.

PockDrug druggability and PPI analysis

PockDrug-Server (http://pockdrug.rpbs.univ-paris-diderot.fr/) is a robust pocket druggability prediction server with references pocket boundary uncertainties to query the druggable of protein pockets.[ Therefore, we used PockDrug-Server to predict drug targets with druggable protein pockets in the overlapped genes between hub genes and DMGs. A druggability probability of more than 0.5 was considered a druggable pocket.

For drug repositioning, we used the Search Tool for the Retrieval of Interacting Genes (STRING) database to construct a PPI network of drug targets with druggable protein pockets. STRING database (http://string-db.org/) is used to calculate protein-protein association for many organisms and can uncover direct and indirect relationships. The 0interaction score should be >0.7. The data were visualized by Cytoscape with version 3.3.0.[ Drug information were obtained from Drugbank (https://go.drugbank.com/).

Functional enrichment analysis

Functional enrichment analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) version 6.8 (https://david.ncifcrf.gov/), a web-accessible program that integrates functional genomic annotations with intuitive graphical summaries.[ The annotated genes carries out gene ontology (GO) analysis, elucidating biological process (BP), molecular function (MF), and cellular component (CC).[ Kyoto Encyclopedia of Genes and Genomes (KEGG, version 90.0; www.kegg.jp) is a common online resource for interpreting biological systems from molecular-level data.[ Thus, DAVID was used to annotate genes in the PPI network and GO and KEGG enrichment analyses were performed to elucidate enrichment of the genes. P < .05 was considered statistically significant for all tests.

Results

WGCNA analysis and module significance calculation

The top 5000 DEGs of thyroid cancer were obtained from the TCGA database, which was used to construct the co-expression network. β = 3 (R2 = 0.9) was set as the soft-threshold to construct the co-expression module (Fig. 1). Ten modules were identified and shown in different colors in dendrogram (Fig. 2A). The eigengene dendrogram and adjacency heatmap were provided in Figure 2C. Then, we calculated the correlation between module eigengenes and the clinical trait of thyroid cancer. As shown in Figure 2B, the red (r = 0.51; P = 3e−38) and blue module (r = 0.51; P = 5e−39) were positively correlated with the occurrence of thyroid tumor, while the green module (r = −.53; P = 2e−24) was negatively correlated with the occurrence of thyroid tumor. According to the correlation between module membership and gene significance (Fig. 2D), the red module had a high correlation (cor = 0.73; P = 9.7e−26). Thus, red module was selected for further analysis.

Figure 1

Determination of the soft-threshold power of the network topology. Analysis of the scale-free fit index for various soft-threshold powers (β) on the left. Analysis of the mean connectivity for various soft-threshold powers in the right.

Figure 2

Gene enrichment and modules identified by WGCNA. (A) Clustering dendrograms of screened genes. Different cluster modules are in different colors. (B) Heatmap of correlation between ME and TC samples. Each row represents a module eigengene, each column represents TC and normal thyroid sample, and each cell contains the corresponding correlation and P value. The matrices are color-coded by correlation according to the color legend. (C) Dendrogram tree and adjacency heatmap of the eigengenes. (D) Scatter plot of the module membership and gene significance in the red module. ME = module eigengenes, TC = thyroid cancer, WGCNA = weighted gene co-expression network analysis,

DMGs in the red module

DNA methylation data of thyroid cancer was obtained from TCGA database. Methylated genes identified by t test and FDR < 0.05 were regarded as differentially methylation. We identified 445 DMGs. Meanwhile, intersection analysis between DMGs and hub genes in the red module found 3 genes, Hairy/enhancer-of-split related with YRPW motif protein 2 (HEY2), TRAF2 and NCK-interacting protein kinase (TNIK), and Low-density lipoprotein receptor-related protein 4 (LRP4) as candidates to be further analyzed for druggability (Fig. 3A). The heatmap of these 3 genes shown in Figure 3B. These 3 genes were hypomethylated in thyroid cancer samples.

Figure 3

Screening of 3 genes from the intersection of DMGs and the red module. (A) Venn diagram between DMGs and hub genes in the red module showing the overlapping of 3 genes, HEY2, TNIK, and LRP4. (B) Heatmap of DNA methylation in HEY2, TNIK, and LRP4 between TC and normal thyroid samples. DMGs = differentially methylated genes, HEY2 = hairy/enhancer-of-split related with YRPW motif protein 2, LRP4 = low-density lipoprotein receptor-related protein 4, TC = thyroid cancer, TNIK = TRAF2 and NCK-interacting protein kinase.

Pocket druggability prediction of HEY2, TNIK, and LRP4 proteins

Druggability is the capacity of a protein to bind to drug-like molecules with high affinity.[ Therefore, it is necessary to assess druggability as the first step of drug target discovery. In our study, we used the PockDrug-Server to predict the pocket druggability of HEY2, TNIK, and LRP4 proteins, and only TNIK had 8 protein pockets (P0–P7; Fig. 4). The protein pockets with an average druggability probability >.5 were considered as druggable pockets. Table 1 showed the parameters of the 8 protein pockets, and the result indicated that P0 (0.9, P = .01), P2 (0.82, P = 0.05), P3 (0.86, P = 0.0), and P6 (0.92, P = 0.02) had the highest probability of druggability, and among them, P6 had the highest druggability probability score. Therefore, we chose TNIK as a possible drug target for thyroid cancer.

Figure 4

The protein pocket of TNIK. P0–P7 represent different protein pockets in TNIK. TNIK = TRAF2 and NCK-interacting protein kinase.

Table 1

Parameters of the 8 protein pockets in TNIK.

Pockets	Vol. Hull.∗	Hydroph. Kyte∗	Polar Res.∗	Aromatic Res.∗	Otyr atom	Nb. Res.∗	Drugg Prob∗	Standard deviation
P0	2418.61	0.22	0.47	0.13	0.01	32	0.9	0.01
P1	388.61	−1.61	0.64	0.09	0.0	11.0	0.05	0.01
P2	444.17	0.04	0.57	0.14	0.03	14.0	0.82	0.05
P3	1161.44	0.22	0.37	0.11	0.0	19.0	0.86	0.0
P4	459.98	−0.46	0.54	0.08	0.0	13.0	0.44	0.04
P5	582.49	−0.54	0.57	0.07	0.0	14.0	0.4	0.04
P6	372.2	0.4	0.4	0.2	0.0	10.0	0.92	0.02
P7	491.79	−1.94	0.89	0.22	0.0	9.0	0.04	0.01

PPI of TNIK and targeted drug analysis

We constructed a PPI network for TNIK using Cytoscape, which consisted of 11 nodes and 39 edges from the STRING database and found that 10 proteins interacted with TNIK (Fig. 5). To determine the biological functions of genes in the PPI network, we performed GO and KEGG pathway enrichment analysis. The genes enriched in the regulation of protein phosphorylation, regulation of mitogen-activated protein kinase (MAPK) cascade, β-catenin/transcription factor 7-like 2 complex, kinase binding, and protein phosphatase binding (Fig. 6). According to KEGG analysis, genes enriched in the MAPK, tumor necrosis factor (TNF) signaling pathways, and amyotrophic lateral sclerosis, among other pathways (Table 2). According to drug repositioning, we searched the targeted drug of proteins in PPI network, and found binimetinib targeted Mitogen-activated protein kinase kinase kinase 1 (MAP3K1). Based on this, we speculated that binimetinib might serve as a targeted drug for thyroid cancer in the future.

Figure 5

PPI network map of TNIK to identify targeted drug genes. TNIK = TRAF2 and NCK-interacting protein kinase. PPI = protein-protein interaction, TNIK = TRAF2 and NCK-interacting protein kinase.

Figure 6

GO enrichment of the genes in the PPI, including BP, MF, and CC. BP = biological process, CC = cellular component, GO = gene ontology, MF = molecular function, PPI = protein-protein interaction.

Table 2

KEGG pathway analysis for genes in PPI.

Pathway ID	Term description	Observed gene count	P value
hsa04010	MAPK signaling pathway	8	2.03E−11
hsa04668	TNF signaling pathway	7	5.84E−12
hsa05014	ALS	5	1.99E−09
hsa04071	Sphingolipid signaling pathway	5	8.38E−08
hsa04210	Apoptosis	5	1.30E−07
hsa04622	RIG-I-like receptor signaling pathway	4	6.48E−07
hsa04920	Adipocytokine signaling pathway	4	6.48E−07
hsa04064	NF-kappa βsignaling pathway	4	1.64E−06
hsa04664	Fc epsilon RI signaling pathway	3	3.35E−05
hsa04912	GnRH signaling pathway	3	6.69E−05
hsa04657	IL-17 signaling pathway	3	7.26E−05
hsa05216	Thyroid cancer	2	0.0007

Discussion

Thyroid cancer is a common endocrine malignancy; however, there are few targeted therapeutic drugs currently available. In this study, we screened genes based on mRNA profiles and DNA methylation level of thyroid cancer in TCGA, and used the PockDrug-Server in combination with the PPI network to identify genes that could target thyroid cancerous treatment. Ten co-expression modules were structured based on the top 5000 DEGs between thyroid cancerous and normal samples by WGCNA, and the red module was determined to have the greatest correlation with thyroid carcinogenesis. Based on interactions with the hub genes and DMGs, we screened 3 genes, HEY2, TNIK, and LRP4. Our heatmap results showed that these 3 genes were hypomethylated in thyroid cancer. DNA hypomethylation is associated with neoplastic progression.[ Therefore, HEY2, TNIK, and LRP4 may play important physiological roles in the occurrence of thyroid cancer. HEY2 is the downstream target of the Notch signaling pathway and the activated Notch signaling pathway inhibits the progression of medullary thyroid cancer.[LRP4 is up-regulated in PTC and activated phosphatidylinositol-3 kinases/serine/threonine-protein kinases-mediated mesenchymal transition promotes PTC metastasis.[

We used PockDrug-Server to analyze and confirm the druggability of protein (HEY2, TNIK, and LRP4) pocket. Protein pocket druggability predicted the affinity of protein pockets to bind drug-like molecules, and was a major criterion in the identification phase of drug target discovery.[ PockDrug-Server is a new bioinformatics online tool for efficiently measuring the druggability of holo-like and apoprotein pockets.[ Compared with other prediction tools for druggability, it synthesizes different pocket estimation methods to provide consistent druggability results.[ Several studies have shown that PockDrug-Server is a widely used online predictive tool. Hamza et al ascertained p53 Ser121 and Val122 mutations as drug targets using PockDrug-Server.[ Trigueiro-Louro et al used PockDrug-Server to predict 3 druggable pockets and 8 new potential hot spot residues in the effector domain of Non-structural 1.[ In addition, PockDrug-Server was used to predict the physicochemical properties of binding pockets of off-targets for gefitinib.[ Our result has shown that only TNIK had 4 druggable protein pockets (Table 1). The score of these 4 druggable protein pockets were all more than 0.8, indicating that they have greatly strong drug binding ability. Consequently, TNIK was predicted as a drug target for thyroid cancer.

TNIK is a member of the Sterile 20 (STE20) serine/threonine protein kinase family and severs as a driver oncogene.[TNIK was an important activator of Wnt signaling pathway to promote tumor progression and invasion in terminal colorectal cancer.[ Moreover, increasing expression of TNIK was the biomarker of poor prognosis in patients with pancreatic, and hepatocellular carcinomas.[ Recently, various classes of small-molecule TNIK inhibitors have developed including NCB-0846 and KY-05009. Some studies have reported that NCB-0846 has high inhibition against TNIK to block Wnt signaling pathway and achieves the purpose to anti-tumor on colorectal cancerous cells.[ While KY-05009 is also a small-molecule of TNIK inhibitor, it suppresses the Transforming growth factor-β1-induced activation of Wnt signaling pathway and epithelial-mesenchymal transition process in lung adenocarcinoma cells.[ However, the potential of these small-molecule TNIK inhibitors has remained unexplored. Thus, it remains to explore drugs targeted to TNIK. In our study, based on the drug repositioning, we constructed the PPI network for TNIK and predicted the underlying targeted drug for thyroid cancer.

Ten protein interacted with TNIK was in the PPI network. KEGG pathway enrichment indicated that most genes were concentrated in the MAPK and the TNF signaling pathways. The MAPK and TNF signaling pathways are common pathogenic signaling networks that influence cellular proliferation, differentiation, development, inflammatory responses, and apoptosis.[ Some small inhibitors of MAPKs signaling pathway effectors can affect the growth of thyroid cancerous cells. Of these 10 genes, TNF, tumor necrosis factor receptor superfamily member 1A (TNFRSF1A), tumor necrosis factor receptor-associated factor 2 (TRAF2), and tumor necrosis factor receptor type 1-associated DEATH domain protein (TRADD) are key effector molecules in the TNF signaling pathway. TRAF2 is able to suppress cell death induced by TNFα, and as death receptors or the target of anti-cancer drugs.[MAP3K1, dual specificity mitogen-activated protein kinase kinase 6 (MAP2K6), dual specificity mitogen-activated protein kinase kinase 3 (MAP2K3), and mitogen-activated protein kinase kinase kinase 5 (MAP3K5) are the effectors of the MAPK signaling pathway.[MAP3K1 is a member of the MAP3K family and the STE superfamily and induces multiple signaling pathways including Wnt/β-catenin signaling pathway and nuclear factor-kappa B pathway to regulate cell suvival and apoptosis.[ While mutation and copy number loss of MAP3K1 were observed in ovarian, prostate, and breast cancer.[

In recent years, targeted therapy has become an attractive strategy for cancer therapy. However, surgery is still preferred for thyroid cancerous treatment. As drug research is a long complex process that requires a great deal of investment, targeted therapies for thyroid cancer have remained underexplored. Drug repositioning is a good selection step to overcome the limitations of traditional methods. It can find novel uses for existing drugs through PPI network to reduce the costs and risks of drug development.[ In our study, based on the PPI network of TNIK, we identified the targeted drugs searching DrugBank and found that binimetinib targeted MAP3K1 which is a potent and selective oral Mitogen-activated protein kinase 1/2 inhibitor. The US Food and Drug Administration approved binimetinib combined with encorafenib for patients with unresectable or metastatic melanoma with Serine/threonine-protein kinase B-rafV600E(BRAFV600E) or BRAFV600K mutation on June 27, 2018. Clinical studies on binimetinib are currently underway for the treatment of thyroid cancer with BRAFV600E mutation (NCT04061980). Therefore, binimetinib is likely to have potential effects in the treatment of thyroid cancer and we hope that the drug would achieve good results in clinical trials.

Owing to the restriction of information collection regarding detailed interaction dynamics and PPI network, the prediction results are not comprehensive. Clinical trials on binimetinib for thyroid cancerous treatment are being performed but its efficacy and safety still require experimentation for verification. Our predicted results provide a theoretical basis for future targeted therapy of thyroid cancer.

Conclusion

In conclusion, with WGCNA and druggability analysis for genes in thyroid cancer, we revealed that TNIK could serve as a potential drug target for the treatment of thyroid cancer. Further studies that explore novel methods to screen targeted genes with druggability and potential drug targets for thyroid cancer are imperative. We hope these findings will contribute to the research on new-targeted therapeutic drugs for thyroid cancer.

Acknowledgments

We thank Omigen, Inc. for their assistance in data analysis.

Author contributions

YFY and BY conceived and designed the study; BY, XXZ, and YHZ carried out data collection and performed data analysis; YFY wrote the paper; BY, XXZ, and YHZ edited the manuscript and provided critical comments. All authors read and approved the final manuscript.

Conceptualization: YiFei Yang, Bin Yu.

Data curation: Bin Yu, Xiuxia Zhang, Yunhua Zhu.

Formal analysis: Bin Yu, Xiuxia Zhang, Yunhua Zhu.

Methodology: Bin Yu.

Visualization: Bin Yu, Xiuxia Zhang, Yunhua Zhu.

Writing – original draft: YiFei Yang, Xiuxia Zhang.

Writing – review & editing: YiFei Yang, Bin Yu, Yunhua Zhu.