Assessment of microRNA differential expression and detection in multiplexed small RNA sequencing data

Abstract

Small RNA sequencing can be used to gain an unprecedented amount of detail into the microRNA transcriptome. The relatively high cost and low throughput of sequencing bases technologies can potentially be offset by the use of multiplexing. However, multiplexing involves a trade-off between increased number of sequenced samples and reduced number of reads per sample (i.e., lower depth of coverage). To assess the effect of different sequencing depths owing to multiplexing on microRNA differential expression and detection, we sequenced the small RNA of lung tissue samples collected in a clinical setting by multiplexing one, three, six, nine, or 12 samples per lane using the Illumina HiSeq 2000. As expected, the numbers of reads obtained per sample decreased as the number of samples in a multiplex increased. Furthermore, after normalization, replicate samples included in distinct multiplexes were highly correlated (R > 0.97). When detecting differential microRNA expression between groups of samples, microRNAs with average expression >1 reads per million (RPM) had reproducible fold change estimates (signal to noise) independent of the degree of multiplexing. The number of microRNAs detected was strongly correlated with the log2 number of reads aligning to microRNA loci (R = 0.96). However, most additional microRNAs detected in samples with greater sequencing depth were in the range of expression which had lower fold change reproducibility. These findings elucidate the trade-off between increasing the number of samples in a multiplex with decreasing sequencing depth and will aid in the design of large-scale clinical studies exploring microRNA expression and its role in disease.

Keywords: coverage; detection; differential expression; microRNA; multiplexing; sequencing.

FIGURE 1.

Sequencing of multiplexed small RNA samples. (A) Using the Illumina TruSeq kit, lung tissue samples were given a unique index, placed in pools consisting of one, three, six, nine, or 12 samples, and sequenced using the Illumina HiSeq 2000. The same six samples were sequenced in the 6-plex, 9-plex, and 12-plex (gray) and used in differential expression analysis. (B) As expected, the total number of reads per sample decreased as more samples were included in the multiplex. (C) Percentage of aligned reads and (D) percentage of reads that aligned with a mismatch were not significantly associated with numbers of samples within a lane indicating that multiplexing does not affect read quality.

FIGURE 2.

Assessment of microRNA expression across different multiplexes. (A) All pairwise Pearson correlations were calculated between the samples in the 1-plex, 3-plex, 6-plex, 9-plex, and 12-plex. The correlation for sample B between the 3-plex and the 12-plex is shown as an example (R = 0.98). (B) For each sample sequenced in different multiplexes, we performed a test based on the hypergeometric distribution to determine if the proportions of reads for an individual microRNA were significantly different for a replicate sample sequenced in two multiplexes than what would be expected by chance. The x-axis shows the average RPM expression for each microRNA across the 3-plex and the 12-plex. The y-axis shows the difference in RPM expression for each microRNA between the 3-plex and the 12-plex. Positive and negative values on the y-axis indicate that the RPM values in the 12-plex were smaller or larger than RPM values in the 3-plex, respectively. Red indicates that the microRNA had a significant difference in the proportions of reads between multiplexes (FDR q < 0.05). (C) The distribution of Pearson correlation coefficients was significantly higher between replicate samples sequenced in different multiplexes (representing technical variability) compared with the correlations between different samples sequencing within the same multiplex (representing biological variability) and between different samples sequenced in different multiplexes (representing technical + biological variability). (D) Likewise, the numbers of microRNAs with a significant difference in the proportion of reads were lower when comparing replicate samples in different multiplexes than when comparing different biological samples. Asterisk indicates P < 0.001 from a Wilcoxon rank-sum test.

FIGURE 3.

Effect of multiplexing on fold change reproducibility. (A) Fold changes between all possible combinations of three versus three samples (n = 10) were calculated within the 6-plex and 12-plex using DESeq (Anders and Huber 2010). In addition, the interaction effect between sample class and multiplex was calculated which determines the degree to which fold changes in the 6-plex and 12-plex are discordant. (B) The difference between the fold change estimates between the two levels of multiplexing is plotted as a function of average expression level. While only nine microRNAs displayed a significant interaction effect on any comparison (black; P < 0.05), the microRNAs with lower average expression tended to show the greatest difference in fold change estimates between the multiplexes. (C) The correlation between fold changes of different multiplexes using all microRNAs was R = 0.51. However, when microRNAs with lower average expression were iteratively removed, the correlation between fold changes rapidly improved. (D) The fold change estimates in the 12-plex are plotted against the fold change estimates in the 6-plex for all microRNAs. MicroRNAs with average RPM > 1 across both plexes are colored black while those with average RPM < 1 across both plexes are gray. These results suggest that the fold change estimates for microRNAs with average RPM > 1 are largely reproducible when the samples are sequenced to greater depths.

FIGURE 4.

Clustering of replicate sequenced in different multiplexes. MicroRNAs with an average RPM > 1 were used to perform hierarchical clustering with average linkage. Replicate samples sequenced in different multiplexes clustered together rather than different biological samples from the same multiplex.

FIGURE 5.

Detection of annotated microRNAs in samples across different levels of multiplexing. (A) The number of microRNAs detected with at least one read is correlated with the log₂ number of aligned reads (R = 0.96) suggesting an exponential relationship. (B) The number of microRNAs detected in each multiplex is shown for microRNAs with average RPM < 1 (left) and average RPM > 1 (right). When sequencing depth is decreased by multiplexing more samples, the number of detected microRNAs with RPM < 1 decreased while the number of detected microRNAs with RPM > 1 remained relatively constant.

FIGURE 6.

Examining the trade-off between having more samples at lower sequencing depth versus having fewer samples with higher sequencing depth. Fold changes ranging from 1.25 to 25 were spiked into randomly selected microRNAs in sample sets of different sizes with different levels of sequencing depth. Sets with more samples had fewer numbers of reads per sample on average. MicroRNAs were binned according to average expression across all samples in the 9-plex. The x-axis is the median log₂ RPM for all microRNAs within a bin. The y-axis is the smallest fold change that had median sensitivity of at least 80% across all microRNAs within that bin. These results indicate that, for microRNAs with higher average expression (RPM > 1), larger samples sizes are better for detecting smaller fold changes compared with smaller sample sizes.