This article is mentioned in:

Introduction

The incidence of thyroid cancer has been increased over the past few years (1–3). Thyroid nodules are one of most prevalent thyroid diseases, detectable by cervical echography in between 50 and 67% of healthy individuals (4). A confirmation study is required to verify the diagnosis as only 5% of these thyroid nodules are malignant (4).

Usually, diagnosing thyroid nodules as benign or malignant is performed by cytological evaluation (5). For this purpose, fine-needle aspiration (FNA) is the most commonly used sample extraction technique, since it is rapid, inexpensive and simple, as shown by Knezević-Usaj et al (6). Subsequent cytological analysis following FNA provides four possible results: non-diagnostic, positive to malignancy or suspicious, indeterminate and benign cytology (7). Indeterminate FNA cytological results are obtained in between 15 to 30% of cases (8–11). Moreover, only between 5 and 15% of cases are malignant, particularly in those with indeterminate cytological results (11). Consequently, patients with FNA indeterminate cytological results undergo diagnostic surgery, even though this has been proven to be unecessary in >50% of cases where patients are later found to have benign disease (5,12).

Microarray-based gene expression profiling studies of thyroid nodules have proposed molecular markers (11,13,14). However, it has been demonstrated in other types of cancer that some biomarkers identified in one cohort may fail to reproduce similar results with a high degree of accuracy in other cohorts (15,16). For example, the accuracy of the Afirma® genomic test for FNA thyroid samples, which is based on the expression of >170 genes and one of the most extensively studied, has been confirmed by some studies (17–20); however, it has been seriously questioned by more recent investigations in terms of its sensitivity, cost-effectiveness, or its ability to complement tests with a high specificity such as the BRAF mutation test (21–24).

Given that those patients with thyroid nodules of indeterminate FNA cytology may undergo unnecessary surgical intervention, and that previously proposed molecular biomarkers cannot be used in these cases, or that the accuracy of a currently available test has been questioned, a molecular, robust, biomarker for FNA, that is simple, and cost-effective, is still required for clinical investigations and practice. In contrast to other authors, in this study, we propose a molecular biomarker designed specifically from FNA indeterminate thyroid samples identified by a bioinformatics approach, which has been validated in six datasets, including four from other authors. We demonstrate that the accuracy of the proposed biomarker is superior to other previously proposed biomarkers for thyroid tumors. The proposed biomarker is composed of 15 genes and has the potential to be easily implemented into clinical practice using common and cost-effective real-time-polymerase chain reaction (RT-PCR) assays.

Data collection methods

Datasets and processing

We used six gene expression micro-array datasets from five different authors (Table I), which we obtained from large microarray repositories. The main inclusion criteria were that the number of samples was >40 and that the study contained histopathological diagnoses. To compare the results from the different datasets and microarray platforms, we transformed the gene expression data to a uniform distribution between 0 and 1, where 0 represents the lowest and 1 the highest expression. Multiple probes assigned to the same gene were averaged if they were correlated using a Pearson coefficient of ≥0.7. The probe with the highest expression was used if duplicate symbols remained. To facilitate future biomarker measurements in clinical practice using RT-PCR, which may use an internal control for normalization (25), we transformed the original Alexander dataset of 173 genes (11), which represent the previously identified Afirma® test, to a dataset of all combinations of gene-by-gene expression differences. This generated a dataset of 2,850 gene expression differences. In preliminary experiments, we observed that differences allowed better prediction than the raw expression measure (data not shown), which is consistent with other observations, where pairs of genes are more accurate predictors than separate genes, as shown by Grate (26).

Table I Characteristics of datasets used.

Table I

Characteristics of datasets used.

Authors/(Refs.) dataset/(use)	ID/Platform	Sample characteristics	No. of benign/malignant samples	Diagnosis
Alexander et al (11) indeterminate (training set)	GSE34289	265 Indeterminate FNA:	180/85	FNA cytology
	Affymetrix	180 Benign after surgery (B)
	Afirma-T (custom) 173 probes	85 Malignant after surgery
Giordano et al (13) (test set)	GSE27155	89 Adenomas/carcinomas:	17/72	Surgical pathology
	Affymetrix	10 Follicular adenomas (B)
	HG_U133A 22,283 probes	7 Oncocytic adenomas (B)
		13 Follicular carcinomas
		8 Oncocytic carcinomas
		51 Papillary carcinomas
Borup et al (14) (test set)	E-MEXP-2442	69 Adenomas/carcinomas:	45/24	Surgical pathology
	Affymetrix	22 Follicular adenomas (B)
	HG U133 Plus 2.0 54,613 probes	12 Microfollicular adenomas (B)
		9 Nodular goiters (B)
		18 Follicular carcinomas
		4 Anaplastic carcinomas
		2 Papillary carcinomas
		2 Normal (B)
Alexander et al (11) determinate (test set)	GSE34289	99 Determinate FNA:	44/55	FNA
	Affymetrix	44 Benign after surgery (B)		cytology
	Afirma-T (custom) 173 probes	55 Malignant after surgery
TCGA (test set) (https://tcga-data.nci.nih.gov/tcga/)	Illumina Hi-Seq	547 Thyroid samples:	57/490	Surgical pathology
	RNA-Seq	490 Papillary cancers
	20,500 probes	57 Benign tissues (B)
Tomás et al (51)	GSE33630	105 Thyroid tumor/non-tumor:	45/60	International
Dom et al (52) (test set)	Affymetrix	11 Anaplastic carcinomas		Pathology Panel of the Chernobyl Tissue Bank
	HG_U133 Plus 2.0 54,675 probes	49 Papillary carcinomas 45 Patient-matched non-tumor controls (B)

[i] B, indicates benign samples; FNA, fine-needle aspiration.

Biomarker identification

To the best of our knowledge, the Alexander dataset is the only data providing details of cytologically indeterminate thyroid samples [Alexander et al (11)]; therefore, we used this 'training' dataset as a gold standard in order to identify the biomarker. To discover combinations of gene differences that together yield the optimal classification of malignant and non-malignant samples, we used GALGO, a genetic algorithm for feature selection (27). GALGO is a feature selection approach based on genetic algorithms coupled with a classifier. Briefly, GALGO first generates a population of random combinations of features. Each combination of this population is evaluated using the accuracy of a classifier and the selected features. The genetic algorithm then selects those combinations with higher accuracy, which are subsequently re-combined and changed replacing a gene difference with another. The process is repeated until a predefined number of cycles yields a highly accurate feature combination. Since this process is stochastic, the specific features may change; thus, GALGO typically performs this procedure multiple times. Subsequently, a representative feature combination is selected based on the number of times each feature is present in the highly accurate combinations and a forward selection procedure. In this study, as proposed by GALGO tutorials (http://bioinformatica.mty.itesm.mx/GALGO), we used 300 combinations having five features to select the representative biomarker. For the classes, we used benign and malignant cytology as the sample class. For classification, we used the nearest centroid (NC) method shown in the study by Dabney (28). The NC method is based on centers per gene per class estimated as the mean of the gene expression values from the samples of the same class. The samples were classified as the class with the minimum Euclidean distance. The GALGO tutorial has further details of the genetic algorithm and the NC classifier (27). This procedure was performed using the subset of the Alexander dataset (GSE34289) that corresponded to indeterminate FNA cytology and post-surgery determinate, which is composed of 188 samples, 131 as benign and 57 as malignant. We used only 57 randomly selected benign samples for training to balance the number of benign samples with malignant samples.

Biomarker evaluation

To evaluate the performance of the proposed biomarkers and those proposed by other authors, we used an NC classifier learning the parameters from the gene expression measurements of their corresponding datasets. Given that the Alexander dataset is the only available data providing details of indeterminate thyroid samples (11), we used this dataset as the gold standard to evaluate the performance of biomarkers in the undetermined samples. For the cytologically-determined samples, we used the other five datasets shown in Table I as Test datasets.

Results

Previously proposed thyroid cancer biomarkers are not robust

To predict thyroid tumor malignancy, we compared the accuracy of four previously described biomarkers involving between 3 and 167 genes (10,11,14,29) evaluated in six datasets using an NC classifier. The average accuracy ranged between 73 and 78% (Fig. 1). However, we observed some issues. Firstly, none of these four biomarkers accurately predicted The Cancer Genome Atlas (TCGA) subtypes; the maximum was 55%. Secondly, the accuracies evaluated in the 265 indeterminate FNA samples (Alexander dataset) were poor. Thirdly, two of the biomarkers needed almost 100 genes or more, which would generate technical and economic difficulties in clinical practice.

Figure 1 Evaluation of proposed and previous biomarkers for thyroid cancer in six datasets. The accuracy, shown on the vertical axis, was estimated using a nearest centroid classifier. Biomarkers are shown on the horizontal axis along with the different datasets indicated by the bars in different shades of gray. Average accuracy is shown below each biomarker. The number of genes per biomarker is shown in parenthesis. The number below some bars corresponds to the number of genes found in the corresponding dataset. For the Tomei biomarker, we did not use BRAF exon 15 status.

Identification of a highly accurate and robust 15-gene biomarker

The application of previously proposed biomarkers was associated with several concerns: low accuracy, lack of robustness, need to screen of a high number of genes, as well as poor performance when used to analyze indeterminate samples. These biomarkers were all identified through the study of thyroid samples with a definitive cytological diagnosis. Therefore, we specifically selected thyroid tumors with indeterminate cytology from the Alexander GSE34289 dataset (11). To facilitate measurements in a clinical laboratory using RT-PCR and to improve accuracy, we used all combinations of gene differences (26,30–32) instead of the 173 gene expression profiles in GSE34289. To select a low number of genes, we used a multivariate search to identify optimal combinations (27). This strategy is based on genetic algorithms coupled with an NC classifier. Finally, to validate the proposed biomarker in silico, we used five additional data-sets (Table I).

The average accuracy (83%) of the proposed biomarker was superior to the other biomarkers (Fig. 1). The proposed biomarker was the most accurate in four of the six datasets, including the TCGA dataset and the cytologically indeterminate samples; it was also highly competitive in the remaining two datasets [Borup (14) and Giordano (13)]. The biomarker identified using GALGO was based on 15 gene differences covering 15 genes (CCND1, CLDN16, CPE, LRP1B, MAGI3, MAPK6, MATN2, MPPED2, PFKFB2, PYGL, PTPRE, SEMA3D, SERGEF, SLC4A4 and TIMP1). This signature seems to be preserved across the six datasets, exhibiting clear differences between malignant and benign samples (Fig. 2). Twelve of these gene differences were statistically altered in four, five or six datasets (Table II). We also tested the signature across available strata. We observed a high degree of accuracy which was independent of gender, age, tumor size and ethnicity (Table III).

Figure 2 Profiles of the proposed biomarker in the six datasets. The heatmaps show samples in columns and gene differences in rows. Purple and pink columns represent benign (B) and malignant (M) samples, respectively. The total number of samples per subtype is shown as 'B' and 'M'. Accuracy (Acc) is indicated. Cells within the heatmap are shaded in colors; red for high expression, white for median expression, and blue for low expression. The difference involving the MAGI3 gene that was not found in Giordano was removed.

Table II Differences in centroids between benign and malignant samples across datasets.

Table II

Differences in centroids between benign and malignant samples across datasets.

Gene difference	Giordano		Borup		Alex. Ind		Alex. Det		TCGA		Tomás		Highlyb significant
	Diff	p-valuea	Diff	p-valuea	Diff	p-valuea	Diff	p-valuea	Diff	p-valuea	B	p-valuea
MPPED2 - CPE	0.89	5	4.15	6	0.35	10	3.67	20	−0.05	8	0.29	20	6
LRP1B - CPE	0.45	2	2.91	10	0.39	13	3.53	23	−0.09	31	−0.26	32	5
PYGL - TIMP1	0.29	2	0.98	2	0.19	12	2.61	19	0.17	14	−0.26	NS	3
SLC4A4 - CPE	0.73	3	3.19	4	0.29	8	3.14	18	−0.1	20	−0.21	28	6
PFKFB2 - CLDN16	0.66	6	0.94	NS	0.44	13	4.47	19	0.17	37	0.00	15	5
MATN2 - CPE	0.66	3	1.53	2	0.21	6	2.55	13	0.24	13	0.01	6	5
MAGI3 - CLDN16	–	–	1.18	2	0.51	12	3.44	19	0.08	15	0.33	20	4
PFKFB2 -CCND1	0.25	NS	−0.22	NS	0.17	8	2.04	19	0.18	15	−0.35	10	4
PFKFB2 - CPE	−0.65	NS	−1.23	2	0.13	2	2.13	9	−0.11	23	−0.47	NS	2
SERGEF - PFKFB2	0.04	NS	0.44	NS	−0.14	7	−1.06	7	−0.19	14	0.26	3	4
CPE - MAPK6	−0.26	NS	0.22	NS	−0.4	10	−4.7	24	0.09	26	0.03	4	4
CPE - SEMA3D	−0.34	NS	−1.17	2	−0.38	9	−4.06	17	−0.01	2	0.47	21	3
CCND1 - LRP1B	−0.5	5	−2.66	9	−0.41	15	−3.43	26	−0.1	13	0.14	35	6
PTPRE - PYGL	−0.37	7	−0.4	NS	−0.1	8	−0.96	9	−0.07	14	0.08	6	5
CLDN16 - SEMA3D	−0.79	11	−2.08	3	−0.67	12	−5.8	18	−0.25	41	−0.0	29	6

a p-values are expressed as -log10 (p-values);

b number of datasets with a p-value <0.001. Ind, indeterminate; Det, determinate; Diff, differences in the centroids; NS, not significant.

Table III Accuracy of the biomarker in groups of samples.

Table III

Accuracy of the biomarker in groups of samples.

Authors/(Refs.) Dataset	Groups	Samples	Accuracy
Alexander et al (11) (indeterminate)	Data available	265	0.82
	Male	61	0.82
	Female	204	0.82
	Age ≤47	103	0.79
	Age >47	162	0.83
Alexander et al (11) (determinate)	Data available	102	0.94
	Male	27	0.91
	Female	75	0.94
	Age ≤47	38	0.91
	Age >47	64	0.95
Borup et al (14)	Data available	69	0.86
	Male	23	0.71
	Female	46	0.89
TCGA (https://tcga-data.nci.nih.gov/tcga/)	Data available	555	0.88
	Male	151	0.89
	Female	404	0.87
	Tumor size ≤3 cm	323	0.90
	Tumor size >3 cm	196	0.86
	Asian	54	0.94
	African-American	25	0.86
	Latino	42	0.94
	Caucasian	330	0.90
Tomás et al (51)	Data available	105	0.93
Dom et al (52)	(no strata available)	–	–
Giordano et al (13)	Data available	89	0.74
	(no strata available)	–	–

Genes in biomarker play important roles in cancer

Remarkably, the majority of the 15 genes that compose the biomarker have been previously associated with cancer. CLDN16 has been shown to be elevated in patients with thyroid papillary cancer (33). LRP1B inactivation has been shown to influence the tumor environment, thereby increasing the growth and invasiveness of thyroid cancer cells (34). SLC4A4 is expressed in low levels in papillary thyroid carcinoma (35). TIMP1, an inhibitor of the metalloproteinases in the extracellular matrix (36), has been shown to be highly expressed in thyroid cancer (37,38). CCND1, which is involved in the inactivation of the retinoblastoma (RB) protein, as well as in the G1-S phase transition within the cell cycle, has been shown to be associated with many tumors (39), including thyroid papillary carcinomas and follicular adenomas and carcinomas, as shown by Seybt et al (40). SEMA3D, a semaphorin that guides migrating cells during developmental morphogenesis and in adult tissues (41), has been shown to have anti-tumorigenic properties (42). The expression levels of CLDN16, LRP1B, SLC4A4, TIMP1, CCND1 and SEMA3D across subtypes in the six datasets analyzed in the present study were consistent with the findings of the above-mentioned studies (Fig. 3). PFKFB2, which is involved in the control of glycolysis, has been shown to be highly expressed in patients with papillary thyroid cancer aged >40 years compared with younger patients (43). However, in the present study, PFKFB2 appeared to be more highly expressed in the benign tumors across the datasets. CPE mutations have been related to deficiencies in thyrotropin-releasing hormone (44), suggesting that it plays important roles in the thyroid gland. CPE has been shown to be associated with tumor growth and metastases in pheocromocytomas and others types of cancer (45). In this study, we found a consistently high expression of CPE in malignant tumors; however, CPE expression levels varied in the benign samples. MATN2 and MPPED2 seem to be highly correlated (r=0.7 in benign samples in Alexander dataset) and highly expressed in the benign thyroid gland (Fig. 3). It is well known that the former is involved in the formation of filamentous networks in the extracellular matrix and the latter displays low metallophosphoesterase activity. MPPED2 has been proposed to play an important role in neuroblastoma tumorigenesis (46) and the increased expression of this gene has been shown to be associated with a good prognosis (46), which is consistent with the higher expression observed in the benign thyroid tumors in the present study. By contrast, MATN2 overexpression has been observed in pilocytic astrocytoma (47). MAPK6 is a member of the Ser/Thr protein kinase family that has been found to be associated with tumor invasion in lung cancer (48). Polymorphisms in MAGI3 and PYGL have been associated with various disorders, MAGI3 with hypothyroidism (49) and PYGL with relapse in leukemia (50). This literature review of the involved genes suggests that the majority play or may play an important role in thyroid tumors.

Figure 3 Expression levels across six datasets of the genes included in the proposed biomarker. Cells within the heatmap are shaded in colors; red for high expression, white for median expression, and blue for low expression.

Alterations in correlation coefficients characterize differences in gene expression

Eight genes are included in only one difference and seven in more than one difference (Fig. 2). Notably, the CPE gene was found in seven gene differences indicating an important contribution within the biomarker. We observed a higher correlation of genes combined with CPE in the benign samples compared with those correlations in the malignant tumor samples from the Alexander indeterminate dataset (Fig. 4A). By contrast, a higher correlation between the malignant samples was not observed. A similar analysis of all gene pairs across the datasets confirmed this trend (Fig. 4B).

Figure 4 Differences in correlation coefficients between benign and malignant samples. (A) Spearman correlation coefficients in benign and malignant samples of the genes involved in each difference for the Alexander Indeterminate dataset. Arrows mark differences >0.3. (B) Spearman correlation coefficients (horizontal axis) across datasets (vertical axis) in benign samples (top section of table) and the decrease in spearman correlation coefficients in malignant samples (bottom section of table). Columns sorted as in (A). Downward gray arrows mark a decrease >0.3 in correlation coefficients, while upward arrows show an increase >0.3 in correlation coefficients (negative decrease). The 'conserved' cells show the number of datasets whose difference in correlation coefficients is >0.3. Ind, indeterminate; Det, determinate.

Discussion

Previously proposed biomarkers were not robust across datasets or indeterminate FNA samples when evaluated under similar conditions in silico. This may be due to characteristics of the samples, microarray technology, or the methodology used for biomarker identification. Whereas other studies have focused on determinate samples in order to identify biomarkers (10,11,13,14,29,51,52), we specifically used indeterminate samples as the training set, and therefore, we captured the particular expression signatures in these samples. We then showed that the signatures were also conserved in five studies of determinate tumors, which validates the proposed signature. For this purpose, we used gene expression differences between pairs of genes and a multivariate search methodology. Notably, the proposed biomarker is more compact and more accurate than other previously proposed biomarkers.

The proposed biomarker was found to be robust when evaluated in other databases and across patient characteristics (tumor size, age, gender and ethnicity). These results suggest that differences in expression are independent of the cohort, methodology, genomic technology and particular characteristics of the cohort and thus, it is highly likely to represent true biological alterations. In the Giordano (13) and Borup (14) databases, the biomarker was capable of classifying, with high accuracy, many cellular types of thyroid cancer. In the case of determinate FNA samples in the Alexander database (11), the performance of the biomarker was higher (96%) than that of the indeterminate ones (87%) even though the latter was used to identify the biomarker. This result suggests that indeterminate samples may contain transitional stages between benign and malignant subtypes.

The differences in gene expression allow for the easier measurement in widely used technologies, such as RT-PCR, thereby facilitating implementation in clinical practice. Surprisingly, most of the gene pairs in differences were highly correlated in the benign tumors and poorly correlated in the malignant tumors. This concurs with observations in prostate (53), colon, lung, pancreatic, cervical and gastric cancers (54) where the tumor correlation distribution is different than in normal counterparts, generally sharper around zero. Notably, these results suggest that differences in correlations may be an important characteristic of tumor transformation, which may be exploited for biomarker identification, cancer prognosis and gene targeting. We hypothesized that these differences between malignant and non-maligant samples were important for the multivariate search and the classifier to select the genes involved in the proposed biomarker. This may explain the high number of occurrences of the CPE gene, which showed the largest differences in correlation coefficients (ten in benign and one in malignant samples).

From the 15 genes identified in our biomarker, which is a subset of those from the study by Alexander et al (11), none are similar to those proposed by Prasad et al (HMGA2, MRC2, and SFN) (55) and Tomei et al (KIT, C21orf4, PDK3/Hs.296031, DDI2, CDH1, LSM7 and TC1) (29), and only two (MATN2 and MPPED2) are included in the genes from the study by Borup et al (14). Thus, the proposed signature combined with the use of gene differences and a NC classifier appears to be distinctive.

The contribution of the multivariate search was important since the other methodologies tested, such as PAM-R (56) and support vector machine recursive feature elimination (SVM-RFE) (57), generated lower accuracies (65 and 75%, respectively) or higher numbers of gene differences (85 and 15, respectively). Besides the multivariate search, we believe that the use of the gene-pair difference was an important factor which enabled us to identify the highly accurate marker. We then showed that the gene-pair difference may also associated with the difference in correlation coefficients between genes and tumor subtypes. Nevertheless, this approach is almost prohibited in large datasets since a dataset of 20,000 genes would generate 200 million differences combinations. Thus, the use of the 173 genes (or a stringent gene filter) was also a critical factor.

As the proposed biomarker was only tested in silico, a validation study is warranted to confirm the potential use of this biomarker in clinical practice. Although we aim to explore this line of research in the near future, the availability of the proposed biomarker and the methodology used may encourage other research groups to test the biomarker or to design better ones.

In conclusion, the proposed biomarker is composed of 15 gene differences involving 15 genes. The majority of the genes have been associated with cancer and some specifically with thyroid cancer in the research literature. Our analysis suggests that the proposed biomarker is more accurate and robust than previous thyroid biomarkers in tumors and indeterminate FNA samples. Measuring the biomarker may be made relatively easy by RT-PCR facilitating implementation. Changes in the gene expression correlations between benign and malignant samples may be associated with tumor progression and may explain the presence and robustness of the gene differences that compose the proposed biomarker.

Abbreviations:

Acknowledgments

The present study was supported by Grupo de Investigación con Enfoque Estratégico en Bioinformática of the Instituto Tecnológico y de Estudios Superiores of Monterrey, CONACyT (Posgrado Nacional 002087 and grant scholarship 339770). We thank the Instituto Tecnológico y de Estudios Superiores of Monterrey, Hospital San José de Monterrey, and the Instituto Mexicano del Seguro Social for supporting this study.

References