Input material

How well transcripts are represented in the fragment library to be sequenced may already constrain how accurate genes and isoforms can be quantified in the later steps. The number of reads generated from a given isoform depends on the isoform length, the presence of potential 5′-to-3′ biases, the presence of a poly-A tail and the presence of variations of transcription start relative to the annotated transcription start site. To overcome the dependence of the read counts on the isoform length, FPKM was introduced as expression measure [18]. FPKM builds on the assumption that the number of reads that are generated from an isoform is proportional to the isoform abundance as well as the isoform length. Additionally it normalizes for sequencing depth. A closer analysis (see e.g. [5,12]) shows that this model oversimplifies the situation since the isoform length has to be replaced by an effective length that considers the length distribution of the fragmented RNA. In our simulation we find that the relative abundance of each isoform in the starting material agrees well with the FPKM of the sequenced reads. However, if we allow for coverage bias introduced in the library preparation, then only 54% of all isoforms have an FPKM that is within a factor of two of the nominal value implied by the relative abundance and length of the isoform. Short isoforms are massively overrepresented while long isoforms are underrepresented. A less dramatic deviation occurs if the transcribed sequences are modeled to have variable transcription start sites and a poly-A tail. In this situation, 1343 isoforms are overrepresented by more than a factor of two. The length dependency of the FPKM bias induced by coverage bias and transcript modifications is visualized in Figure 1. The sharp peak of overrepresented isoforms is due to the fragmentation that generates fragments of an approximate length of 200 bp. In Figure 2a, we show the density distribution of the bias factors induced by library properties for the RefSeq genes.

Open in a separate window

Figure 1

FPKM bias induced by coverage bias and transcript modifications. The points represent RefSeq isoforms and the plot shows the length dependency of the FPKM of the sequenced reads for input material where each isoform has the same abundance and where the isoforms exhibit coverage bias or transcript modifications, respectively. The nominal FPKM implied by the isoform abundance and the nominal isoform length is indicated as horizontal line. The fragmentation step in the library preparation generates fragments of an approximate length of 200 nt and causes the sharp peak of highly overestimated FPKM for isoforms having approximately this length.

Open in a separate window

Figure 2

Biases of the observed FPKM displayed as density distributions. We visualize the effect of the library (a), alignment (b), and counting biases (c) on the observed FPKM. The log2-ratios were constrained to the range [-2, 2] and the density distributions were generated using a band-width of 0.05. The library composition has the biggest impact and causes over- as well as underrepresentation. The aligners do not align all reads correctly but achieve a very high rate for the majority of the genes. If the transcripts exhibit variable transcription start sites and poly-A tails the alignment rate for the short isoforms drops (see the Additional file 1 for the per isoform mapping rates).

We further check in real data for the presence of a transcript length effect and investigate GC bias. We use the microarray and HiSeq RNA-seq (Accession GSE37704) from Trapnell et al.[11]. Both datasets have been generated using the same samples. In Figure 3a, we show a smooth scatter plot comparing the microarray signals and the RNA-seq FPKM values. The highlighted short isoforms show that the RNA-seq data underestimates their expression relative to the microarray. If we look at isoforms that are well detected by the microarray (signal intensity>50) and low abundant in the RNA-seq (FPKM<0.01) we find that the short isoforms (length<350 nt) are overrepresented (significance in Fisher’s Exact test is 3.9e-05 with an odds ratio of 3.8). Figure 3b shows the same data but now highlighting isoforms with extreme GC content. Here, the data allows no conclusions since the microarray signals also do have a dependency on GC content. However, the presence of a GC-bias has been shown by Hansen et al. [19] revealing that isoforms with extreme GC content are underrepresented.

Open in a separate window

Figure 3

Length and GC Bias in microarray and RNA-seq data. The smooth scatterplots compare the expression scores of a single sample measured by RNA-seq and microarray. In 3a we highlight additionally the short isoforms and find that RNA-seq tends to underestimate expression of short isoforms. In 3b we highlight isoforms with extreme GC content. The plot suggests that microarray and sequencing underestimate expression of low GC isoforms, while the high GC isoforms are overestimated by microarray and underestimated by RNA-seq.

The representation biases of the isoforms together with isoform characteristics and more statistics can be found in Additional file 1.

Alignment

We consider an alignment as correct if the alignment coordinates fall within the range of the read-generating isoform. This is sufficient for the read to be counted towards the expression of the isoform. This definition implies that we tolerate if the aligner clips the read ends. We show later that read-end clipping may affect the detection of isoform switching. The medians of the per-isoform alignment rate for different methods and inputs are listed in Table 2. The highest alignment rate (with half the genes having an alignment rate of 99% or more) is achieved if the reads are simulated according to the reference transcripts. Alignment rates drop if the simulation allows for variable read starts and adds poly-A tails. Since short isoforms have proportionally more reads that include the transcript ends, they suffer most from these transcript modifications. Figure 4 shows the length dependency of the alignment rates. The difference in the alignment rates between tophat and STAR roots in the fact that STAR allows more mismatches and supports local alignments. The full table of alignment rates for the individual isoforms is available as Additional file 1. In this comparison we allowed STAR to report up to 30 valid alignments per read. This is different from the default (10). The choice was based on a preliminary study (data not shown) where we found that STAR would produce a low mapping rate for a few large gene families.

Table 2

Median value of the per-isoform mapping rate

Alignment method	Hom. coverage	Inhom. coverage	Mod. TSS+polyA
RSEM	0.990	0.990	0.898
STAR	0.943	0.943	0.858
tophat	0.814	0.814	0.737

Open in a separate window

Using the alignment methods implemented by STAR, tophat, and RSEM, and different types of input material, we compute for each isoform the mapping rate and report the median value for all isoforms.

Open in a separate window

Figure 4

Isoform length dependence of the alignment rate in the presence of variable transcript start sites and poly-A tails. If the transcripts are simulated with poly-A tails and variable transcript start sites, the mapping rate of the short isoforms is low. For the visualization at the log-scale we offset the counts by 1 and a small jitter.

Abundance estimation

We first look at gene level expression estimates, which can be derived either by a counting or an isoform abundance approach where the isoform abundances are summed up to yield the gene abundance. As Table 3 shows, the RSEM approach yields a good expression estimate for more than 92% of the genes while the counting approaches only yield an accurate result for 87 to 89% of the genes. The impact of over-/under-counting is visualized in Figure 5. In the case of human genome 3-4% of the genes are strongly affected. Figure 5 also shows that RSEM’s posterior estimation is able to correct some of the reads that are initially incorrectly assigned. For reference we report all gene level scores in Additional file 2.

Table 3

Fraction of genes where the expression count bias is less than a factor of two

Alignment and count method	Hom. coverage	Inhom. coverage	Mod. TSS+polyA
STAR-countOverlaps	0.892	0.892	0.892
tophat-countOverlaps	0.893	0.893	0.893
STAR-htseq	0.877	0.877	0.877
tophat-htseq	0.883	0.883	0.883
RSEM-count	0.929	0.929	0.929
RSEM-posteriorCount	0.930	0.930	0.930

Open in a separate window

For each gene we estimate expression in terms of read count and report the fraction of genes for which the estimate is less than a factor of two different from the number of correctly aligned reads. We consider the alignments produced by STAR and tophat in combination with subsequent counting of the overlapping reads using HTSeq and the countOverlaps method. The RSEM software performs alignment and expression counting. Here we list the results for the reported counts as well as the results for the posterior estimates of the counts.

Open in a separate window

Figure 5

Comparison of reported expression counts with the number of reads sequenced. The biases are dominated by the way the reads with multiple alignments are taken into account. Multiple counting leads to overestimation while ignoring multi-reads leads to underestimation. The most accurate results are achieved by the posterior estimates of RSEM.

We show the density distributions of the count biases in Figure 2c. The comparison with the library representation and the read alignment biases shows that the counting biases affect the fewest genes.

Additionally, we verify how well RSEM can reproduce the isoform abundances. We find that RSEM does not correctly produce zero counts in the case where a gene has multiple isoforms but only one isoform is expressed. This observation is also reported by Mezlini et al.[20]. It is caused by the fact that RSEM tends to bias the expression estimates of the individual isoforms towards the mean value of the corresponding gene locus.

Differential expression

We further generate samples that exhibit differential expression and analyze the ability to assess the differentially expressed genes. To this end, we generate different conditions with three replicates each. The conditions have 10% of the genes upregulated by a factor of two and 10% of the genes downregulated by a factor of two relative to the baseline condition. We simulate biological noise by adding a log-normal noise to the expression fractions of the replicate samples that are used by Flux as input. In total, we simulate 11 different conditions, and test subsequently for differential expression between pairs of conditions. Based on the set of all comparisons, we assess the receiver-operator-characteristic (ROC) for the different test-methods and expression estimates. Table 4 shows the Area-under-the-curve (AUC) values for the exactTest, log-ratio, and t-test applied to different expression estimates. We get AUC above 0.91 when using the exactTest independent of the expression estimate. Interestingly, using the effect-size, i.e. the mean log-ratio as score for differential expression, the AUC performance is even a minute amount higher, again independent of the expression estimate used. Additionally we have included the t-test applied to the log-expression scores. The assumptions of the t-test are not compatible with the count data and the t-test does not make use of the mean-variance relationship that can be found in counting data. We show it here to demonstrate the benefit obtained by using an appropriate count model (as in edgeR). Figure 6 shows the ROC graphs when assessing differential expression using exactTest and using different expression estimates as input. We find that if one uses HTSeq counts as input, the maximum achievable sensitivity is lower than for the other methods. This is due to the fact that for overlapping gene symbols HTSeq does not count the reads and for those genes differential expression cannot be positive. For each combination of counting and testing method, we have computed the per-gene AUC values and report them in Additional file 2. We consider a gene as not detectable by a given method if its AUC value is below 0.1. Table 5 shows that the number of genes not detectable by an individual method ranges from 0 to 506 where the highest number of undetectable genes is observed when using HTSeq counts. In principle, all genes for which reads have been produced, can be detected as differentially expressed if the countOverlaps approach or the RSEM posterior counts are used.

Open in a separate window

Figure 6

Receiver-operator-characteristics for the detection of differential expression based on p-values computed with edgeR’s exactTest. The Area-under-the-curve (AUC) depends on the type of preprocessing. Highest AUC is achieved with RSEM counts reported by RSEM. Other methods score slightly lower with htseq-based counts reaching a plateau that is below 1.0 due to the fact that htseq would not report counts for overlapping gene symbols.

The comparison demonstrates the power and benefit of the counting model but also shows that if only few replicates--in our case three--are available, the effect-size is a valuable score for differential expression. Corresponding findings were also reported for microarray data [21]. In a related work on comparing hypothesis test based differential expression of counting data, Robles et al. [22] have shown that with 6 replicates there is a remarkable improvement of the power. This power comes mostly from the improved precision on the variance estimate. For RSEM, the distributions of the read counts do not follow strictly a Poisson or negative-binomial model, since the assignment of multi-alignments leads to fractional read counts and the posterior estimation adds additional noise to the estimate. However, we find that the Poisson approximation is reasonable and leads to accurate results.

We additionally evaluate to which extent isoform switching without expression change affects the assessment of the differential expression estimate. We use again conditions with three replicates where only the isoform changes but not the number of transcripts expressed from a given locus. Figure 7 shows that such isoform changes lead to an increased false positive rate (FPR). Isoform switching can lead to significant change in the isoform length and consequently a change in the read count, which is then reported as expression change. The only way to circumvent this is to resort to expression measures that are normalized for isoform length, like FPKM or transcripts per million transcripts sequenced (TPM).

Open in a separate window

Figure 7

Isoform switching increases the False Positive Rate (FPR) of differential expression. If at a gene locus the expressed isoform changes, this leads to an increase in the false positive detection of differential expression based on gene-level counts. For set of significance thresholds we plot the FPRs in the absence and presence of isoform switching. False Positive Detection can only be avoided if expression counts are normalized for isoform length as is the Transcripts-per-Million (TPM) score returned by RSEM.

Read simulation with flux

For simulation we use the Flux read simulator v1.2. We simulate coverage bias by setting the fragmentation substrate to DNA. We choose the default options for the transcript modification to obtain transcripts with variable read starts and poly-A tail. We set both options to NaN and the fragmentation substrate to RNA to produce transcripts that follow exactly the annotated transcripts. We choose as read length 101 bp and use the default error distribution provided by Flux.

Read alignment

For read alignment, we use tophat v2.0.6 and provided tophat with the RefSeq transcript coordinates as a GTF file so that tophat would also align the reads to the known transcripts. Additionally, we used STAR v2.3.0e with the options“--outSAMstrandField intronMotif --outFilterMatchNmin 20--outFilterMismatchNmax 5 --outFilterMismatchNoverLmax 0.05--outFilterMultimapNmax 30″ that add additional strand information for spliced alignments, limit the number of mismatches to 5 and that increase the number multimapping for reads to 30. We use RSEM v1.2.3 for read alignment to the transcriptome and use RSEM also to compute expression estimates.

Expression estimation

We compute expression with HTSeq v0.5.3p9 using the script HTS.scripts.count, as well as by using the countOverlaps method in the Bioconductor package GenomicRanges v1.10.7 and the expression output generated by RSEM.

Differential expression

We compute differential expression with the exactTest method in the Bioconductor package edgeR v3.0.6, and the t-test in the package genefilter v1.40.0. Additionally we compute the log-ratio as an effect size. We did not apply any normalization with respect to sequencing depth, since all samples had the same sequencing depth.