In this retrospective cohort study, the initial training cohort was derived from the Afirma GSC algorithm training subjects, composed of thyroid nodule patients recruited for the Afirma GEC and subsequent GSC training studies (consecutively collected from 2013 to 2016). These thyroid nodules were mostly ITN, mostly histologically benign, and generally very low-risk thyroid cancer when malignant (4, 6). Given a need to train on samples with outcomes of interest (tumor invasion and locoregional lymph node metastases), a subsequent cohort from an integrative interventional endocrinology and endocrine surgery community practice (Memorial Health, Hollywood, FL, USA) with BIII–VI nodule cytology and malignant final thyroid histopathology was incorporated [consecutive fine-needle aspiration (FNA) dates January 2019 to July 2021]. Together, these cohorts (n = 697) constituted the “training cohort” ( Table 1 ).

After the tumor INV and LNM classifiers were locked, independent cohorts from Memorial Health (n = 63, FNA dates August 2021 to October 2022) and Brigham and Women’s Hospital (n = 196, FNA dates July 2017 to June 2023), all sent consecutively for Afirma testing as part of their routine clinical practice for nodules with BIII–VI cytology, were analyzed as the validation cohort ( Table 1 ). These were consecutive samples with local cytology and histopathology interpretations, and treatment decisions based on the local clinician’s discretion with only commercially available Afirma GSC data.

An evaluation cohort of 17,436 consecutive Afirma-resulted ITN GSC suspicious or Bethesda V/VI samples was derived from the Veracyte CLIA laboratory from routine thyroid nodule molecular testing (2017–2020) (18). The INV and LNM classifiers were applied to assess the proportion of samples ruled out for high risk for invasion label and lymph node metastases by Bethesda cytology category, sex, and mutation type (BRAFV600E, RAS, or no detected expressed alteration).

Patients recruited for the Afirma GEC development and validation study provided written informed consent (4). The samples subsequently used for the Afirma GSC algorithm training were approved by institutional-specific review boards, Chesapeake IRB 15.02.0009 (now Advarra IRB, Columbia, MD, USA), and Copernicus Group Independent Review Board VER3-15-067 (now WCG IRB, Princeton, NJ, USA) (6). Patient data (including cytology and histopathology reports) from Memorial Health were collected under WCG IRB protocol # DHF 005-044, and Brigham and Women’s patient data were collected under WCG IRB protocol # DHF 005-077.

A scoring system was applied to the local pathology thyroid histopathology synoptic report ( Table 1 ). For tumor INV, if pathology reported vascular invasion of ≥4 blood vessels (or described extensive vascular invasion) or there was any extrathyroidal extension, the sample was labeled high risk. Otherwise, the sample was labeled low risk for tumor INV. For LNM, if the pathology reported ≥2-mm central lymph node deposits or ≥40% of the central nodes resected as malignant, or if there was lateral lymph node thyroid cancer involvement, the sample was labeled high risk. Otherwise, the sample was labeled low risk. Cases without lymph node dissection (Nx) were assigned the low-risk label, as routine preoperative imaging to assess lymph node disease is recommended (19), and the American Association of Endocrine Surgeons guidelines for the surgical management of thyroid disease do not recommend a routine or prophylactic neck dissection. A central neck dissection is only recommended in selected cases with imaging or clinical (macroscopic) lymph node disease (20). The cut point for cancer features was targeted to where the 2015 ATA guidelines’ risk of structural disease recurrence diagram (Figure 4 in the guidelines) bridged from low- to intermediate-risk cancers (9). Therefore, minor extrathyroidal extension received the high-risk label to clearly delineate ATA low-risk disease. The ATA guidelines utilize an absolute lymph node number involved of >5 to distinguish ATA intermediate-risk cancers from low-risk cancers. Given a concern for labeling tumors with five of five or four of four positive lymph nodes as low risk for metastases, central metastatic lymph node ratio (MLNR) criteria were used for risk assessment, as Nam et al. and Seok et al. reported that central compartment MLNR of >30% and ≥36%, respectively, were significantly associated with recurrence (21, 22). In both studies, MLNR above the thresholds described was statistically significant for thyroid cancer recurrence, whereas overall lymph node yield was not. Samples could be high risk for one category and low risk for another. High-risk and low-risk descriptors were solely for labeling and are not intended to correlate with ATA thyroid cancer pathology risk or risk of recurrence (9).

RNA-seq data were used to generate gene expression counts. Raw sequencing data (FASTQ file) were aligned to the human reference genome assembly 37 (Genome Reference Consortium) using the STAR RNA-seq aligner. Normalized expression levels were obtained using variance-stabilizing transformation (VST) from the DESeq2 package accounting for sequencing depth and gene-wise variability (23).

For sample quality control, quality metrics were evaluated against prespecified acceptance metrics for total numbers of sequenced and uniquely mapped reads, the overall proportion of exonic reads among mapped reads, the mean per-base coverage, the uniformity of base coverage, and base duplication and mismatch rates. All quality control metrics were generated using RNA-SeQC (24). Only samples that passed the quality criteria were used for downstream analysis. For further details, please see Supplementary Methods in Patel et al. (6)

Histopathology scoring labels (low risk vs. high risk) were used to train machine learning models to classify samples into low- and high-risk categories for invasion and LNM outcomes using both genomic and cytology variables ( Supplementary Table S1 , Supplementary Figure S1 ).

For the invasion classifier, features related to cancer pathway activity, genomic alterations, gene expression, and cytology variables were tested. Pathway/signature scores of 430 gene signature/pathway gene sets from MSigDB were calculated for each sample as described before (25, 26). These pathway scores were used as features for model training. The combinations of several machine learning (ML) models including random forest (RF), penalized generalized linear model (glm), support vector machine (SVM), and several feature engineering methods were evaluated ( Supplementary Figure S2 ). Repeated nested fivefold cross-validation (CV) was used for model training, and parameter optimization was used to reduce overfitting and evaluate model performance. Negative predictive value (NPV), the percentage of patients classified as low risk, and score inter-batch reproducibility were the metrics used for selecting the optimal model. The best-performing model was an RF model that used BRAF status, nine cancer pathways/signatures, and cytology group as features ( Supplementary Figure S2 , Supplementary Table S2 ). For the LNM classifier, the combinations of ML models and feature engineering methods were evaluated using the expression of individual genes, genomic alterations, and cytology groups as features ( Supplementary Figures S1 , S3 ). A similar repeated nested fivefold cross-validation approach was used to find the best model. The best-performing model was a penalized glm that uses BRAF status, cytology group, and the expression of 32 differentially expressed genes ( Supplementary Figure S3 , Supplementary Table S2 ). For classifiers’ reproducibility, 18 samples were used, with each sample run in three different runs with three replicates. These nine replicates/samples were used to calculate inter-batch standard deviation (SD). The inter-batch analytical assessment showed that both classifiers’ scores were reproducible with SD < 5% of the 98% score range [1st percentile–99th percentile]. The final models were retrained on the full training cohort, locked, and then tested in the validation and evaluation cohorts while blinded to the histopathology results.

The classifiers’ cut points were determined using the per-sample median of repeated fivefold CV scores, which resulted in both a high rule-out percent and a high NPV.

There were 379 pathology reports from the Afirma GSC training cohort and 318 pathology reports from Memorial Health that were scored for a total of 697 paired Afirma GSC samples with histopathology outcomes for classifier development ( Table 1 ). There were 152 (21.8%) male and 545 (78.2%) female patients aged 9–86 with a mean age of 51 years [interquartile range (IQR): 38–60]. For tumor invasion, 627 were scored as low risk and 70 as high risk. For LNM, 618 were scored as low risk and 79 as high risk. Among those cases labeled high risk for LNM, where only central nodes were positive, the mean number of nodes resected was eight (median 4, range 1–33 [IQR1–3: 2–14]) with a mean MLNR of 0.68. Fifty-five percent of the training cohort was BIII/IV, and 45% was BV/VI. Thirty-three percent of the samples were BRAFV600E classifier positive. Among all training samples, the prevalence of high-risk scores for invasion and LNM on the surgical pathology report was 10.0% and 11.3%, respectively ( Table 1 ).

The validation cohort included 259 patients, 65 (25.1%) male and 194 (74.9%) female patients aged 16–81 with a mean age of 53 years [IQR: 39–62]. Nodules with BIII/IV cytology and classified as GSC suspicious accounted for 91.5% of the samples, and the rest had BV/VI cytology. Thirty (11.6%) were BRAFV600E classifier positive. Eight (3.0%) were scored as high risk for invasion and 7 (2.7%) as high risk for LNM according to the surgical pathology reports ( Table 1 ).

In the training cohort, the INV classifier had, in fivefold cross-validation, a sensitivity (SN) of 90% [80.5–95.9] and a specificity (SP) of 44.8% [40.8–48.8] and was able to rule out 41.3% of the population for high-risk invasion with a 97.6% NPV ( Figure 1A , Table 2a ). In BIII/IV samples (n = 385), 246 (64%) were ruled out for clinically significant invasion with 98% NPV ( Figure 1B , Table 2b ). The rule-out percentage was similar in male (44.7%) and female patients (40.4%) (Fisher’s exact test p = 0.35) ( Table 3a ). In samples with BV/VI cytology (n = 312), 42 (13.7%) samples were ruled out ( Table 3a ).

In the validation cohort, the INV classifier had an SN of 87.5% [47–100] and was able to rule out 49.4% with 99.2% NPV and a specificity of 50.6% [44.2–56.9] ( Figure 1C , Table 2a ). In BIII/IV samples (n = 237), 127 (53.6%) were ruled out with 99% NPV ( Figure 1D , Table 2b ), and in BV/VI samples (n = 22), one sample (4.5%) was ruled out ( Table 3b ). The one false-negative sample ( Figure 1C ) had a lobectomy with final pathology showing a 1.2-cm infiltrative follicular subtype of papillary thyroid carcinoma (IF-PTC) with extrathyroidal extension into the adjacent strap muscle. Completion thyroidectomy was benign.

There was no significant difference in performance when comparing samples with BIII or BIV cytology ( Tables 2c, d ).

In the training cohort, the LNM classifier had, in fivefold cross-validation, an SN of 94% [85.8–97.9] and an SP of 55% [51–59] and ruled out 49.8% of the population for high-risk LNM with a 98.6% NPV in the training cohort ( Figure 2A , Table 2a ). Of Bethesda III/IV samples (n = 385), 294 (76%) were ruled out for high-risk LNM with 98% NPV ( Figure 2B , Table 2b ), and of those with Bethesda V/VI (n = 312), 51 (16.3%) samples were ruled out ( Table 3a ). The rule-out percentage was similar in male (53%) and female patients (49%) (Fisher’s exact test p = 0.41) ( Table 3a ).

In the validation cohort, 44% of the cases had lymph nodes removed, with 14% of those meeting a threshold of at least six nodes removed, which suggested being adequate as a central neck dissection (27). Fifty percent of the local pathology reports did not make any comment regarding lymph nodes, and these were almost exclusively benign cases or non-invasive follicular thyroid neoplasms with papillary-like nuclear features (NIFTP). The LNM classifier had an SN of 100% [59–100] and ruled out 52.5% with 100% NPV and an SP of 54% [44.6–61.6] ( Figure 2C , Table 2a ). In BIII/IV samples (n = 237), 135 (57%) were ruled out with 100% NPV ( Figure 2D , Table 2b ), and in samples with Bethesda V/VI (n = 22), one sample was ruled out (6.7%) ( Table 3b ).

There was no significant difference in performance when comparing samples with BIII or BIV cytology ( Tables 2c, d ).

The initial surgical intervention was assessed, and all were total thyroidectomy (TT) or lobectomy. In the validation cohort, there were 83 TT (32.5%) ( Table 3b ). Of those with ITN, there were 62 TT (26%). Of all samples with a TT, 16 (19%) had a low-risk invasion classifier alone, 13 (16%) had a low-risk metastasis classifier alone, and 11 (13%) had both low-risk classifiers. Of the 40 (48%) tumors with at least one low-risk classifier, 39 were either histologically benign, NIFTP, or ATA low-risk cancers. The one ATA high-risk cancer was an IF-PTC that had BVI cytology, harboring an NRAS:Q61R variant, was >6 cm in size with extensive vascular invasion, and had a correctly assigned INV classifier (not ruled out for INV) and a correctly assigned low-risk LNM classifier (ruled out for LNM) with N0 on final pathology (0/11 nodes).

An evaluation of 17,346 Afirma GSC-suspicious samples with no clinical outcomes was assessed to compare the tumor classifier scores in an unselected consecutive cohort. These samples were from patients with a median age of 54 years [IQR 40.4–65.9] ( Table 1 ). Sex was 75.5% female, and the Bethesda cytology categories were as follows: 67.5% BIII, 23.2% BIV, 4.6% BV, and 4.7% BVI. Overall, 53.3% had a low-risk invasion INV classifier score, and 44.1% had a low-risk LNM classifier score ( Table 3c ). The percentages of samples ruled out by sex, BRAFV600E (an Afirma GSC classifier with a limit of detection of >5% variant allele frequency considered as positive) (28), RAS, and Xpression Atlas (XA) (17) negative mutation status are shown in Table 3c .