Differential analysis of gene regulation at transcript resolution with RNA-seq

Abstract

Differential analysis of gene and transcript expression using high-throughput RNA sequencing (RNA-seq) is complicated by several sources of measurement variability and poses numerous statistical challenges. We present Cuffdiff 2, an algorithm that estimates expression at transcript-level resolution and controls for variability evident across replicate libraries. Cuffdiff 2 robustly identifies differentially expressed transcripts and genes and reveals differential splicing and promoter-preference changes. We demonstrate the accuracy of our approach through differential analysis of lung fibroblasts in response to loss of the developmental transcription factor HOXA1, which we show is required for lung fibroblast and HeLa cell cycle progression. Loss of HOXA1 results in significant expression level changes in thousands of individual transcripts, along with isoform switching events in key regulators of the cell cycle. Cuffdiff 2 performs robust differential analysis in RNA-seq experiments at transcript resolution, revealing a layer of regulation not readily observable with other high-throughput technologies.

Figure 1

Changes in fragment count for a gene does not necessarily equal a change in expression. (a) Simple read-counting schemes sum the fragments incident on a gene's exons. The exon-union model counts reads falling on any of a gene's exons, whereas the exon-intersection model counts only reads on constitutive exons. (b) Both of the exon-union and exon-intersection counting schemes may incorrectly estimate a change in expression in genes with multiple isoforms. The true expression is estimated by the sum of the length-normalized isoform read counts. The discrepancy between a change in the union or intersection count and a change in gene expression is driven by a change in the abundance of the isoforms with respect to one another. In the top row, the gene generates the same number of reads in conditions A and B, but in condition B, all of the reads come from the shorter of the two isoforms, and thus the true expression for the gene is higher in condition B. The intersection count scheme underestimates the true change in gene expression, and the union scheme fails to detect the change entirely. In the middle row, the intersection count fails to detect a change driven by a shift in the dominant isoform for the gene. The union scheme detects a shift in the wrong direction. In the bottom row, the gene's expression is constant, but the isoforms undergo a complete switch between conditions A and B. Both simplified counting schemes register a change in count that does not reflect a change in gene expression.

Figure 2

An overview of the Cuffdiff 2 approach to isoform-level differential analysis of RNA-seq data. (1) The variability in fragment count for each gene across replicates is modeled. (2) The fragment count for each isoform is estimated in each replicate, along with (3) a measure of uncertainty in this estimate arising from ambiguously mapped reads, which are extremely prevalent in alternatively spliced transcriptomes. (4) The algorithm combines estimates of uncertainty and cross-replicate variability under a beta negative binomial model of fragment count variability to estimate count variances for each transcript in each library. (5) These variance estimates are used during statistical testing to report significantly differentially expressed genes and transcripts.

Figure 3

Comparison of Cuffdiff 2 with other expression platforms. (a) Fold changes in multi-isoform gene expression measured by RNA-seq and microarrays before and after HOXA1 knockdown are highly concordant (Spearman correlation = 0.86). (b) Computing gene expression by isoform deconvolution instead of by fold change in gene-level fragment counts improves agreement with microarrays. Genes shown are those where Cuffdiff 2 and intersection-count fold changes are most discrepant (1% tails). (c) Methods for performing differential gene expression analysis with RNA-seq based on the exon-union counting method returned nearly all of the genes returned by Cuffdiff 2. (d) Lung fibroblasts show lower splicing complexity than human ESCs. Complexity is measured by the Shannon-Entropy of the relative isoform abundances of multi-isoform genes (Methods), where zero complexity indicates a gene has only a single detectable isoform. A gene with equally abundant isoforms has maximal complexity. hLF, human lung fibroblast. (e) The contribution of the major isoform of each gene to total gene expression in lung fibroblasts and ESCs.

Figure 4

Accuracy of Cuffdiff 2 over varied experimental designs. (a) Accuracy of the proposed model explored through simulated RNA-seq. Read-length series was generated with single-ended sequencing data. FPKM, fragments per kilobase per million fragments mapped. (b) Squared coefficient of variation versus expression for genes and transcripts as measured by the HiSeq 2000 and MiSeq. Each series is a fit of a generalized additive model to the individual expression, squared coefficient of variation pairs. (c) Significant gene lists returned by Cuffdiff 2 using the HiSeq 2000 and MiSeq compared against expression microarrays (FDR ≤ 1%). (d) Significant isoform lists returned by Cuffdiff 2 when using the HiSeq 2000 and MiSeq (FDR ≤ 1%).

Figure 5

Changes in expression of cell cycle regulatory genes in response to HOXA1 knockdown. (a) GSEA analysis of the knockdown for selected REACTOME gene sets. (b) Cuffdiff 2 reports an increase in CDK2 expression, which is attributable to a single isoform that includes the full activation loop, a feature required for maximal CDK2 activity. KD, knockdown. (c) Cuffdiff 2 reports a decrease in ORC6 expression, which is attributable to a single isoform that includes the full suite of residues required for optimal DNA binding not present in the minor isoforms arising from the gene. (d) Cuffdiff 2 reports an increase in TBX3 attributable to a single isoform lacking an exon situated within the T-box DNA binding domain that is present in a highly similar minor isoform. (e) Cuffdiff 2 reports a decrease in CDC14B attributable to decreases in the two major isoforms. Error bars indicate 95% confidence intervals in expression. (f) Changes in isoform expression reported by Cuffdiff 2 compared against measurements made with isoform-specific qPCR. The black line indicates perfect correspondence between the two platforms. The orange line is a linear regression through all points, and the red line excludes the three major outliers, which target low abundance isoforms, two of which cannot be distinguished from primary transcript or genomic DNA. DE, differentially expressed.

Figure 6

Cell cycle analysis after HOXA1 knockdown. (a) Human lung fibroblasts transfected with scrambled siRNAs and a HOXA1 siRNA pool 48 h after transfection. HeLa cells transfected with scrambled siRNAs (left) and a HOXA1 siRNA pool 72 h after transfection (right). Scale bars, 500 μm. (b) HOXA1 siRNAs disrupt normal cell cycle distribution. Top, histograms of Hoechst-33342 staining in human lung fibroblasts transfected with scrambled siRNAs and a HOXA1 siRNA pool 48 h after transfection (right). Bottom, same as above using HeLa cells 72 h after transfection. Cell cycle phase gates were drawn as approximations of the Dean-Jett-Fox cell cycle modeling algorithm. (c) HOXA1 siRNAs shift cell cycle distribution toward the G1 phase and sub-G1 phase fraction. Percentage of cells in each phase of the cell cycle (average three experiments) using gates from b with scrambled and HOXA1 siRNAs in human lung fibroblasts and in HeLa cells. (d) Representative dot plot of Annexin V (x-axis) versus propidium iodide (PI) (y-axis) analyses of scrambled (left) and HOXA1 siRNA (right) transfected lung fibroblasts (top) and HeLa cells (bottom). Cells populating the upper (late apoptosis) and lower (early apoptosis) right quadrants were designated as apoptotic cells. (e) Percentage of apoptotic cells (average of three experiments) populating upper and lower right quadrants from lung fibroblasts and HeLa cells transfected with scrambled or HOXA1 siRNAs as shown in d. Error bars indicate 95% confidence intervals in expression.