SASeq: A Selective and Adaptive Shrinkage Approach to Identify and Quantify Condition-Specific
Transcripts using RNA-Seq
Tin Chi Nguyen 1 , Nan Deng 1 and Dongxiao Zhu 1,∗

arXiv:1208.3619v2 [q-bio.QM] 23 Feb 2013

1

Department of Computer Science, Wayne State University, 5057 Woodward Avenue, Detroit, MI 48202, USA.
∗
To whom correspondence should be addressed.

Abstract

tially improved both the data quantity and quality to better
identify and quantify condition-specific transcripts. Existing approaches are either based on statistical modeling of read counts,
such as Poisson [19, 20, 21, 22, 23], Gaussian [24], Multinomial [25, 26], and Dirichlet-Multinomial [27], or are based on
mathematical modeling of the relationships between per-base
exonic coverage signal and transcript structures [28, 29]. Despite initial success, significant challenges remain in model overfitting. For example, there are around 130, 000 reference transcripts corresponding to some 20, 000 annotated human proteincoding genes in the Ensembl database (version 65). 59% of these
genes have 4 or more annotated transcripts. However, only 1 major transcript isoform and 1 or 2 minor transcript isoforms are
typically expressed under a single biological condition, making
the total number of condition-specific transcripts per gene less
than 4. Thus, model overfitting is inevitable for the majority of
genes and appropriate model selection is urgently needed.

Identification and quantification of condition-specific transcripts
using RNA-Seq is vital in transcriptomics research. While initial
efforts using mathematical or statistical modeling of read counts
or per-base exonic signal have been successful, they may suffer
from model overfitting since not all the reference transcripts in
a database are expressed under a specific biological condition.
Standard shrinkage approaches, such as Lasso, shrink all the
transcript abundances to zero in a non-discriminative manner.
Thus it does not necessarily yield the set of condition-specific
transcripts. Informed shrinkage approaches, using the observed
exonic coverage signal, are thus desirable. Motivated by ubiquitous uncovered exonic regions in RNA-Seq data, termed as
“naked exons”, we propose a new computational approach that
first filters out the reference transcripts not supported by splicing
and paired-end reads, then followed by fitting a new mathematical model of per-base exonic coverage signal and the underlying
transcripts structure. We introduce a tuning parameter to penalize the specific regions of the selected transcripts that were
not supported by the naked exons. Our approach compares favorably with the selected competing methods in terms of both
time complexity and accuracy using simulated and real-world
data. Our method is implemented in SAMMate, a GUI software
suite freely available from sammate.sourceforge.net.

1

Existing shrinkage approaches, such as IsoLasso [30] and
SLIDE [31], identify condition-specific transcripts by shrinking the abundance proportion parameters associated with each
transcript isoform towards zero using the standard Lasso procedure [32]. They represent the first efforts to address model
overfitting issues in identifying and quantifying condition-specific
transcripts. Although their statistical models differ in details,
they both divide the coding sequence into sub-exons. Read
counts falling into these sub-exons are then used to estimate the
transcript proportions. The approaches are limited since thousands of per-base signal are interpolated to only a dozen of data
points, which are then used to estimate transcript expression.
Furthermore, we argue that a straightforward application of a
Lasso shrinkage approach for identifying condition-specific transcripts may be limited for a variety of reasons. First, choosing
the tuning parameter for the Lasso shrinkage is not straightforward since the number of expressed transcripts is unknown and
varies from gene to gene and from condition to condition. Second, the Lasso approach shrinks all the component coefficients,
i.e. transcript abundance parameters, in a non-discriminative
way and does not necessarily lead to the set of expressed transcripts [32, 33, 34, 35].

Introduction

Identification and quantification of active transcripts using gene
expression data is essential to a wide range of transcriptomics
research. The problem is non-trivial since the observed exonic expression signal may be aggregated from a set of expressed transcripts encoded by the same gene via alternative
splicing mechanisms. Moreover, excessive sequencing errors and
biases exist in the data making the problem even more challenging [1, 2, 3, 4]. In earlier studies, several computational
approaches have been developed to utilize high throughput gene
expression profiling data collected from microarray experiments
such as an Expectation-Maximization (EM) algorithm using
Expressed Sequence Tags (ESTs) [5], and a Nonnegative Matrix Decomposition (NMF) based algorithm [6] using exon and
exon-exon junction microarrays. These approaches represent
the pioneering efforts to tackle the isoform quantification problem. However, their performances were limited by data quantity
(ESTs) and quality (microarrays).
With the widespread proliferation of next generation sequencing, unprecedented opportunities now exist for both genomics [7, 8, 9, 10, 11] and transcriptomics research [12, 13, 14,
15, 16, 17, 18]. Particularly, RNA-Seq technology has substan-

Here we present a new computational approach that exploits
both read counts and per-base exonic coverage signal, which
carry alignment information from different aspects. Therefore,
combining them is expected to yield a more accurate estimation of expression abundance. First, we use the reads that are
compatible for each isoform to filter out potentially inactive transcript isoforms. Second, we develop a new Selective and Adaptive Shrinkage approach (SASeq) to address model overfitting
issues using per-base exonic coverage signal. The rationale of
our informed shrinkage approach is that transcripts with large
1

A paired-end read is compatible with an isoform as long as both
reads are compatible. As opposed to single-end reads, paired-end
reads usually span a longer range and can thus help to identify
longer splicing patterns. From the read-isoform compatibility
information, existing read counts based approaches [19, 20, 21,
22, 23, 24, 25, 26, 27, 30, 31] typically calculate the read counts
falling into sub-exons and then infer the isoform proportions via
iterative optimization algorithms. As the number of isoforms
increases, the read counts falling into the sub-exons might be
insufficient for transcriptome quantification purposes. Therefore,
we use read-isoform compatibility information only to filter out
inactive isoforms as opposed to calculating an exact numerical
value for each isoform proportion.
As the first pre-processing step, we count the number of compatible reads (single-end or paired-end) for each isoform. Any
isoform that has an potential expression abundance lower than
a cutoff is filtered out. The potential expression of an isoform is
calculated as the number of its compatible reads multiplied by
the read length divided by the isoform length. We set the cutoff
as low as 1% of the gene’s average coverage to avoid filtering
out active isoforms. For example, in Figure 2A, the 4th isoform
will be filtered out because it is not supported by any paired-end
read.
As the second pre-processing step, we augment the per-base
exonic coverage signal by assuming a paired-end read originates
from a cDNA fragment. Our goal is to reconstruct the original cDNA fragment sequences and use them to estimate the
expression abundance of each transcript. It is because a transcript abundance is more proportional to the total bases of its
cDNA fragment than to the total bases of its sequenced reads.
As such, when a pair of short reads overlap, the source cDNA
fragment sequence is reconstructed by removing the overlapping
signal. Otherwise, recovering the cDNA fragment sequence may
be more complicated depending on the structure of the compatible isoforms.
Figure 2B and 2C present a showcase example in which two
cDNA fragment sequences can be unambiguously determined.
The 1st and 2nd paired-end reads (from the top) must skip the
2nd exon and span across the other three exons because they
belong to the 2nd transcript. In this case, we extend the short
reads along the exonic regions to recover their cDNA fragment
sequence. The 3rd paired-end read may be from any of the three
transcripts. Thus, the original cDNA fragment may or may not
cover the exons in the middle. In this case, we just use the signal
observed from the reads and do not augment signal. After this
step, the remaining isoforms and the augmented per-base signal
will serve as input for our shrinkage approach, which will be
described in the following sections.

Prevalence of Naked Exons

5000
4000
3000
1000

2000

Gene Number

6000

7000

Genes With Naked Exons
All Genes

10

1

10

2

10

3

4

10

10

5

6

10

Read Number

Figure 1: Naked exons is prevalently observed in RNA-Seq read alignment results. The histogram’s horizonal axis bins each gene based on
the number of mapped reads. The histogram’s vertical axis displays
the number of genes falling into each bin. Blue corresponds to the
total number of genes falling into each bin whereas green corresponds
to a subset of genes with at least 10% naked exons. The presented
data set consists of 80 million reads with length 75. 20, 617 out of the
total 23, 511 expressed genes (88%) have naked exons.

naked exons in a highly expressed gene are likely to be inactive. Further we note the naked exons are prevalently observed
in real-world RNA-Seq read alignment results, such as the one
shown in Figure 1. Instead of shrinking all of the transcript
abundance parameters, only selected regions in selected transcripts that are not supported by the naked exons are penalized.
We employ a tuning parameter to penalize those regions. The
tuning parameter represents the shrinkage level that is adjusted
adaptively proportional to the exonic coverage via an optimization algorithm. We implement the new computational approach
in our GUI software suite, SAMMate, which is freely available
from sammate.sourceforge.net.

2

Methods

RNA-Seq reads can be aligned using many tools including
TopHat [36], SplitSeek [37], MapSplice [38], ABMapper [39],
and SpliceMap [40]. The alignment results are saved in SAM
or BAM format that can be used as input for SASeq. For each
gene, we extract the read-isoform compatibility information and
also build a vector of per-base exonic signal from the alignment
results.

2.1

2.2

An Informed Shrinkage Motivated by Naked
Exons

We have recently discovered a ubiquitous phenomenon that a
large portion of exonic regions have no coverage signal. We denote these uncovered exonic regions as naked exons. We illustrate this phenomenon in Figure 1 using a real-world RNA-Seq
data set, where the transcriptome of a liver tissue sample from
a 37 year-old Caucasian male was sequenced (NCBI SRA accession number ERX011211). We binarize the genes according
to the number of mapped reads. In each bin, we plotted both
the total number of genes (blue) and genes with at least 10%
naked exons (green). Clearly, a vast majority of genes (88%),

Pre-processing: Filtering and Signal Augmentation

Each gene consists of a set of annotated exons whereas each
transcript isoform consists of a subset of those exons. A singleend read is compatible with an isoform if it maps to one exon
or if it contiguously spans across multiple exons of the isoform.
2

Figure 2: An example of transcript filtering and signal augmentation. The upper panel (A) shows an example of how read counts information
can help to filter out inactive isoforms. The isoforms are in blue and red whereas the reads are in green. The junctions of a read are connected
by the green lines while the read pairs are connected by blue dash lines. The 4th isoform is not supported by any read and thus will be excluded.
The lower panels (B, C) show an example of how to recover the original fragments where paired-end mates are separated by a gap. (B) The
1st and 2nd paired-end reads (from the top) is compatible with the 2nd isoform and must skip the 2nd exon while spanning the 1st , 3rd and
4th exons. The 3rd paired-end read may originate from any of the three transcripts. (C) The exonic coverage signal is augmented for the 1st
and 2nd paired-end reads, but not for the 3rd paired-end reads since we can not determine whether the cDNA fragment covers the 2nd and 3rd
exons.

At each exonic position
i, the expected per-base exonic exP
pression signal is r N
S
pj whereas the observed exonic exij
j=1
pression signal is ei . The latter is accumulated from both exonic
and junction read coverage. In addition to exonic reads, junction
reads are indispensable for identifying and quantifying adjacent
exons as they cover the exon-exon junction regions. They are
also essential for identifying short exons whose expression signal
may be only identifiable by junction reads. Moreover, as the read
length increases, a read is more likely to span multiple exons and
hence carries valuable coverage information for our calculation.
We write the observed per-base expression signal as a sum of the
model explained portion and error as follows:

including highly expressed genes, possess naked exons. The similar trend was also observed for many other RNA-Seq data sets,
generated using different Illumina platforms, from NCBI SRA.
Motivated by the ubiquitous naked exons, we developed an informed shrinkage approach that shrinks the proportions of those
isoforms not supported by the naked exons, and preserve active
isoforms from the genes with low expression. The rationale is
that isoforms with naked exons in a highly expressed gene are
more likely to be inactive than those without.
2.2.1

Gene Structure Model with Informed Shrinkage

For each gene, we denote the number of bases in the exonic
regions by M , the number of transcripts by N and the genelevel expression abundance by r. Thus, the vector of observed
per-base signal is given by e = [e1 , e2 , ..., eM ]T , where ei is the
observed read coverage at the ith exonic position. We also denote the vector of transcript proportions by p = [p1 , p2 , ..., pN ]T ,
where pj is the abundance proportion of the j th transcript,
P
N
j=1 pj = 1 and 0 ≤ pj ≤ 1 for all j. The transcript isoform structures are represented using a splicing matrix, S ∈
{0, 1}M ×N of 0’s and 1’s. Each row i represents a single base,
and each column j represents a transcript isoform. Sij =1 if the
j th transcript isoform contributes to the exonic signal observed
at the ith base and Sij =0 otherwise.
Figure 3A illustrates an example gene model without shrinkage. In this example, the 1st transcript is full-length. Thus,
all of the elements of the 1st column of S take the value of 1
(Sj1 = 1 for all j ∈ [1 . . . M ]). In the 2nd column of S, the elements corresponding to the 2nd exon take the value of 0 because
the 2nd transcript skips the 2nd exon. In the 3rd column of S,
the elements corresponding to the 3rd exon take the value of 0
because the 3rd transcript skips the 3rd exon.


e1
 e2 


 ... 
eM


=

=


r(S11 p1 + S12 p2 + . . . + S1N pN )
 r(S21 p1 + S22 p2 + . . . + S2N pN ) 
+ε



...
r(SM 1 p1 + SM 2 p2 + . . . + SM N pN )

 

p1
S11
S12 . . . S1N


 S21
S22 . . . S2N 
 ×  p2  + ε,
r
  ... 
 ...
pN
SM 1 SM 2 . . . SM N


which we can rewrite in vector form as
e = rSp + ε.

(1)

For each gene locus, the relationship between the observed
per-base exonic signal and the latent transcript proportions can
be mathematically modeled as in Equation (1), where ε is the
error term. In other words, we solve the following constrained
linear least square problem:
3

Figure 3: Transcript quantification using the observed per-base exonic expression signal. r represents the gene expression abundance parameter,
p represents the abundance proportions of the three transcript isoforms, e represents the observed per-base exonic expression signal, S represents
the splicing matrix, and Ŝ represents the modified splicing matrix. In this example, the gene has three reference transcripts and four annotated
exons. The 2nd transcript skips the 2nd exon whereas the 3rd transcript skips the 3rd exon. (A) The model without shrinkage. Each scalar
ei corresponds to the expression signal of the ith exonic position and is expected to be the sum of transcript expression covering position i.
(B) The selective and adaptive shrinkage model. A tuning parameter σ is introduced to penalize selected regions (second exon) of the selected
transcripts (the 1st and the third) without an exonic expression signal. The shrinkage level is adaptively adjusted and optimized according to
the exonic expression level over the algorithm iterations.






noted that the above formulation may be sensitive to noise because it penalizes the exonic regions with exactly zero coverage.
As such, to accommodate for noise, we penalize regions having
very low average coverage compared to the overall coverage of
the whole gene. Please refer to the Appendix for more details
about the noise handling.
We implement two algorithms to solve the above optimization problem: One-Step SASeq and Iterative SASeq. One-Step
SASeq calculates the gene expression abundance r by averaging the gene’s exonic signal. It then estimates the transcript
proportions by minimizing the error term. On the other hand,
Iterative SASeq iteratively estimates both the gene expression
abundance and the transcript proportion vector simultaneously.
While more accurate than One-Step SASeq, Iterative SASeq is
more computationally intensive.

min ke − rSpk2
p

subject to 1T p = 1,



0 ≤ p ≤ 1,

(2)

where 1T = [1, 1, 1, ...1].
{z
}
|
N
Our goal is to minimize ε while avoiding overfitting. Biologically, since only a subset of transcripts are expressed under a
single condition, it is important to avoid a possible model overfitting that uses the entire set of transcripts to explain the observed
the exonic signal. For example, in Figure 3B, it is very likely that
the reads originated from the 2nd transcript. Otherwise, there
should be exonic signal from the 2nd exon (blue) shared by the
1st and 3rd transcripts. Simply solving the constrained linear
least squares problem above will include the 3rd transcript to
explain the tower in the last region (green) as seen in Figure 3B.
We develop a selective and adaptive shrinkage approach
through modifying the splicing matrix according to the observed
exonic signal. Denoting σ as a tuning parameter, we propose a
new selective and adaptive shrinkage model with the modified
structure matrix Ŝ as the following:
(
1+σ
if Sij = 1 and ei = 0,
(3)
Ŝij =
Sij otherwise,

2.2.2

One-Step SASeq

In One-Step SASeq, we first calculate r by averaging the gene’s
exonic signal. We then rewrite the objective function in (2)
as (e − r Ŝp)T (e − r Ŝp) = pT (r 2 ŜT Ŝ)p − 2pT (r ŜT e) + eT e.
After dropping the constant term eT e, canceling out r from the
objective function and denoting A = r ŜT Ŝ and b = ŜT e, (2)
takes the following form:

1


min f (p) = pT Ap − pT b,

 p
2

for all i ∈ [1 . . . M ] and j ∈ [1 . . . N ], where σ is set to be proportional with r. The default value of σ is r. It is very likely that
any transcript with large uncovered regions in a highly expressed
gene is inactive. Thus, highly expressed genes are penalized more
for each uncovered base. On the other hand, this approach allows
us to minimally penalize genes with low expression abundance
or small uncovered regions because the modified splicing matrix
is almost equal to the original splicing matrix S. It should be






subject to 1T p = 1,

(4)

0 ≤ p ≤ 1,

where A is symmetrical and positive semidefinite. Please note
that A and b are known and only p is the subject to be optimized. The above optimization problem can be solve by a classic
convex quadratic programming. Please refer to the Appendix for
4

Figure 4: Runtime comparisons of various algorithms. (A) The runtime comparison for single-end data sets. (B) The runtime comparison for
paired-end data sets. For each plot, the horizontal axis represents the number of short reads in increasing order. The vertical axes represent
the algorithm runtimes. For both data sets, RSEM is at the magnitude of hours (the right vertical axis) while all other algorithms are at the
magnitude of minutes (the left vertical axis).

more details about the active-set method for convex quadratic
programming [41].
2.2.3

The algorithm iterates between the steps described in (5), (6)
and (7) until the proportion vector and expression abundance
values converge.

Iterative SASeq

3

In Iterative SASeq, we estimate both the gene-level expression
abundance and transcript proportions at each iteration. Denoting r (k) and p(k) as the expression abundance and transcript
proportion at the kth iteration, the expression abundance r (0)
can be initialized by averaging the per-base signal. Each iteration consists of two steps: the first step recalculates the transcript proportions while the second step updates the expression
abundance.
In the first step of the (k + 1)th iteration, the splicing matrices
are modified using the gene-level expression abundance r (k) as
follows:
(
1 + σ (k)
if Sij = 1 and ei = 0,
(k+1)
(5)
Ŝij
=
Sij otherwise,

3.1

subject to 1T p = 1,






(6)

0 ≤ p ≤ 1.

In the second step, given the new proportion vector p(k+1) from
(6), we have e = r Ŝ(k+1) p(k+1) + ε. To calculate r (k+1) , we
optimize r as follows:
r (k+1) = arg min ke − r Ŝ(k+1) p(k+1) k2 .
r

(7)

By solving (7) for scalar values of r, we obtain the following
result:
r (k+1) =

eT Ŝ(k+1) p(k+1)
.
kŜ(k+1) p(k+1) k2

Simulation Experiments

Using simulation studies, we demonstrate the accuracy
and speed of SASeq by comparing with Cufflinks [22],
RAEM [26] and RSEM [42].
We used FluxSimulator
[http://flux.sammeth.net/index.html] to simulate the whole
transcriptome sequencing experiments. FluxSimulator takes reference transcript sequences as the input, randomly generates
copy numbers for each transcript, fragments them, selects the
ones of right sizes to sequence in silico and finally outputs the
sequencing reads. Computational approaches for isoform identification and quantification can thus be compared in regards to
the simulation ground truth of transcript isoforms and their copy
numbers. FluxSimulator was used to generate 5 data sets consisting of 15 million, 30 million, 50 million, 100 million and 200
million single-end (SE) reads with lengths of 100 from around
130, 000 reference transcripts available in the Ensembl database
(version 65) with various copy numbers. FluxSimulator was also
used to generate 5 data sets with the same amount of short reads
and length for paired-end (PE) reads. Both the simulation data
used for method comparisons as well as SASeq method are freely
available at http://sammate.sourceforge.net/.
Figure 4 shows the time complexity of the different algorithms
across an increasing number of short reads on the same server (4
x Twelve-Core AMD Opteron 2.6GHz, 256GB RAM). For the
single-end data sets the runtime of single-thread SASeq is quite
comparable with the single-thread RAEM [26] and the multithread Cufflinks [22]. It is also much faster than the multi-thread
RSEM [42]. For paired-end data sets the runtimes of SASeq and
RAEM are still comparable and are much faster than RSEM and
Cufflinks. Due to the additional iterations of Iterative SASeq, it
is slightly slower than One-Step SASeq.
We proceed to compare the accuracy of transcript isoform
quantification using the true copy numbers of all of the transcript isoforms from the simulated data sets. For this purpose,
we used the Jensen-Shannon (JS) divergence since it can capture

for all i ∈ [1 . . . M ], and j ∈ [1 . . . N ], where σ (k) is adjusted
according to the expression abundance r (k) as described in (3) .
Denoting A(k+1) = r (k) (Ŝ(k+1) )T Ŝ(k+1) , b(k) = (Ŝ(k+1) )T e, we
recalculate the proportion vector p(k+1) by solving the following
convex QP:

1


p(k+1) = arg min{ pT A(k+1) p − pT b(k) },


p
2

Results

(8)

5

Figure 5: Comparison of transcript quantification accuracy using JS divergence. The upper left panel (A) compares algorithms in terms of
their divergence distributions from the ground truth for single-end data sets. The upper right panel (B) does the same for the paired-end data
sets. The lower left panels (C, D) compare algorithms in terms of their medians from the ground truth for single-end and paired-end data sets,
respectively. Both One-Step and Iterative SASeq algorithms outperform their competitors for both the single-end and paired-end data sets. The
medians of SASeq are distinguishingly lower than of the competing methods. The performance contrasts get sharper as the number of short
reads increases.

cantly smaller variances for the paired-end data sets indicate a
robust behavior of SASeq. There are likely three reasons for
the performance gaps between SASeq and other methods. First,
other approaches do not exploit the ubiquitous phenomenon of
uncovered regions. Second, the signal correction procedure helps
to yield better per-base signal and thus enhances the numeric estimation of expression abundances of active transcripts. Third,
combining available information from both read counts and perbase signal yield a much better result than using each of them
alone. This is translated into a prospectively more powerful algorithm for the forthcoming ultra-high throughput sequencing
data and longer read length. In summary, from Figure 4 and
Figure 5, SASeq demonstrates an impressive accuracy at a reasonable speed.

both linear and nonlinear relationships. Thus, the most accurate
transcript isoform quantification algorithm will yield a vector of
isoform proportions that is least divergent from the vector of
ground truth proportions. In other words, a lower value for
the JS divergence indicates a better performance. In Figure 5A
and 5B, the horizontal axes represent an increasing number of
short reads for the single-end and paired-end data sets, respectively. The vertical axes of Figure 5A and 5B represent the JS
divergence distributions for the aforementioned data sets. Lower
values indicate a better performance. In Figure 5C and 5D, we
also plot the medians of each JS divergence distribution for the
single-end and paired-end data sets, respectively. Please note
that the number of paired-end reads equals half of the number
of short reads shown in Figure 5B and 5D. For example, the
200M data set in Figure 5B and 5D actually has 100 million
paired-end reads.
In general, all the methods perform better for the paired-end
data sets except RAEM, which does not exploit paired-end information. For all data set, SASeq outperforms its competitors
with Iterative SASeq being superior to One-Step SASeq. More
strikingly, the performance gap between SASeq and other algorithms widens as the number of reads increases for both the
single-end and paired-end data sets. Furthermore, the signifi-

3.2

Real RNA-Seq Data

We also demonstrate the practical utilities of SASeq using a real
RNA-Seq data set of six Mutu I wild-type cell line and six mir155 expressing cell line samples (accession number: SRA011001).
The wild-type Mutu I cell line is of type I latency (limited viral
gene expression) B-cell line whereas the mir-155 expression cell
line is introduced by infecting the Mutu I cell line in duplicate
6

Figure 6: Predicted mir-155 targets. (A) The comparison of target prediction using a bitmap. Numbers at the top show the total predicted
number of down-regulated genes of each method. Numbers to the left show the number of genes common in each respective row. For example,
the first row indicates that 115 genes are predicted by all of the methods whereas the second row indicates that 4 genes are uniquely predicted by
SASeq. (B) The isoform-level fold change distributions of 124 targets validated by 3’-UTR assay from the previous study [43]. (C) Quantification
of the gene TAF5L. SASeq predicts both ENST00000366675 and ENST00000366676 as down-regulated in the mir-155 expressing cell line with
0.83 gene-level differential expression. The isoform ENST00000366675 has a mir-155 seed in the 3’-UTR, which was validated by both RT-PCR
and 3’-UTR assay [26]. The numeric values in green show the fold changes predicted by the competing methods and by both validations
(RT-PCR and 3’-UTR assay) for the ENST00000366675 isoform.

predicted fold changes of 124 mir-155 targets in Figure 6B, which
shows that the fold changes predicted by SASeq are more significant than other methods.
Particularly, we demonstrate in Figure 6C that SASeq identifies the down-regulated mir-155 target at the isoform-level,
which has been validated by both RT-PCR and 3’-UTR assay. In Figure 6C, the upper panel shows the exonic expression
signal for the gene TAF5L whereas the lower panel shows the
structures of five reference transcript isoforms and the associated data analysis. A fold change of 0.83 as well as a visual
inspection of the upper panel show no significant differential expression at the gene-level. Furthermore, the target transcript
isoform ENST00000366675 has been validated by both RT-PCR
and 3’-UTR assay with fold changes of 0.38 and 0.30, respectively [26]. Using SASeq we predicted the most accurate fold
change of 0.28 compared to Cufflinks’ (0.58), RAEM’s (0.56),
and RSEM’s (0.52).

with an mir-155 expressing retrovirus (or an empty vector control retrovirus) to achieve a high mir-155 expression [43]. The
goal of the data analysis is to identify the down-regulated transcripts that can be potential mir-155 targets. We picked this
RNA-Seq data set because hundreds of gene expression were
quantified and a large number of 3’-UTR assays have been done
to validate the predicted mir-155 targets [26, 43]. We first performed a per-base sequence quality check using fastQC Software
(http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/). We
then used TopHat [36] to uniquely align short reads to the human
reference genome (hg19/GRCh37). Default settings and unique
alignment were used. Alignment results were saved in a BAM
format. We re-analyzed the data set using SASeq to re-discover
the previously reported differentially expressed transcripts and
to predict novel targets.
We compared the performance of Cufflinks, RAEM, RSEM,
and SASeq in accurately identifying down-regulated transcript
isoforms. We used a “gold standard” data set of mir-155 targets
(genuine down-regulation) validated by 3’-UTR assays [43]. We
pooled the short reads from the six samples within each condition to gain more exonic coverage for the target genes. We
also removed 25% of the transcripts that were least expressed
in the control case to obtain reliable fold changes. We then
conducted isoform-level differential expression analysis to screen
down-regulated transcript isoforms of the 124 mir-155 targets
that were validated by 3’-UTR assay in the previous study [43].
Using the same cutoff value of 0.8 as in the aforementioned studies [26, 43], SASeq is able to predict more down-regulated genes
than other competing methods (Figure 6A). We also plotted the

4

Conclusions

In this paper, we presented a new approach to identify and quantify active transcripts using RNA-Seq. Our contribution is twofold. First, we combine the two types of complementary information available, i.e. the read counts and the per-base signal, to
yield a more accurate result. Second, we develop a novel shrinkage technique to shrink transcript proportions that are not supported by the observed per-base exonic expression signal and to
adaptively adjust the shrinkage level accordingly. Since the everincreasing numbers of reference transcripts in the databases as
7

well as their error rates aggravate the risk of model overfitting,
model selection via shrinkage is a promising avenue for future
research in this area.
The key innovation of our approach is the informed shrinkage, which is fundamentally different from a Lasso shrinkage
approach. First, choosing the value of the tuning parameter for
the standard Lasso shrinkage approach is not straightforward because the number of condition-specific transcripts is unknown.
Furthermore, this number varies from gene to gene and from
condition to condition. On the other hand, SASeq automatically determines the value of its tuning parameter according
to the overall exonic coverage signal of the gene. Finally, for
any gene the number of condition-specific transcripts is usually
very small compared to the number of reference transcripts available from a database. With that in mind, the Lasso approach
shrinks all the reference transcript abundance parameters in a
non-discriminative way and may not necessarily lead to a set of
condition-specific transcripts. On the other hand, SASeq only
penalizes selected regions of the transcripts that are not supported by naked exons, the remaining is more likely to be expressed under a specific biological condition.
Nevertheless, SASeq will prospectively become more powerful
in the near future with the accelerated augmentation of transcript databases and increasing sequencing depths. The former
will give a more complete set of transcripts to select from. The
latter will permit a more accurate selection and shrinkage of
the transcripts and their regions. In addition to working with
transcript databases, SASeq is able to identify and quantify active transcripts from transcriptome assembly outputs in the gtf
format [22, 44, 45, 46, 47] where model misspecification is also
an outstanding issue due to the excessive assembly bias and errors. SASeq is freely available via the GUI software SAMMate
at http://sammate.sourceforge.net/.

5
5.1

Figure 7: The average coverage matrix.

5.2

We use the active-set method [41] to solve the convex QP described above. Our convex QP can be written as:

1


min f (p) = pT Ap − pT b,

 p
2





Noise Handling

To accommodate for noise, we consider the average coverage
of each exonic regions. We first calculate the average coverage
signal over each region and then penalize those regions having
very low average coverage compared to the overall coverage of
the whole gene (Figure 7). Define a = [a1 , a2 , ..., aM ]T as a new
vector, each scalar value ai is the average coverage of the region
that the ith base belongs to. For One-Step SASeq, we change
the Equation (3) as follows:
Ŝij =

1+σ
Sij

if Sij = 1 and ai < α,
otherwise,

(9)

for all i ∈ [1 . . . M ] and j ∈ [1 . . . N ], where α is the new signalnoise cutoff parameter with a default value of r/100.
Similarly, for Iterative SASeq, we change the Equation (5) as
follows:

(k+1)
Ŝij

=

(

1 + r (k)
Sij

if Sij = 1 and ai < α(k) ,

otherwise,

subject to 1T p = 1,

(11)

0 ≤ p ≤ 1,

where A is symmetric and positive semidefinite N × N matrix,
b and p are vectors with N elements. This convex QP has one
equality constraint and 2N inequality constraints.
Algorithm
Active-Set
Method
for
Convex
QP
Set p(0) = [1, 0, . . . , 0]T as the feasible starting point;
Set W (0) to be a subset of the active constraints at p(0) ;
Set k = 0;
loop
Find ∆p(k) to minimize f (p(k) + ∆p(k) ) in the subspace
defined by W (k) ;
if (∆p(k) == 0) then
Compute Lagrange multipliers of inequality-constraints
included in W (k) ;
if (exists a negative Lagrange multiplier) then
Obtain W (k+1) by dropping the corresponding constraint from W (k) ;
Set p(k+1) = p(k) ;
else
return p(k) ;
end if
else
if (p(k) + ∆p(k) is feasible with respect to all constraints)
then
Set p(k+1) = p(k) + ∆p(k) ;
Set W (k+1) = W (k)
else
Set p(k+1) = p(k) + α(k) ∆p(k) where α(k) is chosen to
be the largest value for which all the constraints are
satisfied;
Obtain W (k+1) by adding the blocking constraints to
W (k) ;
end if
end if
Set k = k + 1;
end loop

Appendix

(

Convex Quadratic Programming Algorithm

(10)

for all i ∈ [1 . . . M ], and j ∈ [1 . . . N ], where α(k) is adjusted
according to the expression abundance r (k) as described in (9).

For a given iteration p(k) and a working set W (k) , if the ob8

jective function is not minimized in the subspace defined by
the working set, we compute the step ∆p(k) by solving an
equality-constrained QP, in which the constraints corresponding to the working set W (k) are treated as equalities. We set
p(k+1) = p(k) + ∆p(k) if it is feasible to all constraints. Otherwise, we set p(k+1) = p(k) + α(k) ∆p(k) where the step-length
α(k) is chosen to be the largest value for which all the constraints
are satisfied. The blocking constraint is also added to the working set. If the objective function is minimized (∆p(k) =0) in the
subspace defined by W (k) but one of the Lagrange multipliers
corresponding to the an inequality constraint in the working set
is negative (does not satisfy the Karush-Kuhn-Tucker condition),
we remove this constraint from the working set. Upon reaching
a Karush-Kuhn-Tucker point that minimizes the objective function over its current working set, the algorithm terminates.

[9] Wade, C. M., Giulotto, E., Sigurdsson, S., Zoli, M., Gnerre,
S., Imsland, F., Lear, T. L., Adelson, D. L., Bailey, E.,
Bellone, R. R. et al. (2009) Genome Sequence, Comparative
Analysis, and Population Genetics of the Domestic Horse.
Science, 326(5954), 865–867.
[10] Medvedev, P., Fiume, M., Dzamba, M., Smith, T. and
Brudno, M. (2010) Detecting copy number variation with
mated short reads. Genome Research, 20(11), 1613–1622.
[11] Wu, Y.-W., Rho, M., Doak, T. G. and Ye, Y. (2012) Stitching gene fragments with a network matching algorithm improves gene assembly for metagenomics. Bioinformatics,
28(18), i363–i369.
[12] Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L.
and Wold, B. (2008) Mapping and quantifying mammalian
transcriptomes by RNA-Seq. Nature Methods, 5(7), 621–
628.

Acknowledgments
The authors would like to thank Guorong Xu and Zhansheng
Duan for their initial efforts on this work. The authors also
thank Thair Judeh for his useful comments.

[13] Wilhelm, B. T., Marguerat, S., Watt, S., Schubert, F.,
Wood, V., Goodhead, I., Penkett, C. J., Rogers, J. and
Bähler, J. (2008) Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature,
453(7199), 1239–1243.

References
[1] Roberts, A., Trapnell, C., Donaghey, J., Rinn, J. and
Pachter, L. (2011) Improving RNA-Seq expression estimates
by correcting for fragment bias. Genome Biology, 12(3),
R22–14.

[14] Lister, R., O’Malley, R. C., Tonti-Filippini, J., Gregory,
B. D., Berry, C. C., Millar, A. H. and Ecker, J. R.
(2008) Highly integrated single-base resolution maps of the
epigenome in Arabidopsis. Cell, 133(3), 523–536.

[2] Wu, Z., Wang, X. and Zhang, X. (2011) Using non-uniform
read distribution models to improve isoform expression inference in RNA-Seq. Bioinformatics, 27(4), 502–508.

[15] Wang, Z., Gerstein, M. and Snyder, M. (2009) RNA-Seq: a
revolutionary tool for transcriptomics. Nature Reviews Genetics, 10(1), 57–63.

[3] Jones, D. C., Ruzzo, W. L., Peng, X. and Katze, M. G.
(2012) A new approach to bias correction in RNA-Seq.
Bioinformatics, 28(7), 921–928.

[16] Pepke, S., Wold, B. and Mortazavi, A. (2009) Computation
for ChIP-seq and RNA-seq studies. Nature Methods, 6(11s),
S22–S32.

[4] Li, W. and Jiang, T. (2012) Transcriptome assembly and
isoform expression level estimation from biased RNA-Seq
reads. Bioinformatics, 28(22), 2914–2921.

[17] Tarazona, S., Garcı́a-Alcalde, F., Dopazo, J., Ferrer, A. and
Conesa, A. (2011) Differential expression in RNA-seq: a
matter of depth. Genome research, 21(12), 2213–2223.
[18] Garber, M., Grabherr, M. G., Guttman, M. and Trapnell,
C. (2011) Computational methods for transcriptome annotation and quantification using RNA-Seq. Nature methods,
8(6), 469–477.

[5] Xing, Y., Yu, T., Wu, Y. N., Roy, M., Kim, J. and Lee,
C. (2006) An expectation-maximization algorithm for probabilistic reconstruction of full-length isoforms from splice
graphs. Nucleic Acids Research, 34(10), 3150–3160.

[19] Jiang, H. and Wong, W. H. (2009) Statistical inferences
for isoform expression in RNA-Seq. Bioinformatics, 25(8),
1026–1032.

[6] Anton, M. A., Gorostiaga, D., Guruceaga, E., Segura,
V., Carmona-Saez, P., Pascual-Montano, A., Pio, R.,
Montuenga, L. M. and Rubio, A. (2008) SPACE: an algorithm to predict and quantify alternatively spliced isoforms
using microarrays Genome Biology. Genome Biology, 9(2),
R46.

[20] Wang, X., Wu, Z. and Zhang, X. (2010) Isoform abundance
inference provides a more accurate estimation of gene expression levels in RNA-Seq. Journal of Bioinformatics and
Computational Biology, 8(Suppl 1), 177–192.

[7] Wheeler, D. A., Srinivasan, M., Egholm, M., Shen, Y.,
Chen, L., McGuire, A., He, W., Chen, Y.-J., Makhijani,
V., Roth, G. T. et al. (2008) The complete genome of an
individual by massively parallel DNA sequencing. Nature,
452(7189), 872–876.

[21] Li, B., Ruotti, V., Stewart, R. M., Thomson, J. A. and
Dewey, C. N. (2010) RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4),
493500.

[8] Medvedev, P., Stanciu, M. and Brudno, M. (2009) Computational methods for discovering structural variation with
next-generation sequencing. Nature methods, 6(11 Suppl),
S13–S20.

[22] Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A.,
Kwan, G., van Baren, M. J., Salzberg, S. L., Wold, B. J. and
Pachter, L. (2010) Transcript assembly and quantification
by RNA-Seq reveals unannotated transcripts and isoform
9

switching during cell differentiation. Nature Biotechnology,
28(5), 511–515.

[37] Ameur, A., Wetterbom, A., Feuk, L. and Gyllensten, U.
(2010) Global and unbiased detection of splice junctions
from RNA-seq data. Genome Biology, 11(3), R34+.

[23] Hu, M., Zhu, Y., Taylor, J. M. G., Liu, J. S. and Qin,
Z. S. (2012) Using Poisson mixed-effects model to quantify
transcript-level gene expression in RNA-Seq. Bioinformatics, 28(1), 63–68.

[38] Wang, K., Singh, D., Zeng, Z., Coleman, S. J., Huang, Y.,
Savich, G. L., He, X., Mieczkowski, P., Grimm, S. A., Perou,
C. M. et al. (2010) MapSplice: Accurate mapping of RNAseq reads for splice junction discovery. Nucleic Acids Research, 38(18), e178.

[24] Feng, J., Li, W. and Jiang, T. (2011) Inference of isoforms
from short sequence reads. Journal of computational biology,
18(3), 305–321.

[39] Lou, S.-K., Ni, B., Lo, L.-Y., Tsui, S. K., Chan, T.-F. and
Leung, K.-S. (2011) ABMapper: a suffix array-based tool
for multi-location searching and splice-junction mapping.
Bioinformatics, 27(3), 421–422.

[25] Nicolae, M., Mangul, S., Mandoiu, I. and Zelikovsky, A.
(2011) Estimation of alternative splicing isoform frequencies from RNA-Seq data. Algorithms for Molecular Biology,
6(1), 9.

[40] Au, K. F. F., Jiang, H., Lin, L., Xing, Y. and Wong,
W. H. H. (2010) Detection of splice junctions from pairedend RNA-seq data by SpliceMap. Nucleic acids research,
38(14), 4570–4578.

[26] Deng, N., Puetter, A., Zhang, K., Johnson, K., Zhao, Z.,
Taylor, C., Flemington, E. K. and Zhu, D. (2011) Isoformlevel microRNA-155 Target Prediction using RNAseq. Nucleic Acids Research, 39(9), e61.

[41] Nocedal, J. and Wright, S. J. (2006) Numerical Optimization, Springer Series in Operations Research and Financial
EngineeringSpringer, second edition.

[27] Katz, Y., Wang, E. T., Airoldi, E. M. and Burge, C. B.
(2010) Analysis and design of RNA sequencing experiments
for identifying isoform regulation. Nature Methods, 7(12),
1009–1015.

[42] Li, B. and Dewey, C. N. (2011) RSEM: accurate transcript
quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12(1), 323+.

[28] Bohnert, R. and Ratsch, G. (2010) rQuant.web: a tool for
RNA-Seq-based transcript quantitation. Nucleic Acids Research, 38(suppl 2), W348–W351.

[43] Xu, G., Fewell, C., Taylor, C., Deng, N., Hedges, D., Wang,
X., Zhang, K., Lacey, M., Zhang, H., Yin, Q. et al. (2010)
Transcriptome and targetome analysis in MIR155 expressing cells using RNA-seq. RNA, 16(8), 1610–1622.

[29] Nguyen, T. C., Deng, N., Xu, G., Duan, Z. and Zhu, D.
(2011) iQuant: A fast yet accurate GUI tool for transcript
quantification. In BIBM Workshops pp. 1048–1050.

[44] Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z.,
Thompson, D. A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q. et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference
genome. Nature Biotechnology, 29(7), 644–652.

[30] Li, W., Feng, J. and Jiang, T. (2011) IsoLasso: A LASSO
Regression Approach to RNA-Seq Based Transcriptome Assembly. Journal of Computational Biology, 18(11), 1693–
1707.

[45] Robertson, G., Schein, J., Chiu, R., Corbett, R., Field, M.,
Jackman, S. D., Mungall, K., Lee, S., Okada, H. M., Qian,
J. Q. et al. (2010) De novo assembly and analysis of RNAseq data. Nature Methods, 7(11), 909–912.

[31] Li, J. J., Jiang, C.-R., Brown, J. B., Huang, H. and Bickel,
P. J. (2011) Sparse linear modeling of next-generation
mRNA sequencing (RNA-Seq) data for isoform discovery and abundance estimation. Proceedings of the National
Academy of Sciences, 108(50), 19867–19872.

[46] Zhao, Z., Nguyen, T. C., Deng, N., Johnson, K. M. and
Zhu, D. (2011) SPATA: A seeding and patching algorithm
for de novo transcriptome assembly. In BIBM Workshops
pp. 26–33.

[32] Tibshirani, R. (1994) Regression Shrinkage and Selection
Via the Lasso. Journal of the Royal Statistical Society, Series B, 58, 267–288.

[47] Schulz, M. H., Zerbino, D., Vingron, M. and Birney, E.
(2012) Oases: robust de novo RNA-seq assembly across the
dynamic range of expression levels. Bioinformatics, 28(8),
1086–1092.

[33] Leng, C., Lin, Y. and Wahba, G. (2006) A note on the Lasso
and related procedures in model selection. Statistica Sinica,
16(4), 1273–1284.
[34] Zhao, P. and Yu, B. (2006) On Model Selection Consistency
of Lasso. Journal of Machine Learning Research, 7, 2541–
2563.
[35] Zou, H. (2006) The Adaptive Lasso and Its Oracle Properties. Journal of the American Statistical Association,
101(476), 1418–1429.
[36] Trapnell, C., Pachter, L. and Salzberg, S. L. (2009) TopHat:
discovering splice junctions with RNA-Seq. Bioinformatics,
25(9), 1105–1111.
10

