Normalization of RNA-sequencing data from samples with varying mRNA levels

Abstract

Methods for normalization of RNA-sequencing gene expression data commonly assume equal total expression between compared samples. In contrast, scenarios of global gene expression shifts are many and increasing. Here we compare the performance of three normalization methods when polyA(+) RNA content fluctuates significantly during zebrafish early developmental stages. As a benchmark we have used reverse transcription-quantitative PCR. The results show that reads per kilobase per million (RPKM) and trimmed mean of M-values (TMM) normalization systematically leads to biased gene expression estimates. Biological scaling normalization (BSN), designed to handle differences in total expression, showed improved accuracy compared to the two other methods in estimating transcript level dynamics. The results have implications for past and future studies using RNA-sequencing on samples with different levels of total or polyA(+) RNA.

Figure 1. Experimental design.

Total RNA from 5 developmental stages pre- and post-ZGA was isolated and kanamycin polyA⁺ RNA was used to adjust for differences in RNA yield. PolyA⁺ RNA was isolated and four cDNA libraries were generated to compare qPCR results using different template and primers.

Figure 2. Relative polyA⁺ RNA amounts.

Measurements of polyA⁺ RNA determined by a standard laboratory method (full line) and using trimmed mean of M-values (TMM) (stippled line) display an almost identical pattern during early embryogenesis with an early increase and subsequent decrease. The levels are relative to the 1-cell stage.

Figure 3. cDNA template and primer comparison.

Comparison of RT-qPCR results based on polyA⁺ and total RNA and oligo(dT) and random primers for stat3. The increase pre-ZGA is only detected in the polyA⁺ RNA-based cDNA libraries. PolyA  =  polyA⁺ RNA, Total  =  total RNA, OdT  =  oligo(dT) primers, RP  =  random primers.

Figure 4. Distribution of gene expression values.

Box plot of distribution of transcript counts or values before (not normalized) and after normalization (BSN, RPM and TMM).

Figure 5. Comparison of normalization methods.

Log2-tranformed fold-changes comparing RT-qPCR and RNA-seq data normalized using RPM, TMM and BSN for transcripts increasing pre-ZGA (a), decreasing pre-ZGA (b), decreasing post-ZGA (c) and increasing post-ZGA (d).