Improving Power to Detect Changes in Blood miRNA Expression by Accounting for Sources of Variability in Experimental Designs

The use of miRNAs as blood-based biomarkers is a new field of research for diagnostic and preventative medicine. A limitation of this field is the lack of statistical power to detect true differences between comparison groups, which can contribute to difficulties in validating results. Accounting for sources of variability in the experimental design may increase power in blood biomarker studies. Previously, we demonstrated that by controlling for technical variability in preparation of blood RNA for microarray analysis, we were able to improve power to detect small, yet significant, fold changes in blood transcriptomic data (1). Here, we assess sources of inter- and intraindividual and technical variability for miRNAs found in blood samples and predict how repeated measures can improve power to detect fold changes in miRNA expression.

MicroRNAs have been widely studied as biomarkers for a number of diseases. These small noncoding transcripts regulate translation of RNA by binding to the 3′ untranslated region of target RNA. Overall, miRNAs regulate 30% to 60% of RNA translation to protein usually by downregulation of the transcript (2, 3). Disease status, chemical exposures, and life-style factors have been linked to differences in expression of miRNAs observed between individuals (discussed in refs. 4–6). However, as most reported miRNA expression fold changes are small (∼1.5–2-fold), it is difficult to replicate findings and discover true associations. Therefore, it is critical that important sources of variability are predetermined and controlled for in the experimental design, a priori.

Variability in RNA transcription within subjects over time has seldom been discussed in the literature, particularly for miRNAs. Several transcriptomic studies have shown limited fluctuation in blood RNAs when measured from healthy individuals over several weeks to months (7–11). The proportion of transcripts with high intraindividual variability was attributed to a small number of immunologic genes (i.e., immunoglobulin; ref. 9) or inseparable from technical variability due to poor experimental design (10, 11). This evidence from transcriptomics suggests that there may be similarly small intraindividual variability for miRNA transcription; however, this has not been previously measured.

Other overlooked sources of variability include methods for miRNA quantification and extraction. For example, competing platforms for miRNA microarray and qPCR analysis have shown differences in sensitivity, which suggests variance in miRNA measurements (12–14). For processing of fresh blood samples, miRNAs studied in specific blood partitions (i.e., plasma, red blood cells, platelets, and leukocytes) have attributed certain miRNA expression in plasma and serum to contamination of red blood cells and platelets (15–18). Extraction of miRNAs can also introduce variability as systematic differences can depend on method or kit manufacturer (14, 19–23). Most of these previous studies focused on samples obtained from cell lines and did not thoroughly compare miRNA yield obtained from primary cells.

Here, we hypothesized that there are important sources of interindividual, intraindividual, and technical variability in miRNAs extracted from primary human peripheral blood mononuclear cells (PBMC). We calculated the contributions of these sources of variability using experimental data obtained by qPCR and compared this with estimates obtained from an miRNA microarray dataset from a previously published study. As PBMCs are a popular and noninvasive sample type and can be affected in early stages of disease, it is important to improve methods of preanalytical processing and measuring of PBMC biomarkers for future disease-related research.

Four kits were compared with each other for miRNA yield: miRNeasy kit (Qiagen), mirVana kit (Ambion/Life Technologies), ZR-duet (Zymo Research Corporation), and Trizol (Life Technologies), with the addition of 25 to 250 nmoles of c. elegans oligos spike-ins (cel-39 and cel-54) to each sample during the cell lysis step. The AllPrep kit, which comprises of on-column extraction of both DNA and RNA, was also compared with the miRNeasy Kit. Details are included in Supplementary Methods.

To calculate sources of inter- and intraindividual and technical variability, we measured miRNA expression in PBMCs of healthy volunteers over an 8-month time period. A sample size of 12 healthy subjects was included in the study (exclusion criteria for volunteer subjects were chronic illness or pregnancy at the time of blood draws). Samples were obtained at three time points, roughly 2 to 4 months apart. On the day of collection, blood samples were processed to isolate PBMCs. Data collection for this study was approved by the Internal Review Board within University of California Berkeley's Human Research Protection Program. Informed consent was obtained from all participants.

PBMCs were isolated from fresh whole blood collected in EDTA tubes using the standard Ficoll gradient protocol (24). Upon isolation of the PBMCs, they were immediately washed in PBS, pelleted, and resuspended in aliquots of RNAProtect Cell reagent (Qiagen) and frozen at −80°C until further use. At a later time, PBMC samples were thawed and RNA was extracted by the miRNeasy Kit (Qiagen) as described in Supplementary Methods.

Probe-based miRNA TaqMan Assays (Life Technologies) or PerfeCTa SYBR Green-based qScript microRNA Assays (Quanta BioSciences) were used to quantify miRNA targets of interest. Reaction volumes were proportionately scaled-down from the initial protocols (see Supplementary Methods). Six miRNAs (miR30d, let-7d, miR185, miR130a, miR451, miR342-3p) were chosen based on overlap between miRNAs expressed in the miRNA microarray dataset (used for some of the simulation studies; ref. 25) and miRNAs differentially expressed in PBMCs of Type II diabetics (26). One miRNA (cel-39) was used as an exogenous control for elution variability. Two small RNAs (SNU6 and RNU48) and miR16, frequently used for normalization of miRNA expression, were included as well.

We used a mixed-effects model to determine contributions of different sources of variability in miRNA expression for the healthy volunteer blood samples: time point, within-extraction batch, between-extraction batch, and qPCR replicate variability.

The fixed effect coefficient, β₁, represents a change in miRNA Cq value per unit change in RNA concentration for each sample, which we used as precaution against any effects that were not accounted for by using the same input (μg) of RNA for each RT-PCR reaction. The random-effects terms are defined as the change in expression (Cq value) from baseline levels for each index. Thus, Y_ijklm represents the Cq value for the i^th (i = 1…13) individual at the j^th (j = 1, 2, 3) time point, in the k^th (k = 1, 2) extraction batch, for the l^th (l = 1, 2) within-extraction replicate, for the m^th (m = 1, 2, 3) sample for a given individual on the n^th (n = 1, 2, 3) plate. Thus, β_0im is the random-effects term for the qPCR replicate for a given individual at a given time point. Lastly, the term, e_ijklmn, is defined as the “residual variability,” which includes differences from plate-to-plate and other unaccounted for sources of variability. ANOVA tests were run for each small RNA model to determine which random-effects terms were significant (P < 0.05).

We observed several measurable sources of variability for two miRNAs from the qPCR experiment, miR185 and miR451. These two miRNAs were used to determine if repeated measures would improve detection of fold changes in a theoretical study. The constraints for the theoretical study were a sample size of N = 75 versus 75 subjects in two comparison groups (e.g., disease vs. healthy controls) under two experimental designs conditions. Study 1 had no repeated measures, and Study 2 had four repeated measures for 50% of the subjects for each of the four modeled sources of variability based on our empirical qPCR data—seasonal, between-batch, within-batch, and qPCR replicate. The estimates of the parameters of interest, which is the minimum detectable fold change in the mean level of the miRNAs, are computed using the variance values attributed to the four different sources of variability from the empirical study data of 12 subjects.

Therefore, to obtain estimates of variability for the parameters of interest associated with repeating the empirical study, a clustered bootstrap procedure (using 100 bootstrap samples) was used. This procedure provides both a point estimate and confidence intervals (CI) for these fold-change estimates. Each bootstrap sample consists of data associated with 12 subjects drawn with replacement from the 12 subjects of the empirical study. The data associated with a bootstrap sample were used to estimate the variances attributed to the four sources of variability using linear mixed models (27). The distribution of the minimum detectable fold changes over the 100 bootstrap samples was used to estimate the confidence intervals of this parameter of interest.

A simulation procedure was used to estimate the parameters of interest for Study 1 and Study 2 given the variance estimates from a bootstrap sample, bs. For each miRNA within a given study design, data for 100 studies were simulated (assuming a normal distribution) using the bs variances estimates. The SD of the miRNA expression across these 100 studies provides an estimate of the standard error (SE_i,bs) of the mean level of the miRNA in a bs sample for a given study design. This SE was then used to estimate the minimum detectable fold change, FC_i,bs, with 80% power (corresponding to a 5% family-wise error rate) by the following equation:

where

To expand upon our findings of variability for individually tested miRNAs by qPCR, we examined a publicly available microarray dataset that measured hundreds of blood miRNAs simultaneously at several time points for each subject by Honda and colleagues (25). The study looked at the effects of chronic academic stress on miRNA levels in whole blood of medical students by obtaining measurements two months before, two days before, and one month after an exam for medical practitioners. The GSE49677 series from the Gene Expression Omnibus (28) Agilent-021827 Human miRNA Microarray (V3) is used in this study. The low expressed miRNAs with mean intensity levels <20 were filtered out as indicated in the study (25) and 143 miRNAs remained. The levels of these miRNAs across the four subjects and three time points are normalized using Cyclic Loess (29).

Variability estimates from the academic stress data were used to determine whether repeated measures would increase detection of differentially expressed blood miRNAs between two comparison groups (e.g., disease vs. healthy controls). We calculated estimates of interindividual variability from this study and assumed that the residual (unexplained) variability in blood miRNA levels was due to other sources. For ease of exposition, we assumed the residual to be time point (e.g., seasonal) variability, although it is likely composed of multiple sources.

For our theoretical replicate design simulation study, our sample size was 75 versus 75 subjects in each of two comparison groups (e.g., disease vs. controls). A total of 2,000 markers were evaluated for purposes of multiple testing under realistic omic-level conditions. We assumed that the sample collection for the subjects in the two groups occurred at two different time points. Therefore, seasonal effects on miRNA levels were not blocked in these experimental designs. We varied the proportion of subjects with repeated measures and the number of repeated measures per subject for each of the seven proposed designs shown in Table 1. We used a clustered bootstrap method, similar to the one described for the simulation of qPCR data, to predict minimum detectable fold changes in the mean levels for the 143 miRNAs under the seven theoretical experimental design conditions. Like before, we again provide CIs for our parameters of interest. These CIs are based on repeating the academic stress experiment (25) with four subjects. The academic stress experiment provides estimates of interindividual and residual variability. Residual variability is assumed to be attributed to seasonal variation.

We evaluated miRNA extraction procedures to find the most efficient and accurate method for our downstream applications. We presumed that lower Cq values for a given extraction method would be a proxy for both greater overall yield of all miRNAs and lower technical variability (i.e., between and within a given batch of extractions). We compared four methods: miRNeasy, miRVana, Trizol, and Zymo-Duet (which extracts both RNA and DNA). The miRNeasy Kit had the lowest Cq value for all small RNAs tested (see Supplementary Results and Supplementary Fig. S1A), which is supported by similar previous findings (14). The miRNeasy Kit also slightly outperformed the AllPrep Kit (Qiagen; Supplementary Fig. S1B) and was thus chosen as the extraction method for the downstream measurements of miRNA variability in the volunteer blood samples.

We measured several sources of variability for nine chosen miRNAs from 12 volunteer subject PBMC samples obtained at three time points over the course of several months. We included extraction replicates to account for both between- and within-extraction batch variability. Our results were analyzed using mixed-effects models. The variance attributed to each source of variability for a given miRNA is shown in the stacked bar graph (Fig. 1A and B), and the significance of each term in the model is reported in Table 2.

An ANOVA test for each of the random effects in each miRNA model indicated the significant variability terms (P < 0.05). For all endogenous and exogenous small RNAs tested, a significant proportion of variability was due to replicate qPCR reactions for a given sample. For miR16 and miR451, interindividual variability was also significant. Furthermore, for the qScript SYBR Green assays, miR451 showed a significant “time point” effect, whereas miR16 showed a significant batch-replicate effect. This is the first evidence to suggest time point and extraction variability in miRNA expression. Other transcriptomics studies have also found differences in RNA expression over the course of 1 day to several weeks (7, 8), while longer time points were inconclusive due to confounding of technical variability (10, 11).

The residual variability present for all miRNAs may be due to several sources. Covariate information was not included in this model, such as age, race, gender, and body mass index, which may contribute to the residual variability. There also may be unmeasured technical variability in the processing of samples. For example, Ficoll separation of PBMCs is not 100% efficient, so miRNAs such as miR16 and miR451, both known to be highly expressed in red blood cells (15), could contribute to variability of these miRNA expression levels. SNU6 and the exogenous miRNA, cel-39, demonstrated the highest proportions of residual variability (Fig. 1B). As expected, biologic and technical variability contributed minimally to total variability and overall variability for the spike-in cel-39, and it had the lowest total variability of all miRNAs tested. For SNU6, the source of residual variability remains unknown, but perhaps this small RNA is not ideal to use for normalization of target miRNAs in future studies if sources of variability are not representative of other miRNAs.

The remaining miRNAs did not show significant contributions of variability from the other measured sources. This may be because (i) the sample size was too small too assign statistical significance or (ii) the scaled-down volumes for the TaqMan assay were not sufficient to measure these effects. Comparing miR16 measurements in both assays, the variance term for the qPCR replicates in the model was much smaller for the SYBR Green assay than the TaqMan assay (0.119 vs. 0.499) and is plotted for each individual in Supplementary Fig. S2A and S2B. In addition, three significant random-effects terms were found with the SYBR Green assay for miR16, whereas only two were significant using the TaqMan Kit, providing further evidence that perhaps other sources of variability could be unveiled if qPCR replicate variance was reduced. A less rigorous comparison of the two assays has been made in a previous study; however, the opposite results were found (13). Our results on the performance of miRNA SYBR Green-based qPCR are supported by a very recent study that examined miRNA expression analysis of qScript and several other platforms in much greater detail (30). As many studies use TaqMan-based assays, our results could help explain the lack of reproducibility reported between studies examining miRNAs in the same tissue for the same disease in similar populations.

We used the estimates of interindividual, intraindividual, and technical variability for two of the miRNAs (miR185 and miR451) for further analysis in a theoretical study of 75 versus 75 subjects. We calculated the minimum detectable fold change with 80% power in a study with no repeated measures (Study 1) versus four repeated measures for each of the following: between-batch, within-batch, time point, and qPCR replicates for 50% of the subjects in each group (Study 2). The minimum fold-change estimate for miR451 decreased with repeated measures from 3.77 (95% CI, 1.75–5.16) to 2.38 (95% CI, 0.97–3.36), and the minimum fold-change estimate for miR185 decreased from 4.18 (95% CI, 2.67–6.46) to 2.4 (95% CI, 1.71–3.74; Table 3). For miR185, a marginally significant (90% CI) decrease in fold change was observed with repeated measures from Study 1 to Study 2. Our estimate of detectable differences in fold change for miR185 is similar to previous findings that showed a 1.82-fold difference between disease conditions measured in PBMC samples (26).

We expanded our investigation of repeated measures to examine more miRNAs simultaneously, as is currently done in omics-level studies. We used a previously published miRNA microarray dataset on four medical students over three time points to estimate the variability in 143 miRNAs (25). The interindividual variability of each miRNA from our empirical study (miR342, miR451, miR16, miR185, miR30d, let-7d, miR130a) was compared with results obtained from Honda and colleagues. There was no significant correlation between the two estimates (data not shown). This lack of correlation may be explained by differences in expression variability in each of the sample types (PBMCs vs. whole blood), as high expression for some of these miRNAs has been reported in red blood cells (15), or by the small sample sizes used to estimate variability for the qPCR data and the Honda and colleagues (25) dataset.

From the Honda and colleagues (25) dataset, we simulated data for a theoretical study to demonstrate benefits of including repeated measures. We ran the analyses with 20%, 50%, and 100% of the subjects randomly chosen for repeat sampling, and tested a total of seven different experimental designs summarized in Table 1. The cumulative distribution over the 143 miRNAs for the minimum detectable fold change across the 100 bootstrap samples is plotted for each of the seven replicate designs (Fig. 2A–C). Without repeated measures, a ≥2-fold change could be detected in approximately 24% of miRNAs. Inclusion of repeated measures for 20% of the samples improved the detection rate to 34% and 59% for Designs 1A and 1B, respectively (Fig. 2A). When repeated measures were included for 50% of the samples, the detection rate for a ≥2-fold change improved further to 46% and 69% for Designs 2A and 2B, respectively (Fig. 2B). Only a minimal increase in detection rate was gained beyond this when performing repeated measures on 100% of the samples in Designs 3A and 3B (Fig. 2C).

The estimates of 95% CIs (based on bootstrapping) for these cumulative distribution curves overlapped under each design scenario (Supplementary Fig. S3A–S3C), providing inconclusive evidence of statistically significant improvements of one design over another. To investigate this further, we compared each of the three designs with each other within the 50% repeated-measure parameter (Design 0 vs. 2A vs. 2B) and looked at a range of CIs (i.e., P values) for the minimum detectable fold change of miRNAs in one study versus another. The proportion of miRNAs with detectable fold-change differences in one design versus another at a given P value is plotted for each pair of designs (Fig. 3). From this analysis, Design 2B (with four repeated measures) shows lower detectable fold changes for approximately 20% of the miRNAs (at P < 0.10) than Design 0 (with no repeated measures). Similar comparisons for Design 1B versus 0 showed very few miRNAs with lower detectable fold changes, whereas Design 0 versus 3B showed lower detectable fold changes in approximately 40% of miRNAs (at P < 0.10; see Supplementary Results and Supplementary Fig. S4A and S4B, respectively).

Although inferences can be made from the simulations of this miRNA microarray dataset, there are still several limitations of this study. First, we assumed that the remaining variability after accounting for interindividual differences is due to intraindividual variability over time; however, technical variability from time point to time point is included too. Our estimates of intraindividual variability may be higher than expected due to the lack of technical variability measurements. Also, changes in study participants' stress levels during the collection time points (25) might also lead to overestimates of intraindividual variability. Note that to mimic realistic omic-level conditions, the theoretical experimental design for this study was limited to 75 versus 75 subjects, and we used 2,000 measured endpoint markers (including the 143 miRNAs) for purposes of multiple-hypothesis testing. Altering these parameters by including a larger (or smaller) number of subjects and/or a greater (or reduced) proportion of repeated measures would shift all three curves to the left (or right; Fig. 2). Thus, we consider the improvements of incorporating repeated measures observed herein to underestimate those expected for detecting smaller fold changes in molecular epidemiologic studies with larger sample sizes.

In summary, miRNAs have great potential as reliable blood biomarkers of early effects of the disease state. In measuring miRNA expression, we concluded that variability due to the qPCR reaction replicates generally outweighs other measured sources of variability. In the future, it would be advantageous to either troubleshoot qPCR reaction conditions to reduce this variance (e.g., increase reaction volumes, modify reaction temperatures, etc.) or increase the number of replicates per sample to account for this source of variability. Previous publications on small fold changes in blood miRNA analyzed by qPCR should be viewed with skepticism in light of this finding. In addition, methods of extraction and miRNA quantification must be rigorously tested to maximize yield and interpretable results, as performance can vary by kit and assay. For unavoidable sources of variability, block experimental designs or repeated measures should be implemented. Identifying sources of variability in future omics-level experimental designs and estimating power a priori using these described methods can save precious resources, funding, and time for molecular epidemiologic studies.

No potential conflicts of interest were disclosed.

Conception and design: S.I. Daniels, L. Zhang, M.T. Smith, R. Thomas

Development of methodology: S.I. Daniels, R. Thomas

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.I. Daniels, F.C.M. Sillé, A. Goldbaum, B. Yee, E.F. Key

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.I. Daniels, R. Thomas

Writing, review, and/or revision of the manuscript: S.I. Daniels, F.C.M. Sillé, E.F. Key, M.T. Smith, R. Thomas

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Goldbaum

Study supervision: L. Zhang, M.T. Smith, R. Thomas

This work was supported by grant # P42ES004705 from National Institute for Environmental Health Sciences Superfund Research Program (to M.T. Smith) and a fellowship from the National Science Foundation (to S.I. Daniels).