- 1Department of Thyroid Surgery, The First Affiliated Hospital of Zhengzhou University, Zhengzhou, China
- 2Academy of Medical Science, Zhengzhou University, Zhengzhou, China
- 3School of Computer and Artificial Intelligence, Zhengzhou University, Zhengzhou, China
Background: The incidence of thyroid carcinoma (THCA), the most common endocrine tumor, is continuously increasing worldwide. Although the overall prognosis of THCA is good, patients with distant metastases exhibit a mortality rate of 5-20%.
Methods: To improve the diagnosis and overall prognosis of patients with thyroid cancer, we screened specific candidate neoantigen genes in early- and late-stage THCA by analyzing the transcriptome and somatic cell mutations in this study.
Results: The top five early-stage neoantigen-related genes (NRGs) were G protein-coupled receptor 4 [GPR4], chondroitin sulfate proteoglycan 4 [CSPG4], teneurin transmembrane protein 1 [TENM1], protein S 1 [PROS1], and thymidine kinase 1 [TK1], whereas the top five late-stage NRGs were cadherin 6 [CDH6], semaphorin 6B [SEMA6B], dysferlin [DYSF], xenotropic and polytropic retrovirus receptor 1 [XPR1], and ABR activator of RhoGEF and GTPase [ABR]. Subsequently, we used machine learning models to verify their ability to screen NRGs and analyze the correlations among NRGs, immune cell types, and immune checkpoint regulators. The use of candidate antigen genes resulted in a better diagnostic model (the area under the curve [AUC] value of the early-stage group [0.979] was higher than that of the late-stage group [0.959]). Then, a prognostic model was constructed to predict NRG survival, and the 1-, 3- and 5-year AUC values were 0.83, 0.87, and 0.86, respectively, which were closely related to different immune cell types. Comparison of the expression trends and mutation frequencies of NRGs in multiple tumors revealed their potential for the development of broad spectrum therapeutic drugs.
Conclusion: In conclusion, the candidate NRGs identified in this study could potentially be used as therapeutic targets and diagnostic biomarkers for the development of novel broad spectrum therapeutic agents.
1 Introduction
Thyroid carcinoma (THCA) begins in the thyroid gland. The main types of thyroid cancer include differentiated thyroid cancer (papillary, follicular, and Hürthle cells), medullary cancer, and anaplastic thyroid cancer (aggressive cancer) (1). Most thyroid cancers are differentiated thyroid cancers (DTC) (2) that develop from thyroid follicular cells (3). Thyroid cancer can be diagnosed at an early stage. Most early-stage thyroid cancers are diagnosed when neck lumps or nodules are noticed in the patient (4). In some cases, early thyroid cancers are detected when individuals undergo imaging tests, such as ultrasound or computed tomography scans, for other health problems (5). The prognosis of patients with thyroid cancer is better than that of patients with other cancer types. The 5-year relative survival rate of patients with localized or regional THCA is >90%, whereas that of patients with distant THCA varies according to the THCA type (6, 7).
The neoantigens of tumor-specific mutated genes are recognized by T cells and participate in the immune response of tumor cells (8–11). The antigenicity and immunogenicity of tumors depend on T cell immunoselection, in which tumor cells with strong tumor-specific mutant antigens are eliminated and those with weak (possibly mutated tumor antigens) or no (tumor cells mutated during antigen processing or presentation) antigens survive (10). Immunotherapies that enhance the ability of endogenous T cells to destroy cancer cells have demonstrated therapeutic efficacy in various malignancies, leading to the widespread application of neoantigens in clinical settings (10, 12). However, the clinical relevance of T cells in killing tumor cells remains unclear. Moreover, whether newly produced tumor-specific antigens play a protective or destructive role remains unknown (13). Recent technological innovations have facilitated the detection of tumor-specific mutations using whole exome sequencing (14, 15). To screen candidate neoantigens for tumor-specific mutations, the expression of host genes must be evaluated. Only recurrent mutations that are highly expressed in tumor cells and lowly expressed or even silent in normal cells can be considered candidate neoantigens (16).
In this study, we analyzed the functional, gene, and mutational differences between patients with early- and late-stage thyroid cancer versus normal thyroid cancer. A machine learning model was applied for further selection of features for diagnosis. The top five mutated genes were identified at each stage as potential neoantigens using the RFE method, and the prognostic characteristics of the selected 10 neoantigen-related genes (NRGs) were analyzed using a risk prediction model. The identified NRGs were effective in diagnosis and prognosis assessment. Therefore, NRGs may improve diagnostic accuracy and facilitate the development of novel targeted immunotherapy approaches to improve the outcomes of patients with THCA.
2 Materials and methods
2.1 Data preparation
We downloaded both whole exome sequencing and RNA-seq data of THCA from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Gene expression profiles (Level 3 data were downloaded from the TCGA data coordination center. This dataset showed the gene-level transcription estimates as log2(x+1) transformed RSEM-normalized counts. Clinical information of all patients, including the MNT stage and survival data, was also retrieved. Clinical information was obtained from 506 patients with THCA, of which 572 and 486 underwent transcriptome and exome sequencing, respectively.
2.2 Differentially expressed gene analysis
For RNA-seq analysis, we used 572 samples, including 59 normal tissue samples. Limma (17) (version 3.48.3 for R 4.2.3) was used to select the DEGs in early- and late-stage THCA compared to those in normal tissues. All genes with p< 0.05 and logFC beyond the 95% confidence interval were considered to be differentially expressed. We used the UpSetR package (version 1.4.0) to capture the overlapping genes of each group. PCA was used to analyze the distribution of each group. The Pheatmap package (version 1.0.12) was used to determine the expression of genes in the tumor (early and late stages) and normal groups. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to analyze the relationship between differential genes and tumor pathways in the normal versus tumor group (early and late stages) and early versus late stage.
2.3 Cluster and function analyses
Genes that were overexpressed at any stage were selected as potential neoantigens. These genes are expressed at low levels or are silent in normal tissues, which makes the development of tumor-specific antigens safe and harmless. Visual analysis of the mutant profiles in patients with early and advanced thyroid cancer was performed using the maftools package (18) (version 2.8.05 for R 4.2.3). Furthermore, differential genes with high mutation frequencies were visualized in the early versus late group, and the mutation sites were integrated for candidate overexpressed genes. We mainly use corresponding R packages such as ggplot2 (version 3.4.2) to draw diagrams.
2.4 Screening of candidate neoantigen related genes
Neoantigens should only be present in patients and expressed at low levels in normal samples. Therefore, we downloaded neoantigens from The Cancer Immunome Atlas database (https://tcia.at/). The mutant samples were intersected with the neoantigen fragment genes and the differentially expressed upregulated genes were obtained using TCIA, and candidate NRGs were initially identified. Preliminary candidates for neoantigens in the early- and late-stage groups were identified using the Search Tool for the Retrieval of Interacting Genes/Proteins (https://string-db.org/) and protein–protein interaction (PPI) networks (19). We used Pearson correlations to analyze the associations between NRGs and constructed PPI networks using Cytoscape (version 3.8.2) (20). Boxplot (version 1.3.1) and other R packages are used to generate figures.
2.5 Construction of a diagnostic model
Candidate neoantigens are only overexpressed in the early or late stages of THCA; therefore, they can be considered potential therapeutic targets. To verify the feasibility of the preliminary screening of candidate neoantigens for the diagnostic model, we used caret package (version 6.0-94 for R 4.2.3, https://cran.r-project.org/web/packages/caret/index.html) with cross validation to screen feature combinations for NRG (21). mLR (nnet, version 7.3-18, https://cran.r-project.org/web/packages/nnet/index.html), Dtree (rpart, version 4.1.19, https://cran.r-project.org/web/packages/rpart/index.html), RF (randomForest, version 4.7-1.1, https://cran.r-project.org/web/packages/randomForest/index.html), and SVM (e1071, version 1.7-13, https://cran.r-project.org/web/packages/e1071/index.html), were used to construct four machine learning models. The area under the receiver operating characteristic (ROC) curve and 10 cross-prediction ROC curves were used to evaluate the feasibility of these NRGs. TSNAdb ((http://biopharm.zju.edu.cn/tsnadb/browse/)) used to analyze the mutant protein polypeptides to estimate their binding affinity for HLA alleles.
2.6 Immunoinfiltration analysis of NRGs
CIBERSORT (22), a calculation method used to quantify the immune cell fraction from RNA-seq data, was used to calculate the immune cell infiltration score in THCA. We analyzed the differences in early- and late-stage NRG immune ratios according to the immune cell type and determined the potential correlations between NRG expression and different types of infiltrating immune cells. Finally, we collected inhibitory immune checkpoints (chemokines, receptors, MHC, and immunostimulators) with therapeutic potential from a previous study (23)and determined their potential correlation with NRGs. Plots were generated with boxplot (version 1.3.1) and ggcorrplot (version 0.1.4) packages.
2.7 Construction of a prognostic model for neoantigen-associated genes
As host genes carrying these mutations are all overexpressed in patients with THCA, NRGs serve as potential diagnostic biomarkers or therapeutic targets. Next, we assessed the impact of these genes on patient outcomes. Based on the 10 gene signatures, the NRG score was calculated using the following formula: , where exp indicates the expression level of NRGs. The median risk score was used to divide the patients into high- and low-risk groups.
2.8 Survival analysis
Kaplan–Meier survival analysis was performed using the survival (24) (version 3.5-5 for R 4.2.3). and survminer (25) (version 0.4.9) R packages. ROC curve was used to evaluate the prediction efficiency of genes for 1-, 3-, and 5-year survival with the R package, timeROC (26) (version 0.4), and log-rank test was used for comparison between groups. Wilcox test was used to compare the statistical differences between two groups. p< 0.05 was considered to be statistically significant. R, version 4.2.3 (http://www.r-project.org) is the main software environment for our statistical operation and graphics.
2.9 Comparison of candidate NRG expression levels
To investigate whether the identified NRGs were specific to THCA or common to other types of cancer, we used the GSCALite (27) (http://bioinfo.life.hust.edu.cn/web/GSCALite/) to compare the corresponding expression patterns in other cancer datasets from TCGA database (27). Comparison of NRG expression levels across tumors revealed the potential of shared neoantigens to act as broad spectrum therapeutic agents.
3 Results
3.1 Sample demographic statistics
Basic clinical information, including sex, age, TNM stage, and vital status, of the patients was shown in Table 1. The number of patients older than 55 years with late-stage THCA was higher than of patients in the early-stage (p< 0.001). More T1 stage patients were in the early stage (p< 0.001), and most patients were alive (p = 0.001). However, no significant difference in sex or TNM stage was observed between the early- and late-stage THCA groups (p = 0.093, 0.097, and 0.978, respectively; Table 1).
Table 1 Overview of thyroid carcinoma (THCA) early- and late-stage clinical data from The Cancer Genome Atlas (TCGA) database.
3.2 Distribution of DEGs
Differential analysis was conducted between the normal and other subgroups (tumor group, early, and late stages) and tumors in the early and late stages. As shown in Figure 1A, 1957 and 2112 DEGs achieved based on 2-fold change from normal to early and late stages, respectively. Upregulated and downregulated genes tended to be balanced when comparing tumor patients with normal controls. However, in the comparison of the early and late stages, DEGs were predominantly upregulated. In total, 979 upregulated and 951 downregulated genes were identified in the early-stage group and 978 upregulated and 1161 downregulated genes were identified in the late-stage group. Moreover, 84 genes were upregulated and 660 genes were downregulated in the early stage compared to the late stage. However, in the normal group, compared to the tumor group, the upregulated and downregulated genes were balanced (971,1018; Supplementary Figure 1A). Further analysis of the upset plot revealed 48 genes shared between normal controls and patients with tumors (early stage, late stage, and tumor group) and early versus late stage (Figure 1B). As shown in Figure 1C, the tumor and normal samples were separated into different groups based on the DEGs. There are also some distinctions between the early and late stages. The heatmap was used to visualize gene expression patterns across all samples but could not clearly distinguish between early- and late-stage samples. The heatmap result was consistent with the PCA finding that tumor samples showed diverse patterns compared to normal samples (Figure 1D). KEGG enrichment analysis of DEGs was similar and was mainly related to tumor pathways. These mainly include the MAPK signaling pathway, cytokine-cytokine receptor interaction, focal adhesion, and other pathways. However, in the early and late subgroups, significantly enriched KEGG domains of DEGs were related to immunity, mainly Chagas disease, T cell receptor signaling pathway, primary immunodeficiency, and natural killer cell-mediated cytotoxicity. Moreover, Cytokine-cytokine receptor and Chagas disease were common pathways in the normal versus tumor subgroups and early versus late subgroups (Figure 1E). The differential genes and functions in the early and late stages of tumors corresponding to the normal group were not obvious; however, the differential genes in the early and late stages mainly played a role in immune function.
Figure 1 Analysis of differentially expressed genes (DEGs) in thyroid carcinoma (THCA) in normal and tumor tissues using The Cancer Genome Atlas (TCGA) data. (A) Volcano plot of DEGs of the two stages. The X axis represents the log2 fold-change, and the Y axis represents the negative log transformation of p values. Up- and downregulated genes are indicated by red and blue colors, respectively. Non-DEGs are indicated by grey color. (B) Distribution of samples in the count group of the database. (C) Intersections of DEGs of the two stages were used to cluster the samples. Normal, early-, and late-stage patients are indicated by blue, red, and green, respectively. (D) Heatmap of 47 genes between normal and tumor tissues. (E) Most significant Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment across normal and tumor subgroup DEGs.
3.3 Neoantigen burden in early- and late-stage samples
At the beginning of the whole exome sequencing data analysis, we briefly summarized the mutation profile at the gene level. Analysis of early and late stages revealed that the missense mutation was the dominant variant (Figures 2A, B). At the gene level, the most frequently mutated genes in the early stages were BRAF (56%), NRAS (9%), HRAS (3%), TG (3%), and TTN (3%) (Figure 2A). The most frequently mutated genes in the late stage were BRAF (62%), NRAS (7%), TG (4%), TTN (4%), HRAS (4%), MUC16 (3%), and USP9X (3%; Figure 2B). The mutation patterns of the mutated genes differed significantly between the early and late stages. We then examined all somatic mutations within the top 5% ranking by frequency; these were selected in early stage and late-stage patients using stacked bar plots. Early and late stages mainly involved missense mutations; BRAF, HRAS, and NRAS were the missense mutations in THCA. We also observed frame mutations in TG and nonsense mutations in TTN in the late-stage group (Figure 2C). BRAF and NRAS are proto-oncogenes, and their mutation spots had only one site (Figures 2D, E). We noted that two conserved domains of BRAF were mutated in the early- and late-stage groups and that these two domains were identical. The Pkinase and PKc-like domains had a missense mutation from V640E to K641E in both stages; however, the P530_Q534del event occurred in early thyroid cancer and the T528_P532del event occurred in late thyroid cancer. BRAF, NRAS, TG, TTN, and HRAS were the most frequently mutated genes in the early and late stages of the tumor; however, other genes were found to be more frequently mutated in the late stage.
Figure 2 Landscape of mutations in THCA. Color-coded matrix of individual mutations in the top 10 most mutated genes in early (A)- and late (B)-stage THCA, indicating the number of recurrent mutations in THCA using TCGA data. (C) Stacked bar plot shows the fraction of variant types in the early and late stages. Mutation spots on BRAF and NRAS proteins in early (D) and late (E) stages. The number of mutations in each spot is shown on the Y axis. Mutation types are colored according to the legend.
3.4 Screening of candidate neoantigens
Given the crucial role of NRGs in THCA, candidate neoantigens were screened and functional analyses were performed in this study. As shown in Figures 3A, B, we first screened the number of neoantigens in the early and late stages and then analyzed the number of associated neoantigen proteins. The results showed that The number of neoantigens in the late stage was significantly higher than that in the early stage. However, the number of neoantigen proteins tended to be balanced between early- and late-stage patients. To further screen the number of candidate neoantigens, Venn diagrams were used to screen the number of candidate neoantigens in normal versus other tumor subgroups (early and late stages; Figures 3C, D). Moreover, 42 neoantigens were identified in the early stage compared to the normal group (Figure 3C), and 55 neoantigens were identified in the late stage compared to the normal group (Figure 3D). Using the PPI network diagram, we identified 16 genes that interacted with the candidate neoantigens in the early stage (Figure 3E) and 28 genes that interacted with the candidate neoantigens in the late stage (Figure 3F). GO function analysis was performed according to the key factors of NRGs in THCA, which were mainly involved in the cAMP-mediated signaling pathway (Supplementary Figure 2). We preliminarily screened a number of candidate neoantigens using TCGA; however, the specific functions still need to be used for diagnostic model construction.
Figure 3 Neoantigen-related genes (NRGs) in THCA. Numbers of neoantigens (A) and associated neoantigen proteins (B) in early- and late-stage patients with THCA. Venn diagram indicating the 42 and 55 candidate neoantigen genes screened in normal versus early stage (C) and normal versus late stage (D) patients with thyroid cancer, respectively. The protein–protein interaction (PPI) networks of candidate neoantigen genes in early (E)- and late (F)-stage patients with THCA were found using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) data.
3.5 Diagnostic model construction
As the candidate neoantigens are tumor-specific and the host genes are overexpressed in tumor tissues, host genes may also be used as diagnostic markers. Based on the preliminary screening of 42 genes in the normal versus early stage group, we used cross-validation to further screen the top 5 NRGs in the early versus normal stage (Figure 4A). The five NRGs identified were G protein-coupled receptor 4, chondroitin sulfate proteoglycan 4 [CSPG4], teneurin transmembrane protein 1 [TENM1], protein S [PROS1], and thymidine kinase 1 [TK1] (Supplementary Table 1). We then used the five NRGs to construct the area under the ROC curve using four machine learning models (mLR, Dtree, RF, and SVM); the lowest predictive value was 0.979 in the normal versus early stages (Figure 4B). Similarly, we preliminarily identified the top 5 NRGs (cadherin 6 [CDH6], semaphorin 6B [SEMA6B], dysferlin [DYSF], xenotropic and polytropic retrovirus receptor 1 [XPR1], and ABR) from 55 candidate neoantigens in late versus normal stages using cross-validation (Figure 4C; Supplementary Table 1). Similarly, the area under the curve (AUC) values of mLR, Dtree, RF, and SVM were 0.996, 0.959, 0.999, and 0.993, respectively, based on the ROC curves of the five NRGs selected from normal and late stages (Figure 4D). We further analyzed the mutation frequency of these 10 NRGs in a normal comparison of tumor stages (early and late). Frequencies of ENM1 and PROS1 mutations in the normal and early groups each accounted for 1%, whereas the other NRGs had no mutation frequency. In addition, no mutation frequency was observed in the 10 NRGs in the late stage of the normal contrast (Figure 4E). The results of TSNAdb (http://biopharm.zju.edu.cn/tsnadb/browse/) were used to predict candidate HLA neoantigen peptide fragments based on the NRGs (shown in Supplementary Table 2). The prediction used NetMHCpan4.0 reveals a strong binding affinity between GPR4, TENM1, and ABR with MHC molecules. Similarly, predictions using NetMHCpan2.8 found that TK1 and XPR1 have strong binding affinities with MHC. This data suggests the potential for these candidate NRGs to develop into new antigens, although further experimentation is required for confirmation.
Figure 4 Feature selection and diagnostic model construction. Early (A)- and late (C)-stage feature filtering using the RFE method. Receiver operating characteristic (ROC) curve of 10-fold cross-prediction performed using four machine learning models in NRGs of normal versus early (B)- and late (D)-stage THCA. (E) Mutation frequencies of candidate neoantigens encoding variants in THCA.
3.6 Immunological correlation analysis of NRGs
CIBERSORT (22) was used to analyze the effect of NRGs on the recruitment of immune cells to the recruitment of immune cells. We first analyzed the recruitment of immune cells to thyroid tissue, as shown in Figure 5A. The main immune cells recruited to the thyroid tissue were T cells CD4 memory resting, macrophages M2, macrophages M0, T cells CD8, and T cells regulatory (Tregs). We further investigated whether the recruitment of immune cells in the early and late groups differed from the recruitment of thyroid tissue described above. Macrophage M1 and naïve B cells, monocytes, and activated dendritic cells were differentially expressed in the early and late groups (p< 0.05), what’s more, T cells CD8, Plasma cells, and CD4 memory activated immune cells were significantly differentially expressed in the early and late groups (p< 0.01; Figure 5B). We further analyzed the correlation between NRGs and the type of immune infiltration, and the results are shown in Figure 5C. The types of immune infiltrates of PROS1, CDH6, XPR1, and ARB neoantigens were similar and were mainly positively correlated with dendritic cell activation and resting dendritic cells, and negatively correlated with Eosinophils and Monocytes. Most NRGs were negatively correlated with activated mast cells, NK cells activated, T cells CD8, B cells memory, and other immune cells. The results showed that the NRGs risk groups strongly correlated with different immune cell types. We also analyzed the immunological profiles of candidate neoantigen genes. Our data suggest that most genes were associated with immunostimulatory factors (Figure 6). ABR, CPR1, CDH6, TK1, and PROS1 were positively correlated, whereas CSPG4 and GRPR4 were negatively correlated with immune checkpoint chemokines, chemokine receptors, MHC molecules, immune stimulators, and inhibitors.
Figure 5 Analysis of candidate neoantigens among different immune subtypes. (A) Boxplot shows the proportion of the 22 types of immune cells in patients with THCA. Box plot (B) and correlation analysis (C) were used to analyze the association between NRGs and immune infiltrating cells in early- and late-stage THCA. *p< 0.05; **p< 0.01; ns, not significant.
Figure 6 Correlation between candidate antigen genes and regulatory factors of relevant immune checkpoints. Chemokines (A). Chemokine Receptors (B). MHC Molecules (C). Immunostimulators (D). Immunoinhibitors (E).
3.7 Construction of a prognostic model for NRGs
Next, we constructed a risk model to understand the effect of NRG expression on the prognosis of thyroid cancer and to evaluate the potential of BRGs as diagnostic markers. Univariate analysis based on the NRGs revealed that SEMA6B (p< 0.05, hazard ratio [HR] = 1.7) and TENM1 (p< 0.05, HR = 0.53) were independent prognostic factors (Supplementary Figure 3). We successfully constructed a prognostic model using 10 candidate neoantigen genes. The NRGs score (risk = -0.5957*exp2$TENM1-0.9129*exp2$TK1-0.4365*exp2$DYSF+0.4896*exp2$ABR+0.7024*exp2$SEMA6B) was calculated using the formula, and the high- and low-risk groups were obtained. Furthermore, we found significant differences between the two groups by constructing a high-low-risk model; the prognosis of the high-risk group was significantly worse than that of the low-risk group (p = 0.00072; Figure 7A). The ROC curve showed that the AUC values were 0.83, 0.87, and 0.86 at 1, 3, and 5 years, respectively (Figure 7B). Almost all patients who died belonged to the high-risk group (Figure 7C). Heatmap analysis showed that TK1 and TENM1 expression levels were low in the high-risk group and high in the low-risk group. Genes with high expression in the high-risk group and low expression in the low-risk group included SEMA6B, ABR, and DYSF (Figure 7D). We then analyzed the clinical prediction model of NRGs by integrating the clinical characteristics of patients in the TCGA database. The results showed that the independent prognostic factors in the univariate analysis were risk (p< 0.001, HR = 2.7) and age (p< 0.001, HR = 1.2; Figure 7E). Multivariate analysis revealed that the risk (p< 0.001, HR = 2.0) and age (p< 0.001, HR = 1.0) were independent prognostic factors (Figure 7F). Therefore, we believe that NRGs can predict survival outcomes.
Figure 7 Univariate analysis using NRGs and prognostic model construction. (A) Risk model survival curves for NRGs. (B) Time-dependent ROC curves of NRGs to evaluate the performance of the prediction model. NRG expression levels in the high- and low-risk groups were assessed using the (C) risk score distribution and (D) heatmap. (E) Univariate analysis based on risk, age, gender, and TNM stage. (F) Univariate analysis based on age.
3.8 Analysis and comparison of NRG expression levels and mutation frequencies in various tumors
Expression patterns of NRGs in different tumor tissues were determined using the TCGA database. The candidate NRGs were all highly expressed in thyroid cancer tissues, and candidate neoantigen genes, such as CDH6, ODZ1, and PROS1, were significantly differentially expressed in thyroid cancers. TK1 was expressed at low levels in KICH but highly expressed in 13 other tumor types (Figure 8A). To analyze the differential expression of NRGs, we examined the differential survival of NRGs in multiple tumors. As shown in Figure 8B, candidate neoantigen SEMA6B was associated with a worse prognosis, and the difference was obvious in the high-expression group of THCA. However, ODZ1 expression was associated with a poor prognosis in the low THCA expression group. Candidate neoantigens, TK1 and XPR1, were associated with poor prognosis in the high expression group for most tumors. NRGs mutate repeatedly in most tumors, especially skin cutaneous melanoma and uterine corpus endometrial carcinoma (Figure 8C). However, the frequency of NRG mutations is relatively low in THCA, possibly related to the favorable prognosis of thyroid tumors.
Figure 8 NRG expression levels and mutation frequencies in multiple tumors. Differential expression (A) and survival difference (B) analyses of NRGs in multiple tumors. (C) Mutation frequencies of NRGs in different tumor samples.
3.9 Analysis of functionally related drugs of candidate neoantigen genes
To analyze NRG-related drug functions, we performed NRG pathway function analysis. NRGs are involved in the activation of epithelial–mesenchymal transition (EMT) but inhibit the DNA damage response. CDH6 plays an inhibitory role in related pathways but participates in the activation of apoptosis, EMT, and cell cycle (Figure 9A). Simultaneously, the candidate-associated neoantigens (CSPG4, XPR1, and CDH6) negatively correlated with the known drugs. However, positive and negative correlations between PROS1 and the drugs occurred simultaneously or alternately (Figure 9B). Through functional analysis of candidate antigen gene drugs, we found that antitumor drugs were involved in antigen presentation. Therefore, the screened NRGs can be used as therapeutic targets and diagnostic biomarkers for thyroid cancer.
Figure 9 NRG functions and related drug analysis. (A) Activities of biological pathways identified using NRGs were analyzed using the GSCA database. Red and azure colors represent the percentage of activation and inhibition, respectively. Non-activity is indicated by grey color. (B) NRG-related drugs were obtained using the GSCA database.
4 Discussion
Thyroid cancer is the most common endocrine tumor and widespread cancer in the USA (2, 28). Most thyroid cancers arise from thyroid follicular cells (90%) and are well-differentiated (29, 30). Most tumors are categorized on histological grounds as papillary thyroid cancers or anaplastic thyroid cancers, the latter being associated with a worse prognosis (2, 31). The standard therapeutic approach for all thyroid cancers includes surgery, with radioactive iodine offered to patients with follicular cell-derived thyroid cancers (1, 32, 33). Several preclinical studies have suggested the potential of immunotherapy for the treatment of thyroid cancer (30, 34). Although this general approach is less developed in terms of clinical trial data, several ongoing immunotherapy trials exhibit great clinical potential. However, the efficacy of immunotherapy for specific THCA stages remains unclear (30, 35, 36).
To predict the potential of candidate neoantigen genes for THCA screening, we comprehensively analyzed the data from the TCGA database and combined them with TSNAdb to predict the HLA candidate neoantigen peptide fragments in this study. Our findings revealed that NRGs had predictive correlations with THCA in terms of gene mutation frequency, DEGs, and type of immune cell infiltration. Expression patterns of downregulated genes were significantly different between the early and late stages of THCA. DEG analysis further revealed that the tumor samples were significantly different from the normal samples (Figure 1). Analysis of the top ten DEGs in candidate patients with early- and late-stage cancer revealed that the mutated genes were mainly concentrated in missense mutations (Figure 2). Genomic instability can be used as a marker of malignant tumors and plays an important role in the development and progression of tumors (37, 38). It not only promotes the evolution of tumors but also assumes a high neoantigen load by tumor cells, which is recognized and localized by the immune system (39, 40). Therefore, we preliminarily screened DEGs and analyzed their mutation types to prepare candidate neoantigens in this study. The ROC curve, consisting of four machine learning models (mLR, Dtree, RF, and SVM), was established via the preliminary screening of candidate neoantigen genes in the early and late groups (AUC of the normal and early groups was > 0.979 and that of the normal and late groups was > 0.959). We believe that the screened NRGs can be used as novel diagnostic biomarkers for thyroid cancer (Figure 3). Because neoantigen loading can be used as a potential biomarker for tumor immunotherapy, the accurate and effective prediction of neoantigens as therapeutic targets may aid in the development of personalized cancer vaccines.
Vaccines are commonly used to prevent infections. Interestingly, vaccines, especially neoantigen-targeted vaccines, also show promise for personalized immunotherapy of cancers, such as hepatocellular carcinoma, melanoma, and epithelial ovarian cancer, in preclinical and clinical studies (41–43). Personalized vaccines are designed to trigger tumor-specific T cell responses against neoantigens to expand the endogenous repertoire of tumor-specific T cells and prevent “off-target” damage to non-tumor tissues (44, 45). Here, we analyzed immune cell infiltration in thyroid cancer and found that the majority of immune-infiltrating cells were T cells (Figure 5A). Furthermore, we compared the immune infiltration in patients with early versus late thyroid cancer and found that the expression levels of T cells CD8, plasma cells, and CD4 + memory-activated immune cells were significantly different (p< 0.01; Figure 5B). Further analysis of the infiltration of 10 candidate NRGs revealed that dendritic and T cells had different degrees of positive and negative correlations with multiple NRGs (Figure 5C). Studies have shown that dendritic cells are primarily involved in acquiring, processing, and presenting tumor-associated antigens on MHC molecules in the tumor microenvironment (TME), and provide costimulators and soluble factors to shape T cell responses (46). Peng et al. reported that dendritic cells participate in antigen presentation and mediate the activation and reactivation of tumor-specific T cells through the ability of endogenous T cell compartments to recognize peptide epitopes that display MHC on the surface of malignant tumor cells (10, 47). We identified NRGs whose neoantigen epitopes result from tumor-specific DNA alterations leading to the formation of new protein sequences that can be used for post-treatment immune memory via neoantigen-specific T-cell responses and prevention of cancer recurrence.
NRGs present a new technique for thyroid cancer immunotherapy. Here, we found that NRGs were expressed in patients with multiple cancers and were significantly correlated (Figure 8A). The expression trends and mutation frequencies of the THCA candidate neoantigens were similar to those in KIRC. Personalized mRNA vaccines based on neoantigens can be used in expression strategies for KIRC solid tumors (48). Candidate neoantigens TK1 and XPR1 were expressed in multiple tumors, and the prognosis was poor in the high expression group of most tumors (Figure 8B). Xie et al. reported that TK1 deactivation significantly inhibits the growth of prostate tumors and is closely related to cell cycle regulation (49). Yoko et al. demonstrated that XPR1-dependent phosphate effervescence leads to the toxic accumulation of intracellular phosphate, inducing growth arrest and apoptosis in ovarian clear cell carcinoma cells (50, 51). To determine the correlation between NRGs and targeted drugs, we used the open GSCA database to identify gene-related targeted drugs for further analysis (Figure 9B). CSPG4, a high-molecular-weight melanoma-associated antigen, is negatively associated with multiple targeted drugs, such as cell-surface proteoglycans. Elevated CSPG4 expression is observed in several aggressive tumors, including ovarian cancer (52), osteosarcoma (53), and triple-negative breast cancer (54). In addition, Egan et al. demonstrated that CSPG4 is a potential immunotherapeutic target for ATCs (55). Bortezomib, as a chemotherapy agent, can trigger immunogenic cell death, thereby promoting anti-tumor immunity (56), and our candidate antigen, CSPG4, was strongly associated with Bortezomib. Trametinib, selumetinib, and PD-0325901 are strongly associated with targeting candidate antigens CSPG4 and PROS1, as the more commonly used MEK inhibitors. Zheng et al. found that the MEK inhibitors combined with radiotherapy can enhance antitumor immunity, and offer a new treatment strategy for KRAS mutations in the tumor (57). Therefore, drugs targeting the candidate antigen, CSPG4, can potentially be used to treat thyroid cancers. In conclusion, a comparison of the expression levels of NRGs among various tumors and their survival analysis revealed that NRGs can predict other tumors, indicating their potential for immunotherapy. Therefore, NRGs can be used as molecular markers for patients with thyroid cancer and converted into antigen-related mRNA vaccines for immunotherapy.
In this study, we screened candidate neoantigens as diagnostic features and therapeutic targets and identified specific candidate neoantigen genes in early- and late-stage thyroid cancer by combining transcriptome and somatic cell mutations. Using machine learning models, we found that the identified candidate genes successfully predicted early- and late-stage thyroid cancer. Simultaneously, we constructed a prognostic model and conducted pan-cancer and targeted drug analyses. Our results revealed that some candidate genes potentially act as candidate neoantigens. Therefore, the candidate tumor-specific neoantigens identified in this study can be used for the personalized treatment of patients with various thyroid tumors.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
Conceptualization and Methodology: JW and MJ; Writing – Original Draft Preparation and Formal Analysis: MJ and JL; Writing – Review and Editing: JL and ZL; Investigation: YQ and QL; Supervision and Project Administration: XL. All the authors have reviewed the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This research was supported by the National Natural Science Foundation of China Youth Science Fund Project (No. 81902881), Academic Leader of Young and Middle-aged Health in Henan Province (Thyroid Surgery; No. HNSWJW-2020004), and the Leading Talents Plan in Henan Province (No. ZYYCYU202012116).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1187160/full#supplementary-material
References
1. Carling T, Udelsman R. Thyroid cancer. Annu Rev Med (2014) 65:125–37. doi: 10.1146/annurev-med-061512-105739
2. Cabanillas ME, McFadden DG, Durante C. Thyroid cancer. Lancet (2016) 388(10061):2783–95. doi: 10.1016/S0140-6736(16)30172-6
3. Kreissl MC, Janssen M, Nagarajah J. Current treatment strategies in metastasized differentiated thyroid cancer. J Nucl Med (2019) 60(1):9–15. doi: 10.2967/jnumed.117.190819
4. Singh ON, Iniguez-Ariza NM, Castro MR. Thyroid nodules: diagnostic evaluation based on thyroid cancer risk assessment. BMJ (2020) 368:l6670. doi: 10.1136/bmj.l6670
5. Lincango-Naranjo E, Solis-Pazmino P, El KO, Salazar-Vega J, Garcia C, Ledesma T, et al. Triggers of thyroid cancer diagnosis: a systematic review and meta-analysis. Endocrine (2021) 72(3):644–59. doi: 10.1007/s12020-020-02588-8
6. Banerjee M, Muenz DG, Worden FP, Wong SL, Haymart MR. Conditional survival in patients with thyroid cancer. Thyroid (2014) 24(12):1784–9. doi: 10.1089/thy.2014.0264
7. Fligor SC, Lopez B, Uppal N, Lubitz CC, James BC. Time to surgery and thyroid cancer survival in the United States. Ann Surg Oncol (2021) 28(7):3556–65. doi: 10.1245/s10434-021-09797-z
8. Matsushita H, Vesely MD, Koboldt DC, Rickert CG, Uppaluri R, Magrini VJ, et al. Cancer exome analysis reveals a T-cell-dependent mechanism of cancer immunoediting. Nature (2012) 482(7385):400–4. doi: 10.1038/nature10755
9. Gubin MM, Artyomov MN, Mardis ER, Schreiber RD. Tumor neoantigens: building a framework for personalized cancer immunotherapy. J Clin Invest (2015) 125(9):3413–21. doi: 10.1172/JCI80008
10. Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science (2015) 348(6230):69–74. doi: 10.1126/science.aaa4971
11. Tran E, Ahmadzadeh M, Lu YC, Gros A, Turcotte S, Robbins PF, et al. Immunogenicity of somatic mutations in human gastrointestinal cancers. Science (2015) 350(6266):1387–90. doi: 10.1126/science.aad1253
12. Hinrichs CS, Restifo NP. Reassessing target antigens for adoptive T-cell therapy. Nat Biotechnol (2013) 31(11):999–1008. doi: 10.1038/nbt.2725
13. Singh G, Agarwal AK, Prosek J, Jayadev M, Singh A. Can killers be saviors? Lupus (2017) 26(9):903–8. doi: 10.1177/0961203316688783
14. Adalsteinsson VA, Ha G, Freeman SS, Choudhury AD, Stover DG, Parsons HA, et al. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors. Nat Commun (2017) 8(1):1324. doi: 10.1038/s41467-017-00965-y
15. Hansen MC, Haferlach T, Nyvold CG. A decade with whole exome sequencing in haematology. Br J Haematol (2020) 188(3):367–82. doi: 10.1111/bjh.16249
16. Karasaki T, Nagayama K, Kuwano H, Nitadori JI, Sato M, Anraku M, et al. Prediction and prioritization of neoantigens: integration of RNA sequencing data with whole-exome sequencing. Cancer Sci (2017) 108(2):170–7. doi: 10.1111/cas.13131
17. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007
18. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118
19. Szklarczyk D, Kirsch R, Koutrouli M, Nastou K, Mehryary F, Hachilif R, et al. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res (2023) 51(D1):D638–46. doi: 10.1093/nar/gkac1000
20. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303
21. Jayawardana K, Schramm SJ, Tembe V, Mueller S, Thompson JF, Scolyer RA, et al. Identification, review, and systematic cross-validation of microRNA prognostic signatures in metastatic melanoma. J Invest Dermatol (2016) 136(1):245–54. doi: 10.1038/JID.2015.355
22. Ye L, Zhang T, Kang Z, Guo G, Sun Y, Lin K, et al. Tumor-infiltrating immune cells act as a marker for prognosis in colorectal cancer. Front Immunol (2019) 10:2368. doi: 10.3389/fimmu.2019.02368
23. Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat Med (2018) 24(10):1545–9. doi: 10.1038/s41591-018-0157-9
24. Chiba Y. Kaplan-Meier curves for survivor causal effects with time-to-event outcomes. Clin Trials (2013) 10(4):515–21. doi: 10.1177/1740774513483601
25. Zhou M, Zhao H, Wang Z, Cheng L, Yang L, Shi H, et al. Identification and validation of potential prognostic lncRNA biomarkers for predicting survival in patients with multiple myeloma. J Exp Clin Cancer Res (2015) 34(1):102. doi: 10.1186/s13046-015-0219-5
26. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med (2013) 32(30):5381–97. doi: 10.1002/sim.5958
27. Liu CJ, Hu FF, Xia MX, Han L, Zhang Q, Guo AY. GSCALite: a web server for gene set cancer analysis. Bioinformatics (2018) 34(21):3771–2. doi: 10.1093/bioinformatics/bty411
28. Seib CD, Sosa JA. Evolving understanding of the epidemiology of thyroid cancer. Endocrinol Metab Clin North Am (2019) 48(1):23–35. doi: 10.1016/j.ecl.2018.10.002
29. Roman BR, Morris LG, Davies L. The thyroid cancer epidemic, 2017 perspective. Curr Opin Endocrinol Diabetes Obes (2017) 24(5):332–6. doi: 10.1097/MED.0000000000000359
30. Laha D, Nilubol N, Boufraqech M. New therapies for advanced thyroid cancer. Front Endocrinol (Lausanne) (2020) 11:82. doi: 10.3389/fendo.2020.00082
31. Chmielik E, Rusinek D, Oczko-Wojciechowska M, Jarzab M, Krajewska J, Czarniecka A, et al. Heterogeneity of thyroid cancer. Pathobiology (2018) 85(1-2):117–29. doi: 10.1159/000486422
32. Hadoux J, Pacini F, Tuttle RM, Schlumberger M. Management of advanced medullary thyroid cancer. Lancet Diabetes Endocrinol (2016) 4(1):64–71. doi: 10.1016/S2213-8587(15)00337-X
33. Tuttle RM. Controversial issues in thyroid cancer management. J Nucl Med (2018) 59(8):1187–94. doi: 10.2967/jnumed.117.192559
34. Ferrari SM, Fallahi P, Galdiero MR, Ruffilli I, Elia G, Ragusa F, et al. Immune and inflammatory cells in thyroid cancer microenvironment. Int J Mol Sci (2019) 20(18), 4413. doi: 10.3390/ijms20184413
35. Chang LS, Barroso-Sousa R, Tolaney SM, Hodi FS, Kaiser UB, Min L. Endocrine toxicity of cancer immunotherapy targeting immune checkpoints. Endocr Rev (2019) 40(1):17–65. doi: 10.1210/er.2018-00006
36. Sun J, Shi R, Zhang X, Fang D, Rauch J, Lu S, et al. Characterization of immune landscape in papillary thyroid cancer reveals distinct tumor immunogenicity and implications for immunotherapy. Oncoimmunology (2021) 10(1):e1964189. doi: 10.1080/2162402X.2021.1964189
37. Andor N, Maley CC, Ji HP. Genomic instability in cancer: teetering on the limit of tolerance. Cancer Res (2017) 77(9):2179–85. doi: 10.1158/0008-5472.CAN-16-1553
38. Shi Z, Shen J, Qiu J, Zhao Q, Hua K, Wang H. CXCL10 potentiates immune checkpoint blockade therapy in homologous recombination-deficient tumors. Theranostics (2021) 11(15):7175–87. doi: 10.7150/thno.59056
39. Germano G, Lamba S, Rospo G, Barault L, Magri A, Maione F, et al. Inactivation of DNA repair triggers neoantigen generation and impairs tumour growth. Nature (2017) 552(7683):116–20. doi: 10.1038/nature24673
40. Mardis ER. Neoantigens and genome instability: impact on immunogenomic phenotypes and immunotherapy response. Genome Med (2019) 11(1):71. doi: 10.1186/s13073-019-0684-0
41. Ott PA, Hu Z, Keskin DB, Shukla SA, Sun J, Bozym DJ, et al. An immunogenic personal neoantigen vaccine for patients with melanoma. Nature (2017) 547(7662):217–21. doi: 10.1038/nature22991
42. Cai Z, Su X, Qiu L, Li Z, Li X, Dong X, et al. Personalized neoantigen vaccine prevents postoperative recurrence in hepatocellular carcinoma patients with vascular invasion. Mol Cancer (2021) 20(1):164. doi: 10.1186/s12943-021-01467-8
43. Sarivalasis A, Morotti M, Mulvey A, Imbimbo M, Coukos G. Cell therapies in ovarian cancer. Ther Adv Med Oncol (2021) 13:17546865. doi: 10.1177/17588359211008399
44. Blass E, Ott PA. Advances in the development of personalized neoantigen-based therapeutic cancer vaccines. Nat Rev Clin Oncol (2021) 18(4):215–29. doi: 10.1038/s41571-020-00460-2 www.ajcr.us/ISSN:2156-6976/ajcr0142628
45. Deng Z, Zhan P, Yang K, Liu L, Liu J, Gao W. Identification of personalized neoantigen-based vaccines and immune subtype characteristic analysis of glioblastoma based on abnormal alternative splicing. Am J Cancer Res (2022) 12(8):3581–600.
46. Wculek SK, Cueto FJ, Mujal AM, Melero I, Krummel MF, Sancho D. Dendritic cells in cancer immunology and immunotherapy. Nat Rev Immunol (2020) 20(1):7–24. doi: 10.1038/s41577-019-0210-z
47. Peng Q, Qiu X, Zhang Z, Zhang S, Zhang Y, Liang Y, et al. PD-L1 on dendritic cells attenuates T cell activation and regulates response to immune checkpoint blockade. Nat Commun (2020) 11(1):4835. doi: 10.1038/s41467-020-18570-x
48. Xu H, Zheng X, Zhang S, Yi X, Zhang T, Wei Q, et al. Tumor antigens and immune subtypes guided mRNA vaccine development for kidney renal clear cell carcinoma. Mol Cancer (2021) 20(1):159. doi: 10.1186/s12943-021-01465-w
49. Zhang C, Ma S, Hao X, Wang Z, Sun Z. Methylation status of TK1 correlated with immune infiltrates in prostate cancer. Front Genet (2022) 13:899384. doi: 10.3389/fgene.2022.899384
50. Akasu-Nagayoshi Y, Hayashi T, Kawabata A, Shimizu N, Yamada A, Yokota N, et al. PHOSPHATE exporter XPR1/SLC53A1 is required for the tumorigenicity of epithelial ovarian cancer. Cancer Sci (2022) 113(6):2034–43. doi: 10.1111/cas.15358
51. Bondeson DP, Paolella BR, Asfaw A, Rothberg MV, Skipper TA, Langan C, et al. Phosphate dysregulation via the XPR1-KIDINS220 protein complex is a therapeutic vulnerability in ovarian cancer. Nat Cancer (2022) 3(6):681–95. doi: 10.1038/s43018-022-00360-7
52. Harrer DC, Schenkel C, Berking C, Herr W, Abken H, Dorrie J, et al. Decitabine-mediated upregulation of CSPG4 in ovarian carcinoma cells enables targeting by CSPG4-specific CAR-T cells. Cancers (Basel) (2022) 14(20):5033. doi: 10.3390/cancers14205033
53. Riccardo F, Tarone L, Iussich S, Giacobino D, Arigoni M, Sammartano F, et al. Identification of CSPG4 as a promising target for translational combinatorial approaches in osteosarcoma. Ther Adv Med Oncol (2019) 11:432496765. doi: 10.1177/1758835919855491
54. Hu ZY, Zheng C, Yang J, Ding S, Tian C, Xie N, et al. Co-expression and combined prognostic value of CSPG4 and PDL1 in TP53-aberrant triple-negative breast cancer. Front Oncol (2022) 12:804466. doi: 10.3389/fonc.2022.804466
55. Egan CE, Stefanova D, Ahmed A, Raja VJ, Thiesmeyer JW, Chen KJ, et al. CSPG4 Is a potential therapeutic target in anaplastic thyroid cancer. Thyroid (2021) 31(10):1481–93. doi: 10.1089/thy.2021.0067
56. Galluzzi L, Buque A, Kepp O, Zitvogel L, Kroemer G. Immunogenic cell death in cancer and infectious disease. Nat Rev Immunol (2017) 17(2):97–111. doi: 10.1038/nri.2016.107
Keywords: thyroid carcinoma, neoantigens, machine learning, immune, prognosis
Citation: Jia M, Liang J, Li Z, Qin Y, Li Q, Wang J and Lu X (2023) Screening tumor stage-specific candidate neoantigens in thyroid adenocarcinoma using integrated exome and transcriptome sequencing. Front. Immunol. 14:1187160. doi: 10.3389/fimmu.2023.1187160
Received: 15 March 2023; Accepted: 11 September 2023;
Published: 03 October 2023.
Edited by:
Reid Rubsamen, Case Western Reserve University, United StatesReviewed by:
Paul E. Harris, Columbia University, United StatesScott Burkholz, Flow Pharma, United States
Copyright © 2023 Jia, Liang, Li, Qin, Li, Wang and Lu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiubo Lu, doctorluxiubo@126.com; Jianwei Wang, wjw725@126.com