rSeqDiff: Detecting Differential Isoform Expression from RNA-Seq Data Using Hierarchical Likelihood Ratio Test

<div><p>High-throughput sequencing of transcriptomes (RNA-Seq) has recently become a powerful tool for the study of gene expression. We present rSeqDiff, an efficient algorithm for the detection of differential expression and differential splicing of genes from RNA-Seq experiments across multiple conditions. Unlike existing approaches which detect differential expression of transcripts, our approach considers three cases for each gene: 1) no differential expression, 2) differential expression without differential splicing and 3) differential splicing. We specify statistical models characterizing each of these three cases and use hierarchical likelihood ratio test for model selection. Simulation studies show that our approach achieves good power for detecting differentially expressed or differentially spliced genes. Comparisons with competing methods on two real RNA-Seq datasets demonstrate that our approach provides accurate estimates of isoform abundances and biological meaningful rankings of differentially spliced genes. The proposed approach is implemented as an R package named rSeqDiff.</p></div

Summary of notations.

<p>Summary of notations.</p

Comparisons of rSeqDiff, MATS, Cuffdiff 2 and RT-PCR assays.

<p>(A) Scatter plot of the <i>Δψ</i> values estimated by rSeqDiff (using junction reads only) and MATS. (B) Scatter plot of the <i>Δψ</i> values estimated by rSeqDiff (using junction reads only) and RT-PCR. (C) Scatter plot of the <i>Δψ</i> values estimated by rSeqDiff (using all reads) and RT-PCR. (D) Scatter plot of the log2 fold changes of isoform abundances between ESRP1 and EV estimated by rSeqDiff and Cuffdiff 2. Transcripts classified as model 0, model 1 and model 2 are shown in green, blue and red, respectively. The solid line is the regression line. The dashed line is the y = x line, which represents perfect agreement of the two methods. <i>Δψ</i>: difference of exon inclusion level between ESRP1 and EV; PCC: Pearson Correlation Coefficient; SCC: Spearman Correlation Coefficient.</p

Ranking of the RT-PCR validated genes with relevant neurological functions.

*<p>The RT-PCR result for this gene is not consistent with the exon-based method in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Voineagu1" target="_blank">[19]</a>, therefore this gene is not validated by RT-PCR. rSeqDiff classifies it in model 1.</p>**<p>FAIL: the gene has “an ill-conditioned covariance matrix or other numerical exception that prevents testing” by Cuffdiff 2 <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Website1" target="_blank">[22]</a>.</p

Illustration of the three models.

<p>(A) A hypothetical gene with three exons and two isoforms in blue and red, respectively. (B) Three models characterizing three biological situations of the gene expression patterns between two conditions. The numbers of red and blue bars represent the relative abundances of the corresponding isoforms in the two conditions.</p

The correlation coefficients of the values between RT-PCR and rSeqDiff, MATS and Cuffdiff 2 for the 164 RT-PCR tested exons^*.

*<p>The values from RT-PCR and MATS are directly adapted from <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Shen1" target="_blank">[11]</a>.</p>**<p>Three genes failed to be tested by Cuffdiff 2 (Reported as “FAIL”) are excluded.</p

Summary of hLRT for model selection.

<p>Summary of hLRT for model selection.</p

Examples demonstrating the estimates from rSeqDiff.

<p>(A)-(C) show NRCAM gene. (D)–(F) show BACE1 gene. (G)-(I) show SCIN gene. (A)(D)(G) show the gene structure and coverage of reads mapped to the gene. (B)(E)(H) show enlargement of the parts in the red boxes in (A)(D)(G), respectively, emphasizing the alternative spliced exons. In (B), the red box emphasizes the alternative exon that was validated by RT-PCR assay in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Voineagu1" target="_blank">[19]</a>, and the two red arrows represent the positions of the primers of RT-PCR <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Voineagu1" target="_blank">[19]</a>. (C)(F)(I) show estimated abundances for each gene and its isoforms by rSeqDiff. Values in the brackets are the 95% confidence intervals for the estimates.</p

Models for estimating the exon inclusion level <i>ψ</i> using the junction reads.

<p>(A) The “exon-exon junction model” used by MATS <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0079448#pone.0079448-Shen1" target="_blank">[11]</a>. Exon 1 and 3 are common exons shared by the two isoforms, and exon 2 is the skipped exon unique for the longer isoform. <i>ψ</i>: exon inclusion level; <i>UJC</i>: number of reads mapped to the upstream junction; <i>DJC</i>: number of reads mapped to the downstream junction; <i>SJC</i>: number of reads mapped to the skipping junction. (B) The “two-isoform model” transformed from (A). The abundances of the longer and shorter isoforms are <i>θ₁</i> and <i>θ₂</i>, respectively, which are estimated using the junction read counts (<i>UJC</i>, <i>DJC</i> and <i>SJC</i>).</p

Security notions of group signature schemes.

<p>Security notions of group signature schemes.</p