BACKGROUND
in higher organisms many multi-exon genes undergo alternative splicing  reactions that produce multiple splice variants, often encoding distinct, but related, protein products.
in contrast to the traditional “one gene, one protein” hypothesis, most, if not all, multi-exon genes are now believed to be subject to alternative splicing in human  <cit> , with as isoforms apparently responsible for many of the salient differences between diverse tissue types.
a significant degree of as has also been observed in various plant and animal species, although the precise magnitude and functional relevance of these events is unknown  <cit> .
while it has long been assumed that differential gene splicing plays an important role in determining the phenotypes of organisms, it has been difficult to quantify as using available high throughput methods.
however, recently developed massively parallel sequencing-by-synthesis technologies from illumina, applied biosystems and roche  <dig> life sciences have the potential to revolutionize the study of the transcriptome  <cit> .
it is now possible to produce enough high quality reads in a single run to rival traditional est libraries that have accumulated over a span of decades.
furthermore, the resulting digital counts are more comparable to the “gold standard” quantification method, quantitative pcr , and may overcome many of the shortcomings inherent in hybridization-based microarray gene expression studies.
because the technological potential of these “rna-seq” protocols is well appreciated, and rapidly advancing, methods for accurate estimation of isoform expression levels are an active area of research.
to compute an isoform expression score, the reads that map to each isoform must be converted into a quantitative expression value.
one approach is to count the number of reads that map to each transcript, normalizing against the transcript length and sequencing depth  <cit> .
unfortunately, this technique is often infeasible for as variants because many reads can map to multiple isoforms simultaneously.
recently, lacroix et al.
investigated the theoretical limitations of transcriptome reconstruction and quantification from a combinatorial perspective  <cit> .
their analysis operated under an “exact information hypothesis” whereby the exact abundances of all relevant transcribed regions is provided error-free.
however, this approach ignores the sampling process that actually generates observed data along with the associated measurement error; in practice, statistical approaches are necessary in order to obtain accurate estimates of transcript abundance.
for example, jiang and wong have described a poisson model for isoform quantification, showing how to estimate its parameters with a maximum likelihood approach  <cit> .
other authors have employed more basic  statistical approaches, for example, fisher’s exact test, to compare levels of as between treatments, e.g.
<cit> .
in this paper, we explore the use of rna-seq datasets to address the “isoform expression estimation problem” as defined in  <cit> .
it is assumed that the set of splice variants is known; the goal is to estimate the relative expression levels of these isoforms in a mixture.
obtaining precise estimates is necessary because important tissue-specific differences in as frequently involve a continuum of isoform ratios, rather than all-or-nothing expression  <cit> .
although, currently, the assumption of known isoforms may be limiting in many cases, we will soon be able to construct detailed lists of known isoforms for various organisms and tissue states using high-throughput sequencing .
a key advantage of our method over prior approaches is that our model takes into account the non-uniformity of rna-seq reads along the targeted transcripts.
in addition, our approach can be easily adapted for use with any high-throughput sequencing technology, including those that employ paired reads.
in the following sections we will describe the details of our model, demonstrate its performance on simulated and real data, and outline topics for future research.
methods
model overview
given a set of n unique as isoforms for a gene, g, it is always possible to partition rna-seq reads from g into 2n categories according to what subset of these isoforms each read is compatible with.
for example, consider two as isoforms, t <dig> and t2:
t1: aaaaaaa uuuuuuuuuu cccccccccc
t2: aaaaaaa ---------- cccccccccc
in this example, transcript t <dig> contains a cassette exon containing only “u” nucleotides.
transcript t <dig> skips this exon.
reads aligned to these transcripts can be classified into  <dig> mutually-exclusive subsets:
• subset s1: reads which contain u’s are only compatible with transcript isoform t <dig>
• subset s2: reads which contain a’s followed immediately by c’s  are only compatible with t <dig>
• subset s3: reads which contain only a’s or only c’s are compatible with both t <dig> and t <dig>
in addition, many reads, including reads containing one or more g’s, are not compatible with either t1or t <dig>  but in the following we will disregard these and only consider reads that map to at least one of the known isoforms.
let:
given the subsets introduced above, the following equation describes the probability that an individual read maps to subset si:
in general, we can assume that , the probability an individual read maps to a particular transcript, is dependent on the  frequency, Øj, of that transcript in the transcript mixture.
we will also assume that a given isoform is sampled with probability proportional to its known length.
similarly, , the probability that an individual read maps to subset si, given the read maps to transcript tj, can be worked out using the known transcript sequence and estimates of the distributions for read length and read start position .
let:
yi denote the number of reads compatible with subset si
r denote the total number of reads for the gene
pij =
=
assuming that individual reads are independent and identically distributed , we then have yi ~ binomial, and
for the example shown above, we can express this linear model in matrix form as follows:
because rp <dig> and rp <dig> will always be zero, the rank of this matrix is  <dig>  and both β <dig> and β <dig> are estimable.
although, in general, the number of rows  grows exponentially with the number of possible transcripts, it is possible to either combine or ignore uninteresting categories.
in fact, a full rank design matrix can always be constructed by considering only the n subsets consisting of reads that map to a single isoform.
distribution of read start position and read length
most methods for estimating isoform abundance assume a uniform sampling distribution for reads along the targeted transcripts .
however, it is widely acknowledged that the true distribution for read position deviates substantially from uniformity, and varies with the fragmentation protocol and sequencing technology  <cit> .
consequently, accurate methods for isoform quantification must incorporate this critical information.
we believe that these distributions should be consistent properties of the instrument and experimental protocol.
with millions of reads often available per experiment, it is feasible to determine these distributions with a high level of accuracy.
we used a kernel density approach to estimate read length and read start distributions using the observed empirical distributions observed for well-annotated transcripts .
the read length was estimated in a similar manner, resulting in an average read length of approximately  <dig> nucleotides for the illumina data set, and about  <dig> nucleotides for the  <dig> dataset.
to investigate the relationship between experimental protocol and read distribution, we also created a simple simulation that emulates the process of cdna fragmentation by nebulization.
the similarity between figure 1b and figure 1c suggests that our simulation captures the main properties of the nebulization process.
we anticipate that more detailed models, which incorporate knowledge of the physical processes of fragmentation and sequencing, should be able to describe observed distributions of read length and position even more accurately.
constructing the design matrix
let h denote the bivariate probability mass function describing the probability that a read has start position k and length m, given that this read aligns to a transcript of length l. we compute  for a particular transcript tj using the procedure detailed in figure  <dig>
estimation of β and Φ
given the construction method described above, the design matrix will always be full column rank, so  will always be fully estimable.
each yi ~ binomial.
for computational simplicity, we use the normal approximation to the binomial distribution.
for the normal linear model with a known covariance matrix, the maximum likelihood estimate  obtained using weighted least squares is the best linear unbiased estimator .
in the system described above, the variances are not known, but can be estimated from the data.
in this case, the “feasible weighted least squares” method can be used to approximate the weighted least squares solution.
in cases where a resulting  is not a valid probability, we truncate the estimate at  <dig> or  <dig>  in addition, we ensure that the total probability is one by dividing each   by the sum across all i.
under the commonly employed assumption that a given isoform is sampled with probability proportional to its length, the probability that a given read maps to a transcript isoform can be expressed as follows:
where lj is the length of transcript j in bases.
given our estimates for the βj ’s, the known lengths for each transcript, and the fact that  , the  are uniquely determined and can be obtained by computing the unique  solution  to a  set  of  linear  equations.
the sampling distribution of the resulting  cannot easily be computed analytically, but confidence intervals can be worked out empirically following a procedure based on  <cit> .
multireads
an important consideration when dealing with short length rna-seq reads is that many reads can map to more than one gene.
for example, using  <dig> nucleotide illumina reads, mortazavi, et al.
found that between 13-25% of the reads were “multireads,” matching multiple genes in the mouse genome  <cit> .
two obvious extensions to the procedure described above can be made to accommodate multireads.
one approach is to identify sets of genes and transcripts that share overlapping regions and jointly estimate the abundances of all transcripts in the set.
however, since the size of the design matrix increases rapidly with the number of transcripts jointly estimated, this method may not be ideal when the number of overlapping genes is large.
a second approach is to discard those reads that map to more than one gene, after making an appropriate adjustment to the design matrix.
when computing the probabilities as described in figure  <dig>  each potential read for a given transcript is checked against all other genes in the target dataset.
if the read does not map to any other genes, the probability  is updated as described in figure 2; however, if the read is a multiread, the read is skipped without updating .
this test can be performed efficiently after pre-computing a hash table mapping nucleotide k-mers  to genes which contain that sequence in one or more transcript isoforms.
given this modification to the probabilities in the design matrix, observed reads that map to more than one gene can now be discarded.
in the following research, we employ the first approach for the purpose of simulation and data analysis.
datasets
we tested the algorithm on two publicly available datasets.
dataset  <dig>  consists of two experiments containing approximately  <dig> million  and  <dig> million  illumina reads from floral tissue from two different arabidopsis strains .
dataset  <dig>  includes approximately  <dig> million  and  <dig> million  illumina reads from control and drought-stressed arabidopsis plant tissue.
implementation
the algorithm described above was implemented in java, with matrix computations by the jama matrix library .
data analyses and simulations were also performed using the r statistical programming language .
RESULTS
simulation
to test our method, we performed a simulated rna-seq experiment using the arabidopsis gene models defined in tair  <dig>   <cit> .
a publicly available dataset  was used to estimate the read length and read start position for illumina reads .
the simulation was performed for several of the multi-isoform genes as follows: first, a relative frequency for each of the alternative isoforms was specified as a simulation parameter, along with a predetermined total number of rna-seq reads.
each of these reads was then simulated by first selecting an isoform with probability proportional to its length and concentration in the mixture.
then, a read start position k and read length m along the selected isoform were drawn from the distribution, h. using these read coordinates, the nucleotide sequence along the sampled isoform was determined, and this sequence was compared with other isoforms in the mixture to identify which subset of isoforms the read is compatible with.
the output of one run of the simulation is a list of subsets and the corresponding counts of simulated reads assigned to those subsets.
for each gene, the entire simulation was repeated  <dig> times.
given the simulated datasets, we then used the linear model described above to infer the original isoform concentrations from the simulated subset counts.
we performed the estimation in two different ways for each gene: first, using a design matrix constructed using the same read position distribution used to generate the simulated reads  and second, using a design matrix constructed from a uniform read position distribution.
note that using the incorrect distribution can introduce a severe bias into the estimates and even change the ordering of the isoform expression levels .
for each of the  <dig> replications, we also computed approximate confidence intervals about the estimates.
to evaluate the performance of our approximation, we checked each confidence interval to see if it included the true value for the parameter.
for the at1g <dig> simulation,  <dig> % of the 95% c.i.’s contained the true parameter;  <dig> % of the 90% c.i.’s contained the true parameter; and  <dig> % of the 65% c.i.’s contained the true parameter.
this entire process was repeated for  <dig> different genes.
over several replications, the mean estimates were always approximately equal to the simulated isoform concentrations when the correct distribution was used to construct the design matrix.
for genes with two isoforms, we found that approximately 750- <dig> reads were often needed to obtain a 95% confidence interval with a width of ~20%.
however, the number of reads required varies according to mixture composition, number of isoforms, and the read length and start site distributions.
.
real rna-seq datasets
we next tested our method on two publicly available rna-seq datasets .
for each dataset, reads were mapped to the transcriptome using the soap v <dig> alignment program  <cit> .
tair  <dig> was used to define the tested gene models.
among the  <dig> arabidopsis gene models defined in tair  <dig>   <dig> genes had more than  <dig> isoform.
in particular,  <dig> genes had  <dig> isoforms,  <dig> had  <dig> isoforms,  <dig> had  <dig> isoforms,  <dig> had  <dig> isoforms,  <dig> had  <dig> isoforms,  <dig> had  <dig> isoforms,  <dig> had  <dig> isoforms, and  <dig> had  <dig> isoforms.
for dataset  <dig>  in the col- <dig> sample , among genes with  <dig> or more isoforms,  <dig> genes had more than  <dig> mapped reads ;  <dig> genes had more than  <dig> reads ;  <dig> genes had more than  <dig> reads ; and  <dig> genes had more than  <dig> reads .
similar coverage was observed in the other samples.
differential splicing
we first used a chi-square test of subset counts to identify genes that were differentially spliced between the two conditions.
genes that did not have at least one read mapping to the unique regions of at least two isoforms were not tested.
in dataset  <dig>   <dig> of  <dig> tested multi-isoform genes  showed significant differential splicing between the two treatments.
significance was determined using a p-value cut-off of approximately  <dig> , corresponding to a false discovery rate of  <dig>   <cit> .
in dataset  <dig>   <dig> of  <dig> tested multi-isoform genes  showed differential splicing, according to the same criteria.
after identifying differentially spliced genes, we then used the linear model described above to infer the isoform ratios within each of the treatment samples.
we discarded reads that mapped to more than one gene, using the adjusted design matrices computed as described in the methods section.
in several cases, genes with highly significant as levels , had only very small differences in isoform composition between the two treatments.
these most likely do not represent significant biological differences.
on the other hand, for both datasets  <dig> and  <dig>  we also identified many genes that had large differences in isoform composition between the two treatments.
in dataset  <dig>  for example,  <dig> of the  <dig> genes  had differences of 15% or more in the proportion for the main isoform; similarly, in dataset  <dig>   <dig> of  <dig> genes  showed a 15% or greater difference in the proportion for the main isoform.
figure  <dig> shows histograms illustrating the distribution of differences in main isoform proportions between the two treatments within the two datasets.
in particular, several genes appear to exhibit on/off “switch-like” differential splicing .
full listings of the differentially spliced genes identified in each dataset, along with expression level estimates for the main isoform, can be found in additional files  <dig> and  <dig>
to further characterize the genes identified as exhibiting differential splicing within the two datasets, we performed a go analysis  <cit>  using the online agrigo tool  <cit> .
the program takes as input a list of genes and compares their associated go terms with those corresponding to a relevant background set.
enriched terms are identified and grouped according to the go categories “cellular component,” “biological process,” and “molecular function.” to perform the analysis, we submitted the lists of differentially spliced genes, and selected the arabidopsis tair gene models as the background reference.
we used the hypergeometric statistical test option, with the yekutieli fdr multi-test adjustment method, set at a significance level of  <dig> .
for dataset  <dig>  which compares leaf tissue for two different untreated strains of arabidopsis, most of the significant go terms were in the cellular component category go category.
notably, many genes were associated with terms identifying cell membrane and membrane-bound cellular components, and the go terms “membrane,” “plasma membrane,” “thylakoid membrane,” and “membrane-bounded organelle,” were all significant.
this is of interest since it has been previously shown that as is often used as a mechanism to alter the functional domains of transmembrane proteins  <cit> .
in the molecular function category, the as gene list was enriched for the terms “structural molecule activity”, “structural constituent of ribosome”, “nucleic acid binding” and “ligase activity”.
in the biological process category, enriched terms include “developmental process”, “multicellular organismal development,” and “post-embryonic development”.
a full list of the significant go terms is included in additional file  <dig>  and graphical representations of the hierarchy of significant terms in the three main go categories are available in additional files  <dig>
in dataset  <dig>  which compares drought-stressed arabidopsis plants to untreated controls, the list of differentially spliced genes was significantly enriched for many terms in the biological process category related to stress response, including “response to stress”, “response to stimulus”, “response to cold,” “response to chemical stimulus,” “response to cadmium ion”, etc.
the significant biological process go terms for both datasets are shown in figure  <dig>  enriched cellular component go terms for dataset  <dig> include “membrane-bounded organelle”, “cytoplasm”, and “ribosomal subunit,” among others.
the full lists are available in additional file  <dig>  with graphical representations of the significant terms available in additional files  <dig> and  <dig>
differential gene expression
for each dataset, we also tested for differential gene expression between the two treatments .
for each gene, we computed, within each treatment, i, the proportion of reads for that gene,, relative to the total number of reads in the sample, ni.
the resulting proportions for each treatment were then compared on a gene-by-gene basis using the statistic  where.
since n <dig> and n <dig> are both very large, and the samples are independent, the test statistic, z, is assumed to follow a standard normal distribution.
because our goal was to compare differentially expressed genes with differentially spliced genes, only genes expressing more than one isoform were considered .
the resulting p-values were subsequently corrected for multiple testing using the benjamini-hochberg method, with a fdr of  <dig>   <cit> .
in dataset  <dig>   <dig> % of the tested genes showed differential expression between treatments according to these criteria; in dataset  <dig>   <dig> % of the genes showed differential expression.
in concordance with similar observations for animal systems  <cit> , differential gene expression and differential splicing events appeared to be statistically independent.
in dataset  <dig>  for example, since  <dig> % of the  <dig> multi-isoform genes were differentially expressed, one would expect a random selection from those genes to include, on average, about  <dig> % differentially expressed genes.
accordingly, the set of  <dig> differentially spliced genes included  <dig> differentially expressed genes .
in dataset  <dig>   <dig>  of the differentially spliced genes were also differentially expressed at the gene level.
.
the correlation between as and differential gene expression p-values was - <dig>  for dataset  <dig>  and  <dig>  for dataset  <dig>
multireads and connected components
in dataset  <dig>   <dig> % of the reads mapped to more than one gene; similarly,  <dig> % of the reads were multireads in dataset  <dig>  in the analysis described above, these reads were discarded according to procedure outlined in methods, after an appropriate adjustment to the design matrix.
an alternative technique for accommodating multireads is to jointly estimate all transcripts that share common read-length subsequences.
to determine what genes should be jointly estimated, and to assess the feasibility of this approach for the arabidopsis genome, we constructed a graph with tair  <dig> gene models as the vertices.
in the graph, two vertices are connected with an  edge if their corresponding gene transcripts contain a common  <dig> nucleotide substring.
the graph was then partitioned into its connected components.
a connected component is a subgraph containing a set of vertices that meet the following conditions: 1) each vertex in the set can reach any other vertex in the set along some path in the graph, and 2) each vertex can not reach any vertex outside the set.
thus, the connected components partition the graph into sets of genes that should be jointly estimated using a shared design matrix.
in order to characterize the types of genes contained in the largest components, we used the agrigo tool to analyze the associated go terms.
several of the largest components appear to correspond to gene families and many were significantly enriched for terms such as “base pairing with rna,” and “nucleic acid binding,” which are indicative of a functional significance for shared sequences.
for example, the second largest component contains a family of genes encoding trnas for the amino acid, proline.
the  <dig> members of this component share variations of the following  <dig> nucleotide sequence:
in contrast, there were no significant terms for the very largest component and many of the edges in this  <dig> gene “super-component” originated from various low complexity sequences.
table  <dig> contains descriptions of the  <dig> largest connected components, taken from the gene annotations available in tair.
discussion
we have developed a new method for using rna-seq data to quantify the alternatively spliced isoforms present in a mixture.
comparable existing methods for isoform quantification assume uniform sampling along targeted transcripts; however, real rna-seq datasets rarely exhibit uniform read sampling .
our simulations demonstrate that using the wrong read sampling distribution can lead to incorrect conclusions about the expression of isoforms.
accordingly, we have designed our method specifically to handle this problem.
some recent work has also suggested that read sampling distributions may include systematic biases that are transcript specific; for example, various aspects of rna-seq library preparation may introduce a dependence on the gc content of a given sequence and/or the terminal nucleotide ditags of the sampled sequence  <cit> .
our model, like other recent approaches to isoform quantification, does not take these sequence-specific biases into account.
however, once these biases are better understood, they can potentially be incorporated into our model, for example, in the form of sequence-specific read start and length distributions.
this is a topic that we are actively exploring.
“multireads,” or reads that map to more than one gene, are another important problem for accurate isoform quantification.
in this paper we have suggested two methods for estimation in the presence of multireads.
one approach is to jointly estimate all genes that share overlapping sequences.
an advantage of this approach, over the alternative which discards multireads  is that we make use of all available information.
however, in our connected components analysis, we have shown that the arabidopsis genome includes gene families with  <dig> or more members.
although theoretically possible with our algorithm, a joint estimation of these transcripts, using all available information, is currently infeasible due to the increasing size of the design matrix.
one major limitation of our approach is that it assumes that all transcripts are known, yet the current state of transcriptome annotation is incomplete for most organisms.
because incorrect assumptions regarding potential transcripts in a mixture could lead to erroneous estimates, we are investigating ways to incorporate residual-based diagnostics into our model.
these enhancements would serve to identify the presence of unknown “hidden” isoforms in a mixture and would complement isoform quantification with a mechanism for transcript discovery.
to demonstrate the practicality of our approach, we have applied it to two public arabidopsis rna-seq experiments and the results reveal a high level of differential splicing between strains and treatments.
in the first of the two example datasets, we identified  <dig> genes that exhibit differential splicing between two different arabidopsis strains; in the second example dataset, we found  <dig> genes that were differentially spliced between drought-stressed and untreated arabidopsis plants.
in many of these cases, the predominant isoform differed between treatments, but we also identified a large number of genes that varied more subtly in isoform expression .
in both datasets, differential splicing and gene expression appear to be statistically independent events.
CONCLUSIONS
next-generation high-throughput sequencing promises to revolutionize the study of the transcriptome.
in this paper, we have introduced a new method for using rna-seq data to quantify the alternatively spliced isoforms present in a mixture.
our method, which can incorporate a non-uniform read sampling distribution, is flexible enough to accommodate a variety of sequencing technologies, including those that incorporate paired reads.
we anticipate that the ability to reliably compute quantitative isoform expression values will help researchers to separate true alternative splicing events from spurious transcripts originating from single mis-spliced transcripts — a major problem in genome-wide alternative splicing studies.
based on recent studies  <cit>  and our own observations derived from two independent arabidopsis rna-seq experiments, isoform expression level changes frequently involve a continuum of isoform ratios, in addition to all-or-nothing expression patterns.
furthermore, in the datasets we have examined, isoform expression level changes appear to be independent of gene expression changes.
this suggests the existence of a so far uninvestigated, dynamic layer of the transcriptome.
for this reason, we believe that, in the future, researchers will prioritize genes for functional analysis based not only on observed changes in gene expression levels, but also on changes in alternative splicing.
competing interests
the authors declare that they have no competing interests.
authors' contributions
beh implemented the method and analyzed the data.
beh and sh conceived of the method and study design, and collaborated to prepare the manuscript.
all authors approved the final manuscript.
supplementary material
additional file 1
diffas_lister.xls
additional file 2
diffas_mockler.xls
additional file 3
lister_biological_process.pdf
additional file 4
lister_cellular_component.pdf
additional file 5
lister_molecular_function.pdf
additional file 6
mockler_biological_process.pdf
additional file 7
mockler_cellular_component.pdf
additional file 8
list of all connected components in a graph joining tair  <dig> genes that share 25-mer subsequences.
connectedcomponents.txt
acknowledgements
this research was supported in part by the national science foundation ios- <dig>  and by the ncsu/epa cooperative training program in environmental sciences research, training agreement ct833235-01- <dig> with north carolina state university.
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2010: selected articles from the  <dig> ieee international conference on bioinformatics and biomedicine.
the full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/11?issue=s <dig>

 in higher organisms many multi-exon genes undergo alternative splicing reactions that produce multiple splice variants, often encoding distinct, but related, protein products. in contrast to the traditional “one gene, one protein” hypothesis, most, if not all, multi-exon genes are now believed to be subject to alternative splicing in human, with as isoforms apparently responsible for many of the salient differences between diverse tissue types. a significant degree of as has also been observed in various plant and animal species, although the precise magnitude and functional relevance of these events is unknown. while it has long been assumed that differential gene splicing plays an important role in determining the phenotypes of organisms, it has been difficult to quantify as using available high throughput methods. however, recently developed massively parallel sequencing-by-synthesis technologies from illumina, applied biosystems and roche life sciences have the potential to revolutionize the study of the transcriptome. it is now possible to produce enough high quality reads in a single run to rival traditional est libraries that have accumulated over a span of decades. furthermore, the resulting digital counts are more comparable to the “gold standard” quantification method, quantitative pcr, and may overcome many of the shortcomings inherent in hybridization-based microarray gene expression studies. because the technological potential of these “rna-seq” protocols is well appreciated, and rapidly advancing, methods for accurate estimation of isoform expression levels are an active area of research. to compute an isoform expression score, the reads that map to each isoform must be converted into a quantitative expression value. one approach is to count the number of reads that map to each transcript, normalizing against the transcript length and sequencing depth. unfortunately, this technique is often infeasible for as variants because many reads can map to multiple isoforms simultaneously. recently, lacroix et al. investigated the theoretical limitations of transcriptome reconstruction and quantification from a combinatorial perspective. their analysis operated under an “exact information hypothesis” whereby the exact abundances of all relevant transcribed regions is provided error-free. however, this approach ignores the sampling process that actually generates observed data along with the associated measurement error; in practice, statistical approaches are necessary in order to obtain accurate estimates of transcript abundance. for example, jiang and wong have described a poisson model for isoform quantification, showing how to estimate its parameters with a maximum likelihood approach. other authors have employed more basic statistical approaches, for example, fisher’s exact test, to compare levels of as between treatments, e.g.. in this paper, we explore the use of rna-seq datasets to address the “isoform expression estimation problem” as defined in. it is assumed that the set of splice variants is known; the goal is to estimate the relative expression levels of these isoforms in a mixture. obtaining precise estimates is necessary because important tissue-specific differences in as frequently involve a continuum of isoform ratios, rather than all-or-nothing expression. although, currently, the assumption of known isoforms may be limiting in many cases, we will soon be able to construct detailed lists of known isoforms for various organisms and tissue states using high-throughput sequencing. a key advantage of our method over prior approaches is that our model takes into account the non-uniformity of rna-seq reads along the targeted transcripts. in addition, our approach can be easily adapted for use with any high-throughput sequencing technology, including those that employ paired reads. in the following sections we will describe the details of our model, demonstrate its performance on simulated and real data, and outline topics for future research. methods model overview given a set of n unique as isoforms for a gene, g, it is always possible to partition rna-seq reads from g into 2n categories according to what subset of these isoforms each read is compatible with. for example, consider two as isoforms, t and t2: t1: aaaaaaa uuuuuuuuuu cccccccccc t2: aaaaaaa ---------- cccccccccc in this example, transcript t contains a cassette exon containing only “u” nucleotides. transcript t skips this exon. reads aligned to these transcripts can be classified into mutually-exclusive subsets: • subset s1: reads which contain u’s are only compatible with transcript isoform t • subset s2: reads which contain a’s followed immediately by c’s are only compatible with t • subset s3: reads which contain only a’s or only c’s are compatible with both t and t in addition, many reads, including reads containing one or more g’s, are not compatible with either t1or t but in the following we will disregard these and only consider reads that map to at least one of the known isoforms. let: given the subsets introduced above, the following equation describes the probability that an individual read maps to subset si: in general, we can assume that, the probability an individual read maps to a particular transcript, is dependent on the frequency, Øj, of that transcript in the transcript mixture. we will also assume that a given isoform is sampled with probability proportional to its known length. similarly,, the probability that an individual read maps to subset si, given the read maps to transcript tj, can be worked out using the known transcript sequence and estimates of the distributions for read length and read start position. let: yi denote the number of reads compatible with subset si r denote the total number of reads for the gene pij = = assuming that individual reads are independent and identically distributed, we then have yi ~ binomial, and for the example shown above, we can express this linear model in matrix form as follows: because rp and rp will always be zero, the rank of this matrix is and both β and β are estimable. although, in general, the number of rows grows exponentially with the number of possible transcripts, it is possible to either combine or ignore uninteresting categories. in fact, a full rank design matrix can always be constructed by considering only the n subsets consisting of reads that map to a single isoform. distribution of read start position and read length most methods for estimating isoform abundance assume a uniform sampling distribution for reads along the targeted transcripts. however, it is widely acknowledged that the true distribution for read position deviates substantially from uniformity, and varies with the fragmentation protocol and sequencing technology. consequently, accurate methods for isoform quantification must incorporate this critical information. we believe that these distributions should be consistent properties of the instrument and experimental protocol. with millions of reads often available per experiment, it is feasible to determine these distributions with a high level of accuracy. we used a kernel density approach to estimate read length and read start distributions using the observed empirical distributions observed for well-annotated transcripts. the read length was estimated in a similar manner, resulting in an average read length of approximately nucleotides for the illumina data set, and about nucleotides for the dataset. to investigate the relationship between experimental protocol and read distribution, we also created a simple simulation that emulates the process of cdna fragmentation by nebulization. the similarity between figure 1b and figure 1c suggests that our simulation captures the main properties of the nebulization process. we anticipate that more detailed models, which incorporate knowledge of the physical processes of fragmentation and sequencing, should be able to describe observed distributions of read length and position even more accurately. constructing the design matrix let h denote the bivariate probability mass function describing the probability that a read has start position k and length m, given that this read aligns to a transcript of length l. we compute for a particular transcript tj using the procedure detailed in figure estimation of β and Φ given the construction method described above, the design matrix will always be full column rank, so will always be fully estimable. each yi ~ binomial. for computational simplicity, we use the normal approximation to the binomial distribution. for the normal linear model with a known covariance matrix, the maximum likelihood estimate obtained using weighted least squares is the best linear unbiased estimator. in the system described above, the variances are not known, but can be estimated from the data. in this case, the “feasible weighted least squares” method can be used to approximate the weighted least squares solution. in cases where a resulting is not a valid probability, we truncate the estimate at or in addition, we ensure that the total probability is one by dividing each by the sum across all i. under the commonly employed assumption that a given isoform is sampled with probability proportional to its length, the probability that a given read maps to a transcript isoform can be expressed as follows: where lj is the length of transcript j in bases. given our estimates for the βj ’s, the known lengths for each transcript, and the fact that, the are uniquely determined and can be obtained by computing the unique solution to a set of linear equations. the sampling distribution of the resulting cannot easily be computed analytically, but confidence intervals can be worked out empirically following a procedure based on. multireads an important consideration when dealing with short length rna-seq reads is that many reads can map to more than one gene. for example, using nucleotide illumina reads, mortazavi, et al. found that between 13-25% of the reads were “multireads,” matching multiple genes in the mouse genome. two obvious extensions to the procedure described above can be made to accommodate multireads. one approach is to identify sets of genes and transcripts that share overlapping regions and jointly estimate the abundances of all transcripts in the set. however, since the size of the design matrix increases rapidly with the number of transcripts jointly estimated, this method may not be ideal when the number of overlapping genes is large. a second approach is to discard those reads that map to more than one gene, after making an appropriate adjustment to the design matrix. when computing the probabilities as described in figure each potential read for a given transcript is checked against all other genes in the target dataset. if the read does not map to any other genes, the probability is updated as described in figure 2; however, if the read is a multiread, the read is skipped without updating. this test can be performed efficiently after pre-computing a hash table mapping nucleotide k-mers to genes which contain that sequence in one or more transcript isoforms. given this modification to the probabilities in the design matrix, observed reads that map to more than one gene can now be discarded. in the following research, we employ the first approach for the purpose of simulation and data analysis. datasets we tested the algorithm on two publicly available datasets. dataset consists of two experiments containing approximately million and million illumina reads from floral tissue from two different arabidopsis strains. dataset includes approximately million and million illumina reads from control and drought-stressed arabidopsis plant tissue. implementation the algorithm described above was implemented in java, with matrix computations by the jama matrix library. data analyses and simulations were also performed using the r statistical programming language.  simulation to test our method, we performed a simulated rna-seq experiment using the arabidopsis gene models defined in tair. a publicly available dataset was used to estimate the read length and read start position for illumina reads. the simulation was performed for several of the multi-isoform genes as follows: first, a relative frequency for each of the alternative isoforms was specified as a simulation parameter, along with a predetermined total number of rna-seq reads. each of these reads was then simulated by first selecting an isoform with probability proportional to its length and concentration in the mixture. then, a read start position k and read length m along the selected isoform were drawn from the distribution, h. using these read coordinates, the nucleotide sequence along the sampled isoform was determined, and this sequence was compared with other isoforms in the mixture to identify which subset of isoforms the read is compatible with. the output of one run of the simulation is a list of subsets and the corresponding counts of simulated reads assigned to those subsets. for each gene, the entire simulation was repeated times. given the simulated datasets, we then used the linear model described above to infer the original isoform concentrations from the simulated subset counts. we performed the estimation in two different ways for each gene: first, using a design matrix constructed using the same read position distribution used to generate the simulated reads and second, using a design matrix constructed from a uniform read position distribution. note that using the incorrect distribution can introduce a severe bias into the estimates and even change the ordering of the isoform expression levels. for each of the replications, we also computed approximate confidence intervals about the estimates. to evaluate the performance of our approximation, we checked each confidence interval to see if it included the true value for the parameter. for the at1g simulation, % of the 95% c.i.’s contained the true parameter; % of the 90% c.i.’s contained the true parameter; and % of the 65% c.i.’s contained the true parameter. this entire process was repeated for different genes. over several replications, the mean estimates were always approximately equal to the simulated isoform concentrations when the correct distribution was used to construct the design matrix. for genes with two isoforms, we found that approximately 750- reads were often needed to obtain a 95% confidence interval with a width of ~20%. however, the number of reads required varies according to mixture composition, number of isoforms, and the read length and start site distributions.. real rna-seq datasets we next tested our method on two publicly available rna-seq datasets. for each dataset, reads were mapped to the transcriptome using the soap v alignment program. tair was used to define the tested gene models. among the arabidopsis gene models defined in tair genes had more than isoform. in particular, genes had isoforms, had isoforms, had isoforms, had isoforms, had isoforms, had isoforms, had isoforms, and had isoforms. for dataset in the col- sample, among genes with or more isoforms, genes had more than mapped reads; genes had more than reads; genes had more than reads; and genes had more than reads. similar coverage was observed in the other samples. differential splicing we first used a chi-square test of subset counts to identify genes that were differentially spliced between the two conditions. genes that did not have at least one read mapping to the unique regions of at least two isoforms were not tested. in dataset of tested multi-isoform genes showed significant differential splicing between the two treatments. significance was determined using a p-value cut-off of approximately, corresponding to a false discovery rate of. in dataset of tested multi-isoform genes showed differential splicing, according to the same criteria. after identifying differentially spliced genes, we then used the linear model described above to infer the isoform ratios within each of the treatment samples. we discarded reads that mapped to more than one gene, using the adjusted design matrices computed as described in the methods section. in several cases, genes with highly significant as levels, had only very small differences in isoform composition between the two treatments. these most likely do not represent significant biological differences. on the other hand, for both datasets and we also identified many genes that had large differences in isoform composition between the two treatments. in dataset for example, of the genes had differences of 15% or more in the proportion for the main isoform; similarly, in dataset of genes showed a 15% or greater difference in the proportion for the main isoform. figure shows histograms illustrating the distribution of differences in main isoform proportions between the two treatments within the two datasets. in particular, several genes appear to exhibit on/off “switch-like” differential splicing. full listings of the differentially spliced genes identified in each dataset, along with expression level estimates for the main isoform, can be found in additional files and to further characterize the genes identified as exhibiting differential splicing within the two datasets, we performed a go analysis using the online agrigo tool. the program takes as input a list of genes and compares their associated go terms with those corresponding to a relevant background set. enriched terms are identified and grouped according to the go categories “cellular component,” “biological process,” and “molecular function.” to perform the analysis, we submitted the lists of differentially spliced genes, and selected the arabidopsis tair gene models as the background reference. we used the hypergeometric statistical test option, with the yekutieli fdr multi-test adjustment method, set at a significance level of. for dataset which compares leaf tissue for two different untreated strains of arabidopsis, most of the significant go terms were in the cellular component category go category. notably, many genes were associated with terms identifying cell membrane and membrane-bound cellular components, and the go terms “membrane,” “plasma membrane,” “thylakoid membrane,” and “membrane-bounded organelle,” were all significant. this is of interest since it has been previously shown that as is often used as a mechanism to alter the functional domains of transmembrane proteins. in the molecular function category, the as gene list was enriched for the terms “structural molecule activity”, “structural constituent of ribosome”, “nucleic acid binding” and “ligase activity”. in the biological process category, enriched terms include “developmental process”, “multicellular organismal development,” and “post-embryonic development”. a full list of the significant go terms is included in additional file and graphical representations of the hierarchy of significant terms in the three main go categories are available in additional files in dataset which compares drought-stressed arabidopsis plants to untreated controls, the list of differentially spliced genes was significantly enriched for many terms in the biological process category related to stress response, including “response to stress”, “response to stimulus”, “response to cold,” “response to chemical stimulus,” “response to cadmium ion”, etc. the significant biological process go terms for both datasets are shown in figure enriched cellular component go terms for dataset include “membrane-bounded organelle”, “cytoplasm”, and “ribosomal subunit,” among others. the full lists are available in additional file with graphical representations of the significant terms available in additional files and differential gene expression for each dataset, we also tested for differential gene expression between the two treatments. for each gene, we computed, within each treatment, i, the proportion of reads for that gene,, relative to the total number of reads in the sample, ni. the resulting proportions for each treatment were then compared on a gene-by-gene basis using the statistic where. since n and n are both very large, and the samples are independent, the test statistic, z, is assumed to follow a standard normal distribution. because our goal was to compare differentially expressed genes with differentially spliced genes, only genes expressing more than one isoform were considered. the resulting p-values were subsequently corrected for multiple testing using the benjamini-hochberg method, with a fdr of. in dataset % of the tested genes showed differential expression between treatments according to these criteria; in dataset % of the genes showed differential expression. in concordance with similar observations for animal systems, differential gene expression and differential splicing events appeared to be statistically independent. in dataset for example, since % of the multi-isoform genes were differentially expressed, one would expect a random selection from those genes to include, on average, about % differentially expressed genes. accordingly, the set of differentially spliced genes included differentially expressed genes. in dataset of the differentially spliced genes were also differentially expressed at the gene level.. the correlation between as and differential gene expression p-values was - for dataset and for dataset multireads and connected components in dataset % of the reads mapped to more than one gene; similarly, % of the reads were multireads in dataset in the analysis described above, these reads were discarded according to procedure outlined in methods, after an appropriate adjustment to the design matrix. an alternative technique for accommodating multireads is to jointly estimate all transcripts that share common read-length subsequences. to determine what genes should be jointly estimated, and to assess the feasibility of this approach for the arabidopsis genome, we constructed a graph with tair gene models as the vertices. in the graph, two vertices are connected with an edge if their corresponding gene transcripts contain a common nucleotide substring. the graph was then partitioned into its connected components. a connected component is a subgraph containing a set of vertices that meet the following conditions: 1) each vertex in the set can reach any other vertex in the set along some path in the graph, and 2) each vertex can not reach any vertex outside the set. thus, the connected components partition the graph into sets of genes that should be jointly estimated using a shared design matrix. in order to characterize the types of genes contained in the largest components, we used the agrigo tool to analyze the associated go terms. several of the largest components appear to correspond to gene families and many were significantly enriched for terms such as “base pairing with rna,” and “nucleic acid binding,” which are indicative of a functional significance for shared sequences. for example, the second largest component contains a family of genes encoding trnas for the amino acid, proline. the members of this component share variations of the following nucleotide sequence: in contrast, there were no significant terms for the very largest component and many of the edges in this gene “super-component” originated from various low complexity sequences. table contains descriptions of the largest connected components, taken from the gene annotations available in tair. discussion we have developed a new method for using rna-seq data to quantify the alternatively spliced isoforms present in a mixture. comparable existing methods for isoform quantification assume uniform sampling along targeted transcripts; however, real rna-seq datasets rarely exhibit uniform read sampling. our simulations demonstrate that using the wrong read sampling distribution can lead to incorrect conclusions about the expression of isoforms. accordingly, we have designed our method specifically to handle this problem. some recent work has also suggested that read sampling distributions may include systematic biases that are transcript specific; for example, various aspects of rna-seq library preparation may introduce a dependence on the gc content of a given sequence and/or the terminal nucleotide ditags of the sampled sequence. our model, like other recent approaches to isoform quantification, does not take these sequence-specific biases into account. however, once these biases are better understood, they can potentially be incorporated into our model, for example, in the form of sequence-specific read start and length distributions. this is a topic that we are actively exploring. “multireads,” or reads that map to more than one gene, are another important problem for accurate isoform quantification. in this paper we have suggested two methods for estimation in the presence of multireads. one approach is to jointly estimate all genes that share overlapping sequences. an advantage of this approach, over the alternative which discards multireads is that we make use of all available information. however, in our connected components analysis, we have shown that the arabidopsis genome includes gene families with or more members. although theoretically possible with our algorithm, a joint estimation of these transcripts, using all available information, is currently infeasible due to the increasing size of the design matrix. one major limitation of our approach is that it assumes that all transcripts are known, yet the current state of transcriptome annotation is incomplete for most organisms. because incorrect assumptions regarding potential transcripts in a mixture could lead to erroneous estimates, we are investigating ways to incorporate residual-based diagnostics into our model. these enhancements would serve to identify the presence of unknown “hidden” isoforms in a mixture and would complement isoform quantification with a mechanism for transcript discovery. to demonstrate the practicality of our approach, we have applied it to two public arabidopsis rna-seq experiments and the results reveal a high level of differential splicing between strains and treatments. in the first of the two example datasets, we identified genes that exhibit differential splicing between two different arabidopsis strains; in the second example dataset, we found genes that were differentially spliced between drought-stressed and untreated arabidopsis plants. in many of these cases, the predominant isoform differed between treatments, but we also identified a large number of genes that varied more subtly in isoform expression. in both datasets, differential splicing and gene expression appear to be statistically independent events.  next-generation high-throughput sequencing promises to revolutionize the study of the transcriptome. in this paper, we have introduced a new method for using rna-seq data to quantify the alternatively spliced isoforms present in a mixture. our method, which can incorporate a non-uniform read sampling distribution, is flexible enough to accommodate a variety of sequencing technologies, including those that incorporate paired reads. we anticipate that the ability to reliably compute quantitative isoform expression values will help researchers to separate true alternative splicing events from spurious transcripts originating from single mis-spliced transcripts — a major problem in genome-wide alternative splicing studies. based on recent studies and our own observations derived from two independent arabidopsis rna-seq experiments, isoform expression level changes frequently involve a continuum of isoform ratios, in addition to all-or-nothing expression patterns. furthermore, in the datasets we have examined, isoform expression level changes appear to be independent of gene expression changes. this suggests the existence of a so far uninvestigated, dynamic layer of the transcriptome. for this reason, we believe that, in the future, researchers will prioritize genes for functional analysis based not only on observed changes in gene expression levels, but also on changes in alternative splicing. competing interests the authors declare that they have no competing interests. authors' contributions beh implemented the method and analyzed the data. beh and sh conceived of the method and study design, and collaborated to prepare the manuscript. all authors approved the final manuscript. supplementary material additional file 1 diffas_lister.xls additional file 2 diffas_mockler.xls additional file 3 lister_biological_process.pdf additional file 4 lister_cellular_component.pdf additional file 5 lister_molecular_function.pdf additional file 6 mockler_biological_process.pdf additional file 7 mockler_cellular_component.pdf additional file 8 list of all connected components in a graph joining tair genes that share 25-mer subsequences. connectedcomponents.txt acknowledgements this research was supported in part by the national science foundation ios- and by the ncsu/epa cooperative training program in environmental sciences research, training agreement ct833235-01- with north carolina state university. this article has been published as part of bmc bioinformatics volume supplement 2010: selected articles from the ieee international conference on bioinformatics and biomedicine. the full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/11?issue=s 
in two publicly available rna-seq datasets, we identified several alternatively-spliced genes with switch-like, on/off expression properties, as well as a number of other genes that varied more subtly in isoform expression. given that changes in isoform expression level frequently involve a continuum of isoform ratios, rather than all-or-nothing expression, and that they are often independent of general gene expression changes, we anticipate that our research will contribute to revealing a so far uninvestigated layer of the transcriptome. 1– november ieee international conference on bioinformatics and biomedicine washington, dc, usa in eukaryotes, alternative splicing often generates multiple splice variants from a single gene. given a set of known splice variants, the goal is to estimate the relative abundance of the individual variants. we believe that, in the future, researchers will prioritize genes for functional analysis based not only on observed changes in gene expression levels, but also on changes in alternative splicing. in many cases, genes exhibiting differential expression of alternatively spliced transcripts were not differentially expressed at the gene level. preliminary tests indicate that the model performs well on both simulated and real data.
 in eukaryotes, alternative splicing often generates multiple splice variants from a single gene. here weexplore the use of rna sequencing datasets to address the isoform quantification problem. given a set of known splice variants, the goal is to estimate the relative abundance of the individual variants. methods our method employs a linear models framework to estimate the ratios of known isoforms in a sample. a key feature of our method is that it takes into account the non-uniformity of rna-seq read positions along the targeted transcripts.  preliminary tests indicate that the model performs well on both simulated and real data. in two publicly available rna-seq datasets, we identified several alternatively-spliced genes with switch-like, on/off expression properties, as well as a number of other genes that varied more subtly in isoform expression. in many cases, genes exhibiting differential expression of alternatively spliced transcripts were not differentially expressed at the gene level.  given that changes in isoform expression level frequently involve a continuum of isoform ratios, rather than all-or-nothing expression, and that they are often independent of general gene expression changes, we anticipate that our research will contribute to revealing a so far uninvestigated layer of the transcriptome. we believe that, in the future, researchers will prioritize genes for functional analysis based not only on observed changes in gene expression levels, but also on changes in alternative splicing. 1– november ieee international conference on bioinformatics and biomedicine washington, dc, usa 