Modeling non-uniformity in short-read rates in RNA-Seq data

Abstract

After mapping, RNA-Seq data can be summarized by a sequence of read counts commonly modeled as Poisson variables with constant rates along each transcript, which actually fit data poorly. We suggest using variable rates for different positions, and propose two models to predict these rates based on local sequences. These models explain more than 50% of the variations and can lead to improved estimates of gene and isoform expressions for both Illumina and Applied Biosystems data.

Figure 1

Counts of reads along gene Apoe in different tissues of the Wold data. (a) Brain, (b) liver, (c) skeletal muscle. Each vertical line stands for the count of reads starting at that position. The grey lines are counts in the UTR regions and a further 100 bp. Here introns are deleted and exons are connected into a single piece. Only shown are counts on one strand of the gene; counts on the other strand show similar similarities in different tissues. Nt: nucleotides.

Figure 2

The coefficients of the Poisson linear models in different datasets. The coefficients of the Poisson linear model in the eight sub-datasets when we consider surrounding sequences as 40 nucleotides before and 40 nucleotides after the first nucleotide of a read. Position -1, 0, 1 means the nucleotide before the first nucleotide of a read, the first nucleotide of a read, and the second nucleotide of a read, respectively. Color coding for nucleotides: red, T; green, A; blue, C; black, G. The coefficients for nucleotide T (red) are the base levels, so they are always zero. (a) Coefficients in the Wold data. Shape coding for sub-datasets: rectangle, brain; triangle, liver; circle, skeletal muscle. (b) Coefficients in the Burge data. Shape coding for sub-datasets: rectangle, group 1; triangle, group 2; circle, group 3. (c) Coefficients in the Grimmond data. Shape coding for sub-datasets: rectangle, EB; triangle, ES. Following are examples of how these coefficients should be read. In the Wold brain data, the coefficient of C in the first nucleotide of a read (the blue rectangle at position 0 in (a)) is 0.82. This means that if the nucleotide T is replaced by C, then the sequencing preference will increase to e^0.82 = 2.27 times. Nt: nucleotides.

Figure 3

Fitting counts for the Apoe gene. Black vertical lines represent counts (experimental values or fitted values) along the Apoe gene (with the UTR and a further 100 nucleotides truncated). (a) Counts of reads (true values) in the Wold brain data. This is the same as the central part (black vertical lines) of Figure 1a. (b) Counts of fitted reads using the Poisson linear model. We use the other 99 genes of the top 100 genes to train the linear model, which is then used to predict the counts for Apoe. This prediction has a (cross-validation) R² = 0.54. (c) Counts of fitted reads using MART. We use the other 99 genes of the top 100 genes to train MART, which is then used to predict the counts for Apoe. This prediction has a (cross-validation) R² = 0.69.

Figure 4

Boxplot of R² for unique genes in the Wold brain data. We divided the genes with at least one read into six groups according to their RPKMs: <1, 1 to 5, 5 to 15, 15 to 30, 30 to 100, and >100; each group contains 4,205, 3,320, 2,807, 1,330, 1,094, and 383 genes, respectively. Note that in these data, 1 RPKM stands for 0.034 reads per nucleotide on average, a gene with RPKM >30 is considered to be relatively abundant, and a gene with RPKM <1 is not robust even for transcript detection [7].

Figure 5

Four isoforms of RefSeq gene Clta in mouse. This figure was generated using the CisGenome browser [36]. At the top are shown the base positions in mouse chromosome 4 and exons as grey blocks. On the bottom are shown the four isoforms, with exons zoomed in. The tail of exon 1 of the first isoform is 6 bp less than that of the other three isoforms. The second isoform has 7 exons, while the third isoform misses both exon 5 (54 bp) and exon 6 (36 bp), and the fourth isoform misses exon 6.