BACKGROUND
more than 90% of human genes  <cit>  are estimated to be alternatively spliced which leads a single gene to produce multiple proteins with distinct functions and is implicated in many diseases including cancer  <cit> .
in recent years, there is an increasing interest in the use of alternative splicing in developing diagnostic tools and in identifying new therapeutic targets  <cit> .
microarrays have been widely used to analyze alternative isoforms by combining exon arrays and exon junction arrays to quantify isoform level expression indexes  <cit> .
however, array based techniques are encountering several fundamental problems such as cross hybridization and weak signals in junction probes which are difficult to overcome  <cit> .
ultra high-throughput sequencing of rna has been developed as an approach for transcriptome analysis in several different species and has offered an attractive approach to measure transcription in a comprehensive manner.
rna-seq allows the direct detection of alternative splicing using the reads mapped at splice junctions including the novel splicing without the annotation information.
genome-wide measurements of transcriptomes are increasingly done by rna-seq which provides a far more precise measurement of expression levels of isoforms than other methods  <cit> .
the rapidly-developing rna-seq techniques require substantial algorithmic advances.
several tools and strategies have been proposed to deal with the complex bioinformatics analysis of rna-seq  <cit> .
pepke et al.
<cit>  provided a comprehensive and up-to-data review of multi-layered analyses of rna-seq data.
mortazavi et al.
<cit>  proposed to quantify the gene level expression of a transcript as reads per kilobase per million mapped reads .
further, jiang and wong  <cit>  presented a statistical model to describe how the isoform expression levels were calculated from the number of reads mapped to the annotated exons of a gene.
meanwhile, bohnert et al.
<cit>  also proposed rquant to determine the abundances for each annotated isoforms by minimizing the deviation of the observation from the expected position-wide read coverage.
all these methods assumed that the number and structures of isoforms of each gene are known from the reference genome.
however, as jiang and wong  <cit>  pointed out, the isoform level annotation is very incomplete due to the complexity of the transcriptome and the limitations of previous experimental approaches.
to address this issue, trapnell et al.
proposed cufflinks  <cit>  to identify transcripts as well as to estimate the expression levels of identified transcripts from mapped reads without annotation information.
in essence, cufflinks constructs a covering relation on the read alignments from tophat  <cit> , and find a minimum path cover on the directed acyclic graph for the relation based on dilworth's theorem  <cit>  to construct a parsimonious set of transcripts.
after that, the expression levels of those constructed transcripts are estimated using established known isoform expression estimation methods  <cit> .
therefore, the construction of transcripts in cufflinks is independent of the expression level estimation.
however, the construction of expressed transcripts and expression level estimation are highly associated.
we argue that the determination of parsimonious set of expressed transcripts and expression level estimation should be implemented jointly.
though scripture  <cit>  can also detect the novel isoforms, the issue of parsimonious expressed isoforms is not addressed.
we also notice that feng et al.
proposed isoinfer  <cit>  to identify isoforms using the detected junctions.
the candidate isoforms were constructed by combining the putative exons followed by selecting a minimum best subset from all the enumerated subsets of the candidate isoforms which can explain the observation best.
however, enumerating all possible subsets of the candidate isoform set with a given size is often infeasible.
isoinfer decomposes the large putative exon set into subsets to address this issue which introduces more parameters.
here, we put forward nsmap to infer the structures of isoforms as well as to estimate the expression levels simultaneously.
first, the exons are constructed based on the detected splicing junctions from rna-seq data using tophat.
all the possible isoforms are enumerated by combination of those detected exons.
then nsmap is applied to identify the true expressed isoforms from the large candidate isoform set as well as to estimate the expression levels with a sparsity control term to restrict the number of expressed isoforms.
the assumption behind this sparsity term is that only as few isoforms as possible should be selected to best explain the observed number of reads falling on each exon of a gene.
finally, a model selection step is conducted to select the solution which compromises the fitting of the observation and the number of expressed isoforms best.
in summary, our algorithm allows for discovering the structures of the expressed isoforms of a given gene and for estimating the concentration of each spliced isoform simultaneously without the annotated isoform information, which makes the identification of new, previously unknown, alternatively spliced isoforms possible.
this study will help rna-seq, a next generation sequencing technology, advance to its full potential in comprehensive transcriptome analysis.
RESULTS
data set
we test nsmap on two simulations with simulated expression levels derived from two publicly available mouse rna-seq as described in  <cit> .
we also apply nsmap on three real in-house rna-seq data sets of myelodysplastic syndromes  transcriptome analysis to identify isoforms including novel ones featured in mds disease.
mds are a diverse collection of hematological conditions united by ineffective production of myeloid blood cells and risk of transformation to acute myelogenous leukemia whose frequency and incidence are increasing in the us population  <cit> .
in our application, cryopreserved marrow cells and paraffin embedded marrow clot sections and marrow core biopsies from  <dig> mds patients as well as  <dig> age-matched control sample are being studied by dr. jeff chang's lab at the methodist hospital.
these mds patients have been thoroughly evaluated for clinical/morphologic/immunophenotypic data and characterized clinically by transfusion dependency and pathologically by significant dysplasia, increased blasts, and immunophenotypic aberrancy.
the control sample are obtained from patients without cytopenias .
we specifically selected these controls to be age-matched for the mdss population to control for the possibility of aging-related changes in the expression profile mrna of hematopoietic cells.
then we apply the rna-seq protocol to sequence our samples.
we sequenced the two mds clinical samples and one normal sample using illumina genome analyzer ii.
there are around  <dig> million single-end reads with read length  <dig> bp for each sample.
algorithm summary
nsmap comprise four consecutive steps, starting with junction detection and reads mapping using tophat  <cit>  and followed by the candidate isoforms construction, expressed isoform identification and expression level estimation along whole regularization path and model selection to select the best solution from the whole solution path.
short reads alignment is the first step in understanding next-generation sequencing data and many free alignment software packages are available  <cit> .
here we use tophat to perform the alignment task which can detect the junctions and map a massive amount of reads to the whole genome flexibly and efficiently.
the reference genome sequences are downloaded from ucsc genome database  <cit> .
after read alignment, the next step is to generate the candidate isoforms according to the alignment and splice junctions obtained from tophat.
candidate isoforms construction
based on the alignment results and detected junctions from tophat, the exons can be constructed from segments whose two ends have been detected as junction points.
for example, figure  <dig> mimics a gene with nine exons constructed from the detected junctions in a gene region.
the isoforms are formed from combination of those exons.
each rectangle represents an exon and arrow line indicates that there is a splice junction detected between the two exons.
this can be interpreted as a directed graph with each exon as nodes.
each candidate isoform can be generated by finding a path between two nodes on this graph.
for instance, figure  <dig> displays the possible isoforms with exon  <dig> and  <dig> as source and sink respectively.
there are total four paths connecting exon  <dig> and  <dig> as indicated in figure  <dig>
isoform expression level estimation
given the candidate isoforms constructed from last step, we need to select the expressed isoforms from this candidate set as well as to estimate the express level of selected isoforms.
we adapted the expression level estimation framework of jiang and wong  <cit>  by incorporating a sparsity regularization term.
for a gene g, suppose it has m exons with lengths  and n isoforms with expression Φ = .
assuming each exon can be either included in an isoform or not, we have a set of observations , where  is an index set of events which we are interested in.
each observation xs ∈ x is a random variable representing the number of reads falling into a certain region of interest in gene g. for example, reads falling into certain exon or exon-exon junction.
the natural statistical model of count data is the poisson distribution.
each xs ∈ x follows poisson distribution with a parameter λ which is the expected count or the mean parameter of the poisson distribution.
for instance, the λ for the number of reads falling into exon j is , where n is the total number of mapped reads and cij is  <dig> if isoform i contained exon j and  <dig> otherwise.
for a exon-exon junction event, the λ is , where l is the length of the junction region, j and k are indexes of the two exons involved in the junction being investigated.
in general, λs is a linear function of ϕ <dig> ϕ <dig> ..,ϕn, i.e.
where asi is a known coefficient.
the likelihood function is
maximum likelihood estimation  method can be used to estimate Φ by maximizing above likelihood function  <cit> .
however, in the work of jiang and wong  <cit> , the structure of individual isoform is fetched from the reference genome where only a few of isoforms are known for each gene.
in our application, the structures of expressed isoforms are unknown which need to be selected from the large candidate isoform set where each isoform is formed by combination of detected exons.
to do this, the most important prior information is arguably sparsity; that is, characteristic pattern of each gene in cells are expressed by very few isoforms, even though the total number of possible combinations of exons in a gene can be large.
utilization of this sparsity information is crucial to the success of finding the most relevant isoforms which explain the observed reads best, in the face of "insufficient" data and/or data uncertainty  <cit> .
to incorporate this sparsity to our model, we propose a prior laplace distribution on Φ as below:
laplace distribution assigns higher probability at zero than gaussian distribution to produce sparser solution and has been widely used to encode sparseness prior  <cit> .
hence we can employ maximum a posteriori  instead of maximum likelihood to estimate Φ:
taking logarithm on the above equation, we get
keep in mind that the expression of each isoform cannot be negative, that is ϕi ≥  <dig> , f can be rewritten as:
where ϕi ≥  <dig> for i =  <dig> ...,s and s is the number of events with which we are interested.
maximizing f with respect to Φ is equivalent to minimizing
to that end, the optimization problem is summarized as:
we call this nonnegativity and sparsity constrained map model as nsmap.
laplace prior distribution has a l <dig> norm term to impose the sparsity.
in some cases of our experiment, we find it may not work well in the identification of the correct expressed isoforms due to the similarities between the isoforms.
the ideal sparsity constraint is the l <dig> regularization, defined as the number of non-zero entries in a vector.
unfortunately, the computation of l <dig> regularization is intractable to solve.
l <dig> norm is popular because of its intrinsic convex property.
some recent works suggest that non-convex regularizations such as lp regularization with  <dig> < p <  <dig> have better performances in parameter selection and sparse signal construction  <cit>  than l <dig> regularization.
this inspires us to use a stronger sparsity constraint with lp,  <dig> <p <  <dig> norm to solve our problem.
we will implement a simulation study to compare the performances of sparsity regularizations with different p. in our experiment, we set p =  <dig> .
so the objective function is finally given by:
where  and .
then the expression levels of all the candidate isoforms can be obtained by optimizing equation  with a given t. when t → +∞, the solution  equals zero without expressed isoforms.
with the decreasing of t, some elements of  will be non-zero to become expressed isoform.
so we have to select an optimal t. before doing this, we will calculate the solution path which consists of solutions corresponding to different values of t. blasso  <cit>  is adapted to approximate this solution path efficiently .
model selection from solution path
after getting the solution path which consists of solutions corresponding to different values of t, we have to select the best solution  from this solution path which makes a good balance between the number of expressed isoforms and the fitting function  in equation .
the solution path will first be grouped into subsets with increasing model size according to the number of expressed isoforms of each solution .
here the model size means the number of expressed isoforms.
then the solution with minimal  in each group is selected as the best solution for each model size.
starting from model size  <dig>  we compare the  of current model size with the next model size which have one more expressed isoforms.
if  is significantly improved by the larger model size, the model with larger model size will be updated to the current model and compared with the remaining models with larger model size than the updated model.
otherwise the solution  corresponding to the current model size is selected as the final solution.
in this way, the solution with smaller number of expressed isoforms will be selected preferably.
see the methods for more detail.
simulation on the whole genome with different lp norm in equation
in the absence of rna-seq data from samples for which we have ground truth isoform quantities, we conduct simulations to validate our method and evaluate its performance in terms of isoform identification and expression levels estimation with p =  <dig> and p =  <dig> .
we first derive expression level of each isoform from the mouse brain and liver rna-seq data sets described in  <cit>  using the poisson model of  <cit>  with the mouse ucsc genes as the annotated reference genome database.
those derived expression levels are employed as the ground truth to perform the following simulations.
which exon or splice region the read will fall on is determined by uniformly sampling proportional to the simulated ground truth rpkm and the same mapped reads number of each gene in the real data.
after uniformly sampling and counting reads number, we can identify isoform structures and their concentrations from the simulated rna-seq reads data and evaluate the performance of nsmap with the simulated ground truth.
we only consider isoforms with expression levels no less than  <dig> rpkm as expressed isoforms.
in this simulation, constructed isoforms containing more than half of the total exons are eligible as the candidate isoforms.
the accuracy of nsmap largely depends on the critical step of candidate isoform generation.
so, we check whether our candidate isoforms generated from the exon reads and exon junctions can cover the true expressed isoforms.
from table  <dig> we see that more than 96% of expressed isoforms in the two data sets are included in the candidate isoforms generated by our strategy which finds paths on a graph constructed from exons and splice junctions.
isoforms included in
we further compute the fraction of isoforms for which the estimates are significantly consistent with the simulated ground truth .
we refer to this statistic as the positive fraction  <cit> .
given the positive isoforms which are within the 5% deviation from the ground truths, overall specificities of the two simulations are also calculated by setting the truly expressed isoforms with rpkm no less than  <dig> as positive isoforms and the false reported isoforms by nsmap and truly expressed isoforms with rpkm <  <dig> as negative isoforms.
the results are summarized in table  <dig>  we can see that nsmap achieved near 80% overall positive fraction in both simulations.
and highly expressed isoforms  have high positive fraction.
while for isoforms with low rpkm, the accuracy is not very good because the number of reads falling on this isoform is small which may not be sufficient to capture the exons and splice junctions to form this isoform.
this phenomenon is consistent with other studies  <cit>  that the estimation on the low abundance isoform is not very accurate.
comparing the performances of nsmap with p =  <dig>  and p =  <dig>  we notice that lp norm with p =  <dig>  provides better results than p =  <dig> in the simulations.
to explain this observation, we conduct the following experiment and give a mathematical interpretation of the two lp norms.
demonstrations of features of lp norm with p =  <dig>  and p = 1
first, we exemplify our observation in the above simulation that lp norm with p =  <dig>  is better than p =  <dig> through a simulation using the structure of gene eif5a.
in the ucsc genome annotation, gene eif5a has  <dig> exons and  <dig> known isoforms.
we set both rpkm expression levels of the fifth and sixth isoforms as  <dig> because the two isoforms express most from the estimation using the method of  <cit>  on the mouse brain and liver rna-seq data in  <cit> .
the simulation is performed based on the expression levels and structures of the two expressed isoforms under the uniform distribution for reads.
based on the simulation data, total  <dig> candidate isoforms is constructed.
the indexes of the two truly expressed isoforms in the candidate set are  <dig> and  <dig>  respectively.
from the upper row of figure  <dig>  we can see that the two expressed isoforms are selected correctly when p = 1/ <dig>  while in the case of p =  <dig>  additional two false unexpressed isoforms  are also selected.
then, the lower row of figure  <dig> shows the mathematical properties of the two norms.
l1/ <dig> norm regularization has more similar sparsity property with l <dig> norm than l <dig> norm because the points with high probabilities in the laplace-like distribution with l1/ <dig> norm are more focusing around the axes than those in the laplace distribution.
this means l1/ <dig> norm regularization imposed stronger sparsity than l <dig> norm.
the simulation depicts the stronger sparsity constraint by setting p = 1/ <dig> is superior to laplace priori distribution with p =  <dig> in our application.
in the following experiments, we select the lp norm with p = 1/ <dig>
comparison with isoinfer on two rna-seq data sets of mds samples
here we compare nsmap with isoinfer which have similar ideas in isoform construction from putative exons and minimum expressed isoform set selection.
however, isoinfer selects a subset as expressed isoforms from candidate isoform set by enumerating all the possible subsets of the candidate isoform set.
in nsmap, the selection of expressed isoforms is embedded into the isoform expression level estimation framework by incorporating a sparsity control term.
a transcript can be constructed from all the exon-intron boundaries as well as the transcription start site  and poly-a site  of an isoform.
the exon-intron boundaries can be inferred from rna-seq using alternative splicing detection tool, such as tophat and splicemap  <cit> .
the tss and pas represent the start and end expressed segments of a transcript, respectively.
theoretically, any expressed segments can be the tss or pas which will introduce many false short isoforms to make isoform inference difficult.
we prefer to retrieve the tss-pas from the ucsc known isoform table as the starts and ends of predicted transcripts and to identify isoforms within the regions of known genes whose functions and pathways are intensively studied.
isoinfer  and nsmap will use those tss-pas and the detected junctions using tophat to infer the expressed isoforms from rna-seq.
because it is infeasible to validate all the predicted isoforms, we evaluate the two methods by comparing their predictions with ucsc known isoform data set.
the performance of the method is measured by sensitivity and precision.
here we use hg <dig> known human isoforms data set downloaded from ucsc which contains  <dig>  transcripts.
a known isoform is identified if it is in the prediction result of a method.
sensitivity is defined as the number of identified isoforms divided by the number of all known isoforms from ucsc data base.
precision is defined as the number of identified known isoforms divided by the number of predicted isoforms by the method.
however, isoinfer selects a subset as expressed isoforms from candidate isoform set by evaluating all the possible subsets of the candidate isoforms.
in nsmap, the selection of expressed isoforms is embedded into the isoform expression level estimation framework by incorporating a sparsity control term.
in this way, the selection of expressed isoforms is automatic and more efficient than testing all the possible sub sets of the candidate isoforms.
we also notice that more known isoforms are predicted in mds sample  <dig> than mds sample  <dig>  because there are more reads are mapped in mds sample  <dig>  this observation tells us that deeper sequencing will improve the performances of isoinfer and nsmap.
here, the sensitivities of the two methods are very low.
the reason is that we compare the predicted isoforms with the large ucsc known isoform set.
some of the ucsc known isoforms may not express in the sample.
so the effective sensitivities will be larger than those numbers.
this comparison against ucsc known isoform data set does not mean all the predicted isoforms without annotation are false.
especially those predicted novel isoforms with high rpkm are promising to be true novel isoforms.
for example, if we select the predicted isoforms with expression level larger than  <dig> rpkm, in mds sample  <dig> and  <dig>   <dig> out of  <dig> and  <dig> out of  <dig> predicted isoforms by nsmap will be annotated in the ucsc known isoforms table.
and the lowly expressed isoforms have higher probability to be false positive or artifacts due to the insufficient reads for capturing the true structure of an isoform.
so we need to set an expression level threshold to refine the predicted isoforms.
this issue is addressed in the following section.
example of identified isoforms using nsmap and expression level threshold selection
here we exemplify isoforms of a gene tcf <dig> estimated by nsmap.
tcf <dig> encodes protein transcription factor  <dig> which localizes to the nucleus and displays dna-binding and transactivation activities.
tcf <dig> has five exons and alternative splicing results in two known transcript variants encoding different isoforms in the ucsc reference genome ).
the rna-seq data is first processed by tophat to detect the splicing junctions.
figure  <dig> shows the splicing junctions of tcf <dig> on mds sample  <dig> detected by tophat.
then nine exons are constructed based on detected junctions in figure  <dig>  nsmap predicts three isoforms with expression level  <dig> ,  <dig>  and  <dig>  respectively ).
compared with the annotated isoforms of tcf <dig> ), we confirm that the top two highly expressed predictions are the same with the annotations.
this example demonstrates the ability of nsmap in isoform identification.
however, nsmap still predicts an additional lowly expressed isoform  which is not annotated in the ucsc known isoform.
higher expression level of an isoform means more reads will be sampled from this isoform to capture the structure of this isoform.
so the predicted isoforms with low expression levels have more probability to be artifacts.
we need an expression level threshold t to exclude those lowly expressed isoforms.
the selection of t is a trade-off between sensitivity and specificity.
the smaller the value t is, the higher sensitivity but lower specificity the result will achieve.
because there is no ground truth for the predicted novel isoforms, we set the predicted isoforms which are annotated by the ucsc known isoforms as positive, others as negative.
to determine this value, we explore the roc curve of the prediction result of mds sample  <dig> where the tcf <dig> in figure  <dig> comes from.
the roc curve and selected threshold of the prediction result of nsmap on the mds sample  <dig> are depicted in figure  <dig>  the optimal threshold t in this roc curve is  <dig> .
the intuition of this optimal threshold is that at this point, the increase rate of sensitivity equals to the increase rate of .
when the threshold is infinite, both sensitivity and  are  <dig>  with the decrease of threshold t, both sensitivity and  will increase.
but the increase rate of sensitivity is larger than that of  until t reduces to the point  <dig> .
after that point, the increase rate of  will be larger than sensitivity.
actually is the false positive rate.
we prefer higher sensitivity and lower false positive rate.
so we select this point as the optimal point.
in the predicted isoforms of gene tcf <dig> in figure  <dig>  the two annotated isoforms have expression levels  <dig>  and  <dig> , respectively, which are above the selected threshold t =  <dig> .
while the expression level of the un-annotated isoform is  <dig>  which is obviously under this threshold.
in this way, we perform a further refinement of the prediction result.
in the following real data analysis, we only consider nsmap predicted isoforms which are larger than the optimal expression level threshold in each sample.
clinical mds sequencing data analysis
the goal is to use our nsmap to identify known and novel isoforms which may be related with mds.
we apply nsmap on the alignment results of the three data sets from tophat to identify the expressed isoforms and their expression levels.
the predicted isoforms are compared with the ucsc annotated reference genes to distinguish the known and novel isoforms.
for the known isoforms, we select the differentially expressed isoforms with fold change greater than  <dig> on both disease samples compared with the control sample.
we also select mds featured novel isoforms which are detected in both mds samples but are absent in the normal sample.
finally, we get total  <dig> differentially expressed known isoforms and  <dig> novel isoforms with rpkm over  <dig>  the two selected isoform sets were fed into ingenuity pathways analysis   <cit>  respectively.
ipa is a database search tool for finding function and pathway for specific biological states.
here we used ipa to explore the pathways enriched in the two selected isoform sets.
the top  <dig> enriched canonical pathways of the two analyses are listed in table  <dig> where mitochondrial dysfunction is both enriched in the known and novel differentially expressed isoform sets which is known to closely related with mds  <cit> .
the oxidative phosphorylation and mitochondrial dysfunction are enriched by both differentially expressed isoform sets which mean the novel differentially expressed isoforms have similar biological functions with the differentially expressed known isoforms.
that observation indicates the prediction of novel isoforms is consistent with the known isoforms.
discussion
through simulations that closely modelled real data, we confirm our method's effectiveness for experiments in both mouse brain and liver rna-seq data.
we also compare nsmap with isoinfer to show that nsmap has comparable performance in identifying known isoforms from rna-seq reads.
finally, we apply nsmap on our mds rna-seq data analysis and find some differentially expressed known isoforms and novel isoform candidates which involve in some mds related pathways.
recently, lacroix et al.
<cit>  showed that unique solution cannot be guaranteed theoretically in isoform identification from short sequence reads.
for example, in our case, all possible transcript isoforms are enumerated according to the detected junction reads.
among these isoforms, one truly expressed isoform may be linear combinations of the other isoforms in terms of exon arrangements.
then the solution of this case is not unique.
the assumption of nsmap to address this issue is that the solution which employs as few expressed isoforms as possible to explain the most observation is preferred.
though this assumption is identical to the assumption made by cufflinks, the implementation of this assumption in nsmap is totally different from cufflinks.
cufflinks constructed a parsimonious set of transcripts followed by the expression level estimations of those constructed transcripts using established expression level estimation model.
however, nsmap enumerates all the possible isoforms formed by the combinations of identified putative exons from tophat and incorporates a prior distribution into the expression level estimation model to control the number of expressed isoforms.
that means the identification of expressed isoforms and the expression level estimations of those identified isoforms are done jointly in nsmap.
paired-end sequencing can dramatically improve the accuracy of isoform level expression estimation which is becoming ubiquitous.
recently, salzman et al.
<cit>  proposed "insert length model" to extend jiang and wong's single-end sequencing work  <cit>  to paired-end sequencing analysis by modeling the insert length distribution.
so we can use this idea to handle paired-end sequencing data in our current framework.
paired-end sequencing can dramatically improve the accuracy of isoform level expression estimation which is becoming ubiquitous.
in paired end sequencing, only the fragments in a specified range will be selected.
several papers have used this information by modeling the fragment length distribution to improve isoform deconvolution problem  <cit> .
salzman et al.
<cit>  proposed "insert length model" to extend jiang and wong's single-end sequencing work  <cit>  to paired-end sequencing analysis by modeling the fragment length distribution.
so we can use this idea to handle paired-end sequencing data in our current framework.
in paired-end sequencing, salzman et al.
defines asi = qn for an event s where the mate reads are mapped into two specified positions on genome.
here fsi is the length of corresponding fragment on the i-th transcript and q is the probability of observing a fragment with length f. in practice, q can be approximated by the empirical probability mass function computed from all the mapped paired-end reads.
in order to reduce the number of events, the minimal sufficient statistics is used to group the events into minimal categories for computational purpose.
in this way, we can incorporate the paired-end information into our model by redefining asi to address paired-end sequencing data.
currently, nsmap uses the tss and pas retrieved from ucsc known isoforms.
we will extend it to identify tss and pas from rna-seq by the following scheme.
if the start point of a putative exon is not a junction point, this putative exon can be regarded as tss.
and if the end point of a putative exon is not a junction point, this putative exon can be regarded as pas.
here junction point means this point is the start or end of a splicing junction.
as our primary motivation is to design a method to identify the isoform structure without annotated reference isoform genome, the usefulness of nsmap is largely dependent on the expression levels of true isoforms and splicing junction detection.
we believe that the accuracy of this approach will increase significantly as the sequencing technology evolves such as paired-end sequencing technique and generates longer sequences with less noise and higher throughput.
CONCLUSIONS
in this paper, we propose a statistical model nsmap for rna-seq data analysis which can be used to identify and quantify isoforms simultaneously without isoform annotations from reference genome.

 more than 90% of human genes are estimated to be alternatively spliced which leads a single gene to produce multiple proteins with distinct functions and is implicated in many diseases including cancer. in recent years, there is an increasing interest in the use of alternative splicing in developing diagnostic tools and in identifying new therapeutic targets. microarrays have been widely used to analyze alternative isoforms by combining exon arrays and exon junction arrays to quantify isoform level expression indexes. however, array based techniques are encountering several fundamental problems such as cross hybridization and weak signals in junction probes which are difficult to overcome. ultra high-throughput sequencing of rna has been developed as an approach for transcriptome analysis in several different species and has offered an attractive approach to measure transcription in a comprehensive manner. rna-seq allows the direct detection of alternative splicing using the reads mapped at splice junctions including the novel splicing without the annotation information. genome-wide measurements of transcriptomes are increasingly done by rna-seq which provides a far more precise measurement of expression levels of isoforms than other methods. the rapidly-developing rna-seq techniques require substantial algorithmic advances. several tools and strategies have been proposed to deal with the complex bioinformatics analysis of rna-seq. pepke et al. provided a comprehensive and up-to-data review of multi-layered analyses of rna-seq data. mortazavi et al. proposed to quantify the gene level expression of a transcript as reads per kilobase per million mapped reads. further, jiang and wong presented a statistical model to describe how the isoform expression levels were calculated from the number of reads mapped to the annotated exons of a gene. meanwhile, bohnert et al. also proposed rquant to determine the abundances for each annotated isoforms by minimizing the deviation of the observation from the expected position-wide read coverage. all these methods assumed that the number and structures of isoforms of each gene are known from the reference genome. however, as jiang and wong pointed out, the isoform level annotation is very incomplete due to the complexity of the transcriptome and the limitations of previous experimental approaches. to address this issue, trapnell et al. proposed cufflinks to identify transcripts as well as to estimate the expression levels of identified transcripts from mapped reads without annotation information. in essence, cufflinks constructs a covering relation on the read alignments from tophat, and find a minimum path cover on the directed acyclic graph for the relation based on dilworth's theorem to construct a parsimonious set of transcripts. after that, the expression levels of those constructed transcripts are estimated using established known isoform expression estimation methods. therefore, the construction of transcripts in cufflinks is independent of the expression level estimation. however, the construction of expressed transcripts and expression level estimation are highly associated. we argue that the determination of parsimonious set of expressed transcripts and expression level estimation should be implemented jointly. though scripture can also detect the novel isoforms, the issue of parsimonious expressed isoforms is not addressed. we also notice that feng et al. proposed isoinfer to identify isoforms using the detected junctions. the candidate isoforms were constructed by combining the putative exons followed by selecting a minimum best subset from all the enumerated subsets of the candidate isoforms which can explain the observation best. however, enumerating all possible subsets of the candidate isoform set with a given size is often infeasible. isoinfer decomposes the large putative exon set into subsets to address this issue which introduces more parameters. here, we put forward nsmap to infer the structures of isoforms as well as to estimate the expression levels simultaneously. first, the exons are constructed based on the detected splicing junctions from rna-seq data using tophat. all the possible isoforms are enumerated by combination of those detected exons. then nsmap is applied to identify the true expressed isoforms from the large candidate isoform set as well as to estimate the expression levels with a sparsity control term to restrict the number of expressed isoforms. the assumption behind this sparsity term is that only as few isoforms as possible should be selected to best explain the observed number of reads falling on each exon of a gene. finally, a model selection step is conducted to select the solution which compromises the fitting of the observation and the number of expressed isoforms best. in summary, our algorithm allows for discovering the structures of the expressed isoforms of a given gene and for estimating the concentration of each spliced isoform simultaneously without the annotated isoform information, which makes the identification of new, previously unknown, alternatively spliced isoforms possible. this study will help rna-seq, a next generation sequencing technology, advance to its full potential in comprehensive transcriptome analysis.  data set we test nsmap on two simulations with simulated expression levels derived from two publicly available mouse rna-seq as described in. we also apply nsmap on three real in-house rna-seq data sets of myelodysplastic syndromes transcriptome analysis to identify isoforms including novel ones featured in mds disease. mds are a diverse collection of hematological conditions united by ineffective production of myeloid blood cells and risk of transformation to acute myelogenous leukemia whose frequency and incidence are increasing in the us population. in our application, cryopreserved marrow cells and paraffin embedded marrow clot sections and marrow core biopsies from mds patients as well as age-matched control sample are being studied by dr. jeff chang's lab at the methodist hospital. these mds patients have been thoroughly evaluated for clinical/morphologic/immunophenotypic data and characterized clinically by transfusion dependency and pathologically by significant dysplasia, increased blasts, and immunophenotypic aberrancy. the control sample are obtained from patients without cytopenias. we specifically selected these controls to be age-matched for the mdss population to control for the possibility of aging-related changes in the expression profile mrna of hematopoietic cells. then we apply the rna-seq protocol to sequence our samples. we sequenced the two mds clinical samples and one normal sample using illumina genome analyzer ii. there are around million single-end reads with read length bp for each sample. algorithm summary nsmap comprise four consecutive steps, starting with junction detection and reads mapping using tophat and followed by the candidate isoforms construction, expressed isoform identification and expression level estimation along whole regularization path and model selection to select the best solution from the whole solution path. short reads alignment is the first step in understanding next-generation sequencing data and many free alignment software packages are available. here we use tophat to perform the alignment task which can detect the junctions and map a massive amount of reads to the whole genome flexibly and efficiently. the reference genome sequences are downloaded from ucsc genome database. after read alignment, the next step is to generate the candidate isoforms according to the alignment and splice junctions obtained from tophat. candidate isoforms construction based on the alignment results and detected junctions from tophat, the exons can be constructed from segments whose two ends have been detected as junction points. for example, figure mimics a gene with nine exons constructed from the detected junctions in a gene region. the isoforms are formed from combination of those exons. each rectangle represents an exon and arrow line indicates that there is a splice junction detected between the two exons. this can be interpreted as a directed graph with each exon as nodes. each candidate isoform can be generated by finding a path between two nodes on this graph. for instance, figure displays the possible isoforms with exon and as source and sink respectively. there are total four paths connecting exon and as indicated in figure isoform expression level estimation given the candidate isoforms constructed from last step, we need to select the expressed isoforms from this candidate set as well as to estimate the express level of selected isoforms. we adapted the expression level estimation framework of jiang and wong by incorporating a sparsity regularization term. for a gene g, suppose it has m exons with lengths and n isoforms with expression Φ =. assuming each exon can be either included in an isoform or not, we have a set of observations, where is an index set of events which we are interested in. each observation xs ∈ x is a random variable representing the number of reads falling into a certain region of interest in gene g. for example, reads falling into certain exon or exon-exon junction. the natural statistical model of count data is the poisson distribution. each xs ∈ x follows poisson distribution with a parameter λ which is the expected count or the mean parameter of the poisson distribution. for instance, the λ for the number of reads falling into exon j is, where n is the total number of mapped reads and cij is if isoform i contained exon j and otherwise. for a exon-exon junction event, the λ is, where l is the length of the junction region, j and k are indexes of the two exons involved in the junction being investigated. in general, λs is a linear function of ϕ ϕ..,ϕn, i.e. where asi is a known coefficient. the likelihood function is maximum likelihood estimation method can be used to estimate Φ by maximizing above likelihood function. however, in the work of jiang and wong, the structure of individual isoform is fetched from the reference genome where only a few of isoforms are known for each gene. in our application, the structures of expressed isoforms are unknown which need to be selected from the large candidate isoform set where each isoform is formed by combination of detected exons. to do this, the most important prior information is arguably sparsity; that is, characteristic pattern of each gene in cells are expressed by very few isoforms, even though the total number of possible combinations of exons in a gene can be large. utilization of this sparsity information is crucial to the success of finding the most relevant isoforms which explain the observed reads best, in the face of "insufficient" data and/or data uncertainty. to incorporate this sparsity to our model, we propose a prior laplace distribution on Φ as below: laplace distribution assigns higher probability at zero than gaussian distribution to produce sparser solution and has been widely used to encode sparseness prior. hence we can employ maximum a posteriori instead of maximum likelihood to estimate Φ: taking logarithm on the above equation, we get keep in mind that the expression of each isoform cannot be negative, that is ϕi ≥, f can be rewritten as: where ϕi ≥ for i =...,s and s is the number of events with which we are interested. maximizing f with respect to Φ is equivalent to minimizing to that end, the optimization problem is summarized as: we call this nonnegativity and sparsity constrained map model as nsmap. laplace prior distribution has a l norm term to impose the sparsity. in some cases of our experiment, we find it may not work well in the identification of the correct expressed isoforms due to the similarities between the isoforms. the ideal sparsity constraint is the l regularization, defined as the number of non-zero entries in a vector. unfortunately, the computation of l regularization is intractable to solve. l norm is popular because of its intrinsic convex property. some recent works suggest that non-convex regularizations such as lp regularization with < p < have better performances in parameter selection and sparse signal construction than l regularization. this inspires us to use a stronger sparsity constraint with lp, <p < norm to solve our problem. we will implement a simulation study to compare the performances of sparsity regularizations with different p. in our experiment, we set p =. so the objective function is finally given by: where and. then the expression levels of all the candidate isoforms can be obtained by optimizing equation with a given t. when t → +∞, the solution equals zero without expressed isoforms. with the decreasing of t, some elements of will be non-zero to become expressed isoform. so we have to select an optimal t. before doing this, we will calculate the solution path which consists of solutions corresponding to different values of t. blasso is adapted to approximate this solution path efficiently. model selection from solution path after getting the solution path which consists of solutions corresponding to different values of t, we have to select the best solution from this solution path which makes a good balance between the number of expressed isoforms and the fitting function in equation. the solution path will first be grouped into subsets with increasing model size according to the number of expressed isoforms of each solution. here the model size means the number of expressed isoforms. then the solution with minimal in each group is selected as the best solution for each model size. starting from model size we compare the of current model size with the next model size which have one more expressed isoforms. if is significantly improved by the larger model size, the model with larger model size will be updated to the current model and compared with the remaining models with larger model size than the updated model. otherwise the solution corresponding to the current model size is selected as the final solution. in this way, the solution with smaller number of expressed isoforms will be selected preferably. see the methods for more detail. simulation on the whole genome with different lp norm in equation in the absence of rna-seq data from samples for which we have ground truth isoform quantities, we conduct simulations to validate our method and evaluate its performance in terms of isoform identification and expression levels estimation with p = and p =. we first derive expression level of each isoform from the mouse brain and liver rna-seq data sets described in using the poisson model of with the mouse ucsc genes as the annotated reference genome database. those derived expression levels are employed as the ground truth to perform the following simulations. which exon or splice region the read will fall on is determined by uniformly sampling proportional to the simulated ground truth rpkm and the same mapped reads number of each gene in the real data. after uniformly sampling and counting reads number, we can identify isoform structures and their concentrations from the simulated rna-seq reads data and evaluate the performance of nsmap with the simulated ground truth. we only consider isoforms with expression levels no less than rpkm as expressed isoforms. in this simulation, constructed isoforms containing more than half of the total exons are eligible as the candidate isoforms. the accuracy of nsmap largely depends on the critical step of candidate isoform generation. so, we check whether our candidate isoforms generated from the exon reads and exon junctions can cover the true expressed isoforms. from table we see that more than 96% of expressed isoforms in the two data sets are included in the candidate isoforms generated by our strategy which finds paths on a graph constructed from exons and splice junctions. isoforms included in we further compute the fraction of isoforms for which the estimates are significantly consistent with the simulated ground truth. we refer to this statistic as the positive fraction. given the positive isoforms which are within the 5% deviation from the ground truths, overall specificities of the two simulations are also calculated by setting the truly expressed isoforms with rpkm no less than as positive isoforms and the false reported isoforms by nsmap and truly expressed isoforms with rpkm < as negative isoforms. the results are summarized in table we can see that nsmap achieved near 80% overall positive fraction in both simulations. and highly expressed isoforms have high positive fraction. while for isoforms with low rpkm, the accuracy is not very good because the number of reads falling on this isoform is small which may not be sufficient to capture the exons and splice junctions to form this isoform. this phenomenon is consistent with other studies that the estimation on the low abundance isoform is not very accurate. comparing the performances of nsmap with p = and p = we notice that lp norm with p = provides better results than p = in the simulations. to explain this observation, we conduct the following experiment and give a mathematical interpretation of the two lp norms. demonstrations of features of lp norm with p = and p = 1 first, we exemplify our observation in the above simulation that lp norm with p = is better than p = through a simulation using the structure of gene eif5a. in the ucsc genome annotation, gene eif5a has exons and known isoforms. we set both rpkm expression levels of the fifth and sixth isoforms as because the two isoforms express most from the estimation using the method of on the mouse brain and liver rna-seq data in. the simulation is performed based on the expression levels and structures of the two expressed isoforms under the uniform distribution for reads. based on the simulation data, total candidate isoforms is constructed. the indexes of the two truly expressed isoforms in the candidate set are and respectively. from the upper row of figure we can see that the two expressed isoforms are selected correctly when p = 1/ while in the case of p = additional two false unexpressed isoforms are also selected. then, the lower row of figure shows the mathematical properties of the two norms. l1/ norm regularization has more similar sparsity property with l norm than l norm because the points with high probabilities in the laplace-like distribution with l1/ norm are more focusing around the axes than those in the laplace distribution. this means l1/ norm regularization imposed stronger sparsity than l norm. the simulation depicts the stronger sparsity constraint by setting p = 1/ is superior to laplace priori distribution with p = in our application. in the following experiments, we select the lp norm with p = 1/ comparison with isoinfer on two rna-seq data sets of mds samples here we compare nsmap with isoinfer which have similar ideas in isoform construction from putative exons and minimum expressed isoform set selection. however, isoinfer selects a subset as expressed isoforms from candidate isoform set by enumerating all the possible subsets of the candidate isoform set. in nsmap, the selection of expressed isoforms is embedded into the isoform expression level estimation framework by incorporating a sparsity control term. a transcript can be constructed from all the exon-intron boundaries as well as the transcription start site and poly-a site of an isoform. the exon-intron boundaries can be inferred from rna-seq using alternative splicing detection tool, such as tophat and splicemap. the tss and pas represent the start and end expressed segments of a transcript, respectively. theoretically, any expressed segments can be the tss or pas which will introduce many false short isoforms to make isoform inference difficult. we prefer to retrieve the tss-pas from the ucsc known isoform table as the starts and ends of predicted transcripts and to identify isoforms within the regions of known genes whose functions and pathways are intensively studied. isoinfer and nsmap will use those tss-pas and the detected junctions using tophat to infer the expressed isoforms from rna-seq. because it is infeasible to validate all the predicted isoforms, we evaluate the two methods by comparing their predictions with ucsc known isoform data set. the performance of the method is measured by sensitivity and precision. here we use hg known human isoforms data set downloaded from ucsc which contains transcripts. a known isoform is identified if it is in the prediction result of a method. sensitivity is defined as the number of identified isoforms divided by the number of all known isoforms from ucsc data base. precision is defined as the number of identified known isoforms divided by the number of predicted isoforms by the method. however, isoinfer selects a subset as expressed isoforms from candidate isoform set by evaluating all the possible subsets of the candidate isoforms. in nsmap, the selection of expressed isoforms is embedded into the isoform expression level estimation framework by incorporating a sparsity control term. in this way, the selection of expressed isoforms is automatic and more efficient than testing all the possible sub sets of the candidate isoforms. we also notice that more known isoforms are predicted in mds sample than mds sample because there are more reads are mapped in mds sample this observation tells us that deeper sequencing will improve the performances of isoinfer and nsmap. here, the sensitivities of the two methods are very low. the reason is that we compare the predicted isoforms with the large ucsc known isoform set. some of the ucsc known isoforms may not express in the sample. so the effective sensitivities will be larger than those numbers. this comparison against ucsc known isoform data set does not mean all the predicted isoforms without annotation are false. especially those predicted novel isoforms with high rpkm are promising to be true novel isoforms. for example, if we select the predicted isoforms with expression level larger than rpkm, in mds sample and out of and out of predicted isoforms by nsmap will be annotated in the ucsc known isoforms table. and the lowly expressed isoforms have higher probability to be false positive or artifacts due to the insufficient reads for capturing the true structure of an isoform. so we need to set an expression level threshold to refine the predicted isoforms. this issue is addressed in the following section. example of identified isoforms using nsmap and expression level threshold selection here we exemplify isoforms of a gene tcf estimated by nsmap. tcf encodes protein transcription factor which localizes to the nucleus and displays dna-binding and transactivation activities. tcf has five exons and alternative splicing results in two known transcript variants encoding different isoforms in the ucsc reference genome ). the rna-seq data is first processed by tophat to detect the splicing junctions. figure shows the splicing junctions of tcf on mds sample detected by tophat. then nine exons are constructed based on detected junctions in figure nsmap predicts three isoforms with expression level, and respectively ). compared with the annotated isoforms of tcf ), we confirm that the top two highly expressed predictions are the same with the annotations. this example demonstrates the ability of nsmap in isoform identification. however, nsmap still predicts an additional lowly expressed isoform which is not annotated in the ucsc known isoform. higher expression level of an isoform means more reads will be sampled from this isoform to capture the structure of this isoform. so the predicted isoforms with low expression levels have more probability to be artifacts. we need an expression level threshold t to exclude those lowly expressed isoforms. the selection of t is a trade-off between sensitivity and specificity. the smaller the value t is, the higher sensitivity but lower specificity the result will achieve. because there is no ground truth for the predicted novel isoforms, we set the predicted isoforms which are annotated by the ucsc known isoforms as positive, others as negative. to determine this value, we explore the roc curve of the prediction result of mds sample where the tcf in figure comes from. the roc curve and selected threshold of the prediction result of nsmap on the mds sample are depicted in figure the optimal threshold t in this roc curve is. the intuition of this optimal threshold is that at this point, the increase rate of sensitivity equals to the increase rate of. when the threshold is infinite, both sensitivity and are with the decrease of threshold t, both sensitivity and will increase. but the increase rate of sensitivity is larger than that of until t reduces to the point. after that point, the increase rate of will be larger than sensitivity. actually is the false positive rate. we prefer higher sensitivity and lower false positive rate. so we select this point as the optimal point. in the predicted isoforms of gene tcf in figure the two annotated isoforms have expression levels and, respectively, which are above the selected threshold t =. while the expression level of the un-annotated isoform is which is obviously under this threshold. in this way, we perform a further refinement of the prediction result. in the following real data analysis, we only consider nsmap predicted isoforms which are larger than the optimal expression level threshold in each sample. clinical mds sequencing data analysis the goal is to use our nsmap to identify known and novel isoforms which may be related with mds. we apply nsmap on the alignment results of the three data sets from tophat to identify the expressed isoforms and their expression levels. the predicted isoforms are compared with the ucsc annotated reference genes to distinguish the known and novel isoforms. for the known isoforms, we select the differentially expressed isoforms with fold change greater than on both disease samples compared with the control sample. we also select mds featured novel isoforms which are detected in both mds samples but are absent in the normal sample. finally, we get total differentially expressed known isoforms and novel isoforms with rpkm over the two selected isoform sets were fed into ingenuity pathways analysis respectively. ipa is a database search tool for finding function and pathway for specific biological states. here we used ipa to explore the pathways enriched in the two selected isoform sets. the top enriched canonical pathways of the two analyses are listed in table where mitochondrial dysfunction is both enriched in the known and novel differentially expressed isoform sets which is known to closely related with mds. the oxidative phosphorylation and mitochondrial dysfunction are enriched by both differentially expressed isoform sets which mean the novel differentially expressed isoforms have similar biological functions with the differentially expressed known isoforms. that observation indicates the prediction of novel isoforms is consistent with the known isoforms. discussion through simulations that closely modelled real data, we confirm our method's effectiveness for experiments in both mouse brain and liver rna-seq data. we also compare nsmap with isoinfer to show that nsmap has comparable performance in identifying known isoforms from rna-seq reads. finally, we apply nsmap on our mds rna-seq data analysis and find some differentially expressed known isoforms and novel isoform candidates which involve in some mds related pathways. recently, lacroix et al. showed that unique solution cannot be guaranteed theoretically in isoform identification from short sequence reads. for example, in our case, all possible transcript isoforms are enumerated according to the detected junction reads. among these isoforms, one truly expressed isoform may be linear combinations of the other isoforms in terms of exon arrangements. then the solution of this case is not unique. the assumption of nsmap to address this issue is that the solution which employs as few expressed isoforms as possible to explain the most observation is preferred. though this assumption is identical to the assumption made by cufflinks, the implementation of this assumption in nsmap is totally different from cufflinks. cufflinks constructed a parsimonious set of transcripts followed by the expression level estimations of those constructed transcripts using established expression level estimation model. however, nsmap enumerates all the possible isoforms formed by the combinations of identified putative exons from tophat and incorporates a prior distribution into the expression level estimation model to control the number of expressed isoforms. that means the identification of expressed isoforms and the expression level estimations of those identified isoforms are done jointly in nsmap. paired-end sequencing can dramatically improve the accuracy of isoform level expression estimation which is becoming ubiquitous. recently, salzman et al. proposed "insert length model" to extend jiang and wong's single-end sequencing work to paired-end sequencing analysis by modeling the insert length distribution. so we can use this idea to handle paired-end sequencing data in our current framework. paired-end sequencing can dramatically improve the accuracy of isoform level expression estimation which is becoming ubiquitous. in paired end sequencing, only the fragments in a specified range will be selected. several papers have used this information by modeling the fragment length distribution to improve isoform deconvolution problem. salzman et al. proposed "insert length model" to extend jiang and wong's single-end sequencing work to paired-end sequencing analysis by modeling the fragment length distribution. so we can use this idea to handle paired-end sequencing data in our current framework. in paired-end sequencing, salzman et al. defines asi = qn for an event s where the mate reads are mapped into two specified positions on genome. here fsi is the length of corresponding fragment on the i-th transcript and q is the probability of observing a fragment with length f. in practice, q can be approximated by the empirical probability mass function computed from all the mapped paired-end reads. in order to reduce the number of events, the minimal sufficient statistics is used to group the events into minimal categories for computational purpose. in this way, we can incorporate the paired-end information into our model by redefining asi to address paired-end sequencing data. currently, nsmap uses the tss and pas retrieved from ucsc known isoforms. we will extend it to identify tss and pas from rna-seq by the following scheme. if the start point of a putative exon is not a junction point, this putative exon can be regarded as tss. and if the end point of a putative exon is not a junction point, this putative exon can be regarded as pas. here junction point means this point is the start or end of a splicing junction. as our primary motivation is to design a method to identify the isoform structure without annotated reference isoform genome, the usefulness of nsmap is largely dependent on the expression levels of true isoforms and splicing junction detection. we believe that the accuracy of this approach will increase significantly as the sequencing technology evolves such as paired-end sequencing technique and generates longer sequences with less noise and higher throughput.  in this paper, we propose a statistical model nsmap for rna-seq data analysis which can be used to identify and quantify isoforms simultaneously without isoform annotations from reference genome. 
most existing methods employ rna-seq to quantify the expression levels of already annotated isoforms from the reference genome. a nonnegativity and sparsity constrained maximum aposteriori model has been proposed to estimate the expression levels of isoforms from rna-seq data without the annotation information. in contrast to previous methods, nsmap performs identification of the structures of expressed isoforms and estimation of the expression levels of those expressed isoforms simultaneously, which enables better identification of isoforms. then, we apply nsmap on two rna-seq data sets of myelodysplastic syndromes samples and one normal sample in order to identify differentially expressed known and novel isoforms in mds disease. nsmap provides a good strategy to identify and quantify novel isoforms without the knowledge of annotated reference genome which can further realize the potential of rna-seq technique in transcriptome analysis. novel study on isoform inference and estimation purely from rna-seq without annotation information is desirable. however, the current reference genome is very incomplete due to the complexity of the transcriptome which hiders the comprehensive investigation of transcriptome using rna-seq. in the simulations parameterized by two real rna-seq data sets, more than 77% expressed isoforms are correctly identified and quantified.
 the development of techniques for sequencing the messenger rna enables it to study the biological mechanisms such as alternative splicing and gene expression regulation more deeply and accurately. most existing methods employ rna-seq to quantify the expression levels of already annotated isoforms from the reference genome. however, the current reference genome is very incomplete due to the complexity of the transcriptome which hiders the comprehensive investigation of transcriptome using rna-seq. novel study on isoform inference and estimation purely from rna-seq without annotation information is desirable.  a nonnegativity and sparsity constrained maximum aposteriori model has been proposed to estimate the expression levels of isoforms from rna-seq data without the annotation information. in contrast to previous methods, nsmap performs identification of the structures of expressed isoforms and estimation of the expression levels of those expressed isoforms simultaneously, which enables better identification of isoforms. in the simulations parameterized by two real rna-seq data sets, more than 77% expressed isoforms are correctly identified and quantified. then, we apply nsmap on two rna-seq data sets of myelodysplastic syndromes samples and one normal sample in order to identify differentially expressed known and novel isoforms in mds disease.  nsmap provides a good strategy to identify and quantify novel isoforms without the knowledge of annotated reference genome which can further realize the potential of rna-seq technique in transcriptome analysis. nsmap package is freely available at https://sites.google.com/site/nsmapforrnaseq. 