Next-Generation Sequencing Identifies a Highly Accurate miRNA Panel That Distinguishes Well-Differentiated Thyroid Cancer from Benign Thyroid Nodules

Patients and biopsies

The study protocol was reviewed and approved by the Independent Ethics Committee, Hadassah-Hebrew University Medical Center (Jerusalem, Israel, Helsinki committee, protocol #HMO-0375). Adult patients who were referred to elective thyroid surgery between January and December 2015 were offered to participate in the study. The investigators obtained written informed consent from patients with a preoperative diagnosis of benign goiter, PTC, or indeterminate cytology. Once the thyroid gland was surgically removed, an ex vivo FNAB was performed as described previously (9). In brief, in patients with thyroid cancer, the nodule was palpated and 2 to 3 passes were performed using a 19-guage needle. In patients with benign thyroid disease, the biopsy was performed from normal appearing thyroid tissue. All preoperative diagnoses were confirmed by standard histopathology. The “validation group” consisted of patients with preoperative indeterminate cytology report (Bethesda III or IV). In these patients, the ultrasound report was carefully reviewed to confirm the exact location of the thyroid nodule that was biopsied before surgery. Ex vivo biopsies were performed from the very same nodule as described above. Samples were coded without identification of the group classification (benign, malignant, validation). Data were collected on each patient including age, gender, indications for surgery, detailed preoperative ultrasound report, preoperative cytology diagnosis, and detailed final histopathology diagnosis. All data were recorded on a Case Report Form designed for this study.

RNA extraction and small RNA cDNA library preparation and sequencing

All cell samples were collected from the aspiration biopsy needle by air injection and irrigation into a collecting tube containing 1 mL of RNAlater (ThermoFisher Scientific). Tubes were stored at 4°C until further processing. Total RNA was extracted with TRIzol, following Polytron homogenization using disposable tips. Concentration was measured with Qubit 2.0 Fluorometer and RNA BR assay reagents (Thermo Fisher Scientific), and RNA was stored at −80°C.

Total RNA was subjected to in-house multiplexed small RNA cDNA library preparation, which entailed ligation of barcoded 3′ adapters to 24 different samples, pooling of samples, ligation of a 5′ adapter, reverse transcription, and PCR, as described previously (10), with modifications allowing multiplexing of several 24-sample libraries on a single HiSeq lane, namely, 48 to 120 small RNA samples per lane (11), using TruSeq-like PCR primers. Libraries were sequenced on an Illumina HiSeq 2500 sequencer, and the information obtained was analyzed by an automated computer pipeline to decode and annotate small RNA reads (12).

Statistical analysis

Statistical procedures on count data were based on the “DESeq2” R/Bioconductor package for analysis of differential expression in RNA sequencing experiments (13), as described previously (14), and complemented by calculations with “edgeR” (15) and “samr” (16). Differential expression analyses were performed upon a matrix of read counts with samples as columns and miRNA as rows, after filtering out low-coverage miRNA (retaining miRNA with at least 3 reads in at least 1/3 of study samples). To overcome the difficulty of multiple testing, the DESeq2 package uses an adjusted P value for which we determined that less than 0.05 would be regarded as significant. The samr package utilizes a Q-value to address multiple testing, which we also determined to be significant if <0.05. Prior to any analysis based on miRNA read counts, we applied the “sva” package (17) to remove batch biases from the data. In addition, prior to analyses other than differential expression, miRNA read counts were transformed using DESeq2′s variance stabilizing transformation, to generate normally distributed data. To describe the extent of technical variability (namely differences between miRNA profiles generated from repeated aspiration biopsies from the same nodule or normal tissue), we calculated correlation distances between samples based on miRNA abundance profiles, and compared the distribution of distances derived from all possible pairwise tests to the distribution of distances from technical repeat pairwise tests, using empirical cumulative distribution function plots and a Mann–Whitney U test. When testing for differential expression, we first aggregated miRNA counts from all technical repetitions of each specimen (separately for tumor and adjacent normal tissue), such that each biological sample was represented once, by the sum of counts of its technical repetitions. In the differential expression models, patient was included as a blocking factor, being that often malignant and benign tissues were obtained from the same specimen. The ability of differentially expressed miRNA to discriminate malignant from benign/normal tissue was evaluated in the discovery set of samples by calculating the c-statistic, equivalent to the area under the ROC curve. To maximize the discrimination ability of miRNA, we devised a score for each sample (“miRNA score”), assembled by summing the standardized levels (z-values) of all significantly upregulated miRNA, and the negative of the z-values of all significantly downregulated miRNA. In the validation set of samples, z-values were calculated by relating to the normal mean and SD values from the discovery set. ROC curve was used to test the discriminatory characteristics of the miRNA score in the validation set. Sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and accuracy were calculated at an optimal cutoff point. Plots were generated with “ggplot2.”

Availability of data

miRNA deep sequencing count data for all study samples, along with partial metadata, are available as a spreadsheet file online. Demultiplexed sequencing read files were deposited at NCBI's Gene Expression Omnibus (GEO record GSE116196).