Thyroid tumor samples

10 follicular thyroid carcinoma (FTC) and 10 follicular thyroid adenoma (FA) samples were selected for miRNA profiling by HTS. Hematoxylin and eosin (HE) stained histology sections of all tumor samples were independently rated by two experienced pathologists to guarantee the precise tumor classification. Approval of the ethics commission of the medical faculty of the University of Leipzig and written consent of the respective patients was obtained.

RNA isolation and library preparation

Total RNA was extracted from 20-50 mg snap-frozen thyroid tissue using 1 mL Trizol reagent. Reverse transcription was done with miScript Reverse Transcription kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocols in a 10 μL reaction using 1 μg of total RNA. 500 ng of total RNA prepared from the tumor samples was used in the small RNA protocol with the TruSeq™ Small RNA sample prepkit v2 (Illumina) according to the instructions of the manufacturer. The barcoded libraries were size restricted between 140 and 165 bp, purified and quantified using the Library Quantification Kit - Illumina/Universal (KAPA Biosystems) according to the instructions of the manufacturer.

High-throughput sequencing

A pool of ten libraries was used for cluster generation at a concentration of 10 nM using an Illumina cBot. Sequencing of 50 bp was performed with an Illumina HighScan-SQ sequencer using version 3 chemistry and flowcell according to the instructions of the manufacturer.

High-throughput sequencing data analysis

For naming of the sequences the following nomenclature in concordance with miRBase 18 was used. According to miRBase the entry sequence represents the most frequent isoform of a miRNA. A hairpin miRNA or pre-miRNA sequence is the precursor of up to 2 canonical mature miRNA sequences (denoted as 5p or 3p). Isoforms of such a miRNA sequence vary around the miRBase entry with a certain offset at the 5′ end or 3′ end or a slight difference in nucleotide sequence (miRNA point mutations and SNPs).

High-throughput sequencing data was analyzed according to the three pipelines presented in Figure 1. The aim of the different analysis pipelines is to measure the expression of miRNA with different approaches: the classical approach - alignment of reads, and two new approaches - counting of isoforms which is independent of any database entries and seed analysis which focuses analysis to the bases that are most likely functionally relevant.

At first, we removed adapters sequences from raw reads (50 bases, fastq files) of all 20 samples (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39112) using Cutadapt software [30]. From the remaining sequences only those 15–27 bases long were kept for further analysis because they most likely contain mature miRNA sequences. Sequences with median quality lower than phred = 30 were discarded from the analysis. To monitor the success of the adapter cutting and filtering we used the FastQC program (Bioinformatics Group at the Babraham Institute, UK).

At this point further analysis was split into pathways: (i) mapping based pipeline A intersecting known miRNAs coordinates from miRBase after alignment to the human genome and (ii) pipeline B and C which count identical reads first and then annotate individual reads as miRNA isoform when they contain a miRBase entry sequence or map to hairpin miRNA sequences. Moreover sequences that did not succeed in this annotation step were annotated as uncharacterized small RNA sequences if they align to the human genome. Pipelines B and C aim for de novo analysis of thyroid tumor miRNA and small RNA profiles and therefore include all sequences independent of a homology to known miRNA molecules.

where:

R_miR – number of reads mapped to a particular miRNA reference in the sample,

R_all – total number of reads mapped in the sample,

RPM – normalized miRNA abundance value.

Out of the available normalization methods including DESeq [35] and TMM [36] we have chosen RPM which is a simple method that has been frequently used in current miRNA studies [37, 38]. Moreover, we normalized our data with DESeq and compared these results to RPM normalization (Additional file 1).

6 miRNA sequences representing low, medium and high abundance in HTS data analysis pipeline A were used for correlation with microarray, qPCR and HTS data from analysis pipelines B and C (Tables 1 and 2, Figure 2). Low and high abundance miRNAs were selected from the lower and upper quartile of expression, respectively. Medium abundance miRNAs stem from the interquartile range. Differences between abundance levels are significant as presented in the Table 1 and Figure 2A.

miRNA	HTS (A)		qPCR target sequence
	Abundance class	RPM
hsa-miR-27b-5p	Low1 (L1)	33	agagcuuagcugauuggugaac
hsa-miR-708-3p	Low2 (L2)	26	caacuagacugugagcuucuag
hsa-miR-30e-5p	Medium1 (M1)	3639	uguaaacauccuugacuggaag
hsa-miR-146a-5p	Medium2 (M2)	1002	ugagaacugaauuccauggguu
hsa-miR-181a-5p	High1 (H1)	28011	aacauucaacgcugucggugagu
hsa-miR-143-3p	High2 (H2)	36537	ugagaugaagcacuguagcuc

qPCR (cDNA)	***0.99	*0.90	*0.86	*0.68	**0.96
-	qPCR (library)	*0.87	*0.90	*0.73	***0.99
-	-	Array-signal	*0.87	*0.82	*0.89
-	-	-	HTS (A)	**0.93	**0.95
-	-	-	-	HTS (B)	*0.82
-	-	-	-	-	HTS (C)

*p < 0.05,

**p < 0.01,

***p < 0.001.

Pipelines B and C of HTS data analysis also start with trimmed raw sequences but count the number of reads of each small RNA isoform present in the fastq files. Their main goal is to identify isoform specific and miRNA-seed specific (functional) abundance. This counting approach allows inclusion of sequences that do not exactly match miRBase entry sequences or have an offset at the 5′ or 3′ end which could be cell type or organ specific. Moreover uncharacterized small RNAs could be detected and evaluated. For the purpose of counting a memory efficient Perl script was implemented (Sequence_counter.pl). The counting algorithm is similar to approaches proposed in the literature [35, 39], however our Perl script is designated for string operations which is more time efficient in some cases. We further analyzed our data with miRDeep2 [39] to find novel miRNAs in our data set. Because this analysis did not lead to the detection of novel miRNA species we concentrated on the evaluation of miRNA isoforms obtained from our counting algorithm (pipelines B and C). Counting introduces a huge amount of read variants that might contain sequencing errors or represent sequences of unknown origin, function and small RNA classification. Therefore, we propose a median abundance filter to remove RNA artifacts from the thyroid tissue data. Only sequences detected in at least half of the 20 samples will pass this filter and were considered in downstream analysis. Counts of sequences were also normalized to RPM values assuming:

RmiR – number of sequence counts in the sample,

Rall – total number of reads in the sample.

Pipeline B stops at this point resulting in RPM values obtained for each detected isoform. Normalized data were used for the correlation analysis of 6 miRNAs and the comparison to published FTC-FA markers.

Microarray data analysis

miRNA data were obtained from the Illumina scanner in the form of normalized log2 signal intensities. Probe sequences were re-annotated with information from miRBase version 18 according to http://www.ebi.ac.uk/arrayexpress/files/A-GEOD-8179/A-GEOD-8179.adf.txt. 1145 normalized microarray intensities were used for correlation analysis and to compare dynamic ranges. Microarray data were obtained for 20 samples analyzed with HTS, without technical replicates.

qPCR methodology

High-throughput sequencing data

We analyzed 4.7 ± 1.5 million raw sequence reads per sample (count and RPM numbers are given as mean or mean ± standard deviation) from an Illumina small RNA sequencing run using three pipelines (Figure 1). After adapter cutting 6.6*10^5 ± 4.5*10^5 reads were accepted for further analysis. Direct mapping (closest family member sequencing pipeline A) of trimmed reads to human genomic sequences using bowtie2 resulted in 9.9*10^5 ± 8.3*10^5 hits. These hits intersect with 1101 ± 153 out of 4446 genomic coordinates for miRbase entries as determined by annotation using BEDtools. Our results correspond to 91.51 ± 3.36% of all reads per sample aligned to the human genome (hg19) compared to 61.09 ± 14.45% aligned to mature miRBase sequences. Additionally, 69.23 ± 15.79% of all reads per sample mapped to hairpin miRBase sequences. The percentage of reads that comprise uncharacterized small RNA sequences is 22.3 ± 3.2%. Additional file 3 lists the ten most abundant small RNA sequences together with their genomic origin.

In isoform analysis (pipeline B) we counted individual sequence isoforms and normalized these counts to sequence abundance as reads per million (RPM). We also compared RPM to DESeq normalization (Additional file 1). Pearson correlation of data normalized with the two algorithm is 0.97. One of the 20 samples had a higher number of total reads compared to all other which resulted in a separate line in the graph (see upper diagonal line in Additional file 1). Pipeline B resulted in 18855 small, various RNA sequences that were detected in at least half of the 20 thyroid samples. Furthermore, these small sequences aligned to 502 ± 88 hairpin precursor miRBase sequences. Using miRDeep2 [39] none of the remaining uncharacterized small RNAs were identified as new miRNAs. miRBase v18 annotation was attributed to those isoform groups that contain an exact miRNA entry sequence. Seed analysis (pipeline C) merged the 18855 isoforms from pipeline B into 5189 isoform groups with identical seed sequences (i.e. the first 8 bases of each miRNA).

Comparison of high-throughput sequencing and microarray

First we compare average fluorescent intensities detected in microarray analysis with the number of reads mapped to all probe sequences of the miRNA expression BeadChip v2 (Illumina) as reference. This analysis follows the closest family member sequencing pipeline but uses microarray probe sequences and their annotation instead of human genome sequences. In Figure 3 we show the comparison for all 20 thyroid samples and 1145 individual sequences. We detect a strong saturation of microarray signal intensities compared to a wide dynamic range of the HTS data. Moreover, miRNA sequences that were not present in the HTS data show a wide spectrum of signal intensities on the array (see the most left column of dots in Figure 3). We also compared microarray signal intensities to the results of HTS pipeline B and C (Additional file 4). The saturation effect of the microarray is also seen in the plot for the pipeline B data (Additional file 4A). However, comparing microarray signal intensities to pipelines B and C data (Additional file 4) does not reveal a meaningful correlation.

Validation with qPCR

Table 1 shows the 6 mature miRNAs we selected for qPCR validation based on the level of prevalence calculated with the closest family member sequencing pipeline A. This method outputs RPM normalized read numbers after alignment to human genome sequences with bowtie2 mapping software. Our selection includes miRNAs with low, medium and high prevalence in all HTS data samples and therefore covers a wide dynamic range. miRNAs were selected from groups of prevalence at random, assuming they must be abundant in all the samples. Selected miRNAs were investigated with 2 additional high-throughput sequencing analysis approaches (pipelines B and C) using the set of 20 samples. Figure 2 summarizes analysis of HTS with pipelines A-C (2A-C), microarray (2D) and qPCR (2E and F). As seen in the direct comparison of microarray and HTS data, the dynamic range of RPM is wider than the array signal intensities. This is most prominent for the medium and high prevalence sequences (Figure 2, M1-H2). A comparison of HTS data (pipeline A, closest family member sequencing, Figure 2A) with qPCR using cDNA reverse transcribed from total RNA (Figure 2E) shows a very good concordance for 5 out of 6 miRNAs.

Correlation of small RNA qPCR and high-throughput sequencing data

To compare HTS analysis pipelines with experimental data from microarray and qPCR data we calculated Pearson correlation for all three miRNA expression profiling techniques. To detect the miRNA levels with the lowest experimental noise (which might appear e.g., due to different RT efficiencies), we used qPCR not only with the cDNAs prepared from the original RNA samples, but also with the HTS library DNAs as targets. Table 2 gives the correlation coefficients of the pair wise comparisons between methods for the expression of the selected miRNAs in 20 samples (Table 1). We observed a good correlation between the qPCR methods and microarray data with a slightly better coefficient for qPCR validation starting from reverse transcribed total RNA (r = 0.90, p < 0.05). The correlation of qPCR and HTS data strongly depends on the analysis pipeline. For both qPCR methods (cDNA and library DNA) seed analysis pipeline C (r = 0.96 and 0.99) performs much better than the closest family member pipeline A (r = 0.86 and 0.90). Overall, the highest correlation for the HTS analysis with pipeline C is detected in the comparison to qPCR validation of library DNA (r = 0.99, p < 0.001). qPCR validation starting from reverse transcribed total RNA with seed analysis pipeline C data resulted in a very similar coefficient (r = 0.96, p < 0.01). In general, qPCR and high-throughput sequencing data are highly correlated when the appropriate analysis method (seed analysis pipeline C) is applied.

Limits of qPCR validation

qPCR from cDNA template resulted in unexpected high expression values for hsa-miR-30e (Figure 2E-M1) compared to HTS analysis with pipeline A (3A). Strikingly, DNA from the sequencing libraries (Figure 2F) as template for qPCR led to results identical to cDNA. Therefore, the differences for hsa-miR-30e were not explained by the different templates, which rules out a methodological bias of the library preparation. These differences are better explained by the highly diverse isoform distribution (Figure 4) of hsa-miR-30e-5p (30e) and sequence similarities to hsa-miR-30a-5p (30a), 30c-5p (30c) and 30d-5p (30d).

Isoform sequencing pipeline B detects a total of 466 isoforms which share the same seed sequence but contain offsets at both the 5′ and the 3′ end and a high number of substitutions. The seed sequence (first 8 bases) of these isoforms matches miRBase entry sequences for different miRNAs of the hsa-miR-30 family. In principal, 30a, 30d and 30e are more related since the attribution of an individual isoform to either of these miRNAs could be based on two single nucleotide mismatches (C versus U) in the miRBase entry sequences. The distance between 30a, 30d and 30e to 30c is a little higher (4 bases). Very likely qPCR based on the target sequence for 30e will at least also partially amplify targets from 30a and 30d. Concordance of qPCR results with seed analysis pipeline C which summarizes reads from members of the hsa-miR-30 family with the same seed sequence will therefore be higher.

Evaluation of published miRNA FTC-FA markers

Table 3 summarizes sequencing results for 22 canonical mature 5p or 3p miRs listed in miRBase which originate from 12 hairpin pre-miRNA transcripts recently published as FTC – FA markers [27–29]. We tested their differential expression in our sample set of 10 FTC and 10 FA. Our HTS analysis pipelines detect 11 to 14 of the 22 miRNAs, depending on the applied analysis pipeline. When we consider a fold-change of at least 1.25 as increase/up-regulation and 0.8 as decrease/down-regulation our analysis allows a call in 9 to 12 of the 22 miRNAs (Table 3). Calls from the three HTS pipelines are concordant for 7 miRNAs. Best concordance of sequencing with published data is seen for the human hairpin miRNAs 144, 146, 197, 221 and 222. For hsa-miR-144 both 5p and 3p sequences were confirmed as differentially expressed. This result corroborates to DESeq - the most frequently used RNA-Seq expression profiling algorithm [35].

In contrast to previously published markers, the 6 miRNAs evaluated in this study do not show any significant differences between the classes of malignity as presented in the Additional file 5.

miRNA expression profiling guidelines

In a summary of the results we propose a guidance for miRNA expression profiling techniques. For experiments with aim of differential expression measurements we propose the high-throughput sequencing analysis pipeline based on seed analysis. This method will grant a functional importance of findings and the highest likelihood of qPCR validation success. In case of clinical studies, applied to samples with small cytological differences the best approach is to apply the isoform sequencing analysis pipeline. Such a method will grant the most sensitive differentiation option and most accurate use of HTS data. It also allows de novo detection of differentially expressed isoforms. Finally, for experiments with known, significant biological differences between sample classes cost of analysis per sample is an issue so microarray techniques will still be an acceptable solution. It is important to underline that among many microarray platforms only those with the newest annotations (miRBase > 17) should be considered.

Availability and implementation

Additional files 7, 8 and 9 provide raw counts for analysis pipelines A, B and C, respectively.

Data access for publication:

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39112

Software implemented for the purposes of publication: http://www.zmnieo.home.pl/autoinstalator/joomla/index.php/bioinformatics.html.