BACKGROUND
understanding of the sugarcane transcriptome is limited due to the complexity in gene copy number, repetitive content, and heterozygosity in the genome  <cit> .
it is not clear how many transcript isoforms result from the alternative splicing in this potentially very complex transcriptome.
sugarcane is a polyploid hybrid between saccharum officinarum and s. spontaneum, and each sugarcane hybrid has its own unique chromosome set , containing up to  <dig> copies of each gene and a total ~ <dig>  predicted genes  <cit> .
therefore, it is expected that sugarcane transcripts represent transcription of genes/homoelogues that are not only unique to the progenitor genomes but also transcription from alternate splicing, which we collectively refer to as transcript isoforms.
sugarcane genes such as those in the sucrose phosphate synthase  gene families  <cit> , invertase genes  <cit>  and sucrose synthase family  <cit>  have been shown to be comprised of many isoforms.
most of the sugarcane studies, including transcriptome studies found in the literature, i.e.
in  <cit>  and  <cit> , have been based on sorghum genomic/transcript sequences  <cit>  which have the highest gene synteny and orthologous alignment with the sugarcane genome  <cit> ; sugarcane expressed sequence tags   <cit> ; saccharum officinarum gene indices - sogi v <dig>   <cit>  representing ~90% of the estimated genes in s. officinarum  <cit> ; and other resources reviewed in  <cit> .
studies based on these databases have provided useful information on the sugarcane transcriptome, while a whole genome sequence is not yet available.
however, it is thought that there are still many sugarcane genes missing in these databases  <cit>  and in addition, the full-length  sequences of distinct transcript isoforms are not included.
use of these transcript databases for rna-seq analysis leading to the identification differentially expressed genes does not provide information on the corresponding isoform/s or the homoelogue/s contributing to the differential expression.
there is a need to construct fl transcript sequences including such isoforms to facilitate analysis of isoform differential expression, and also to extend our understanding of the sugarcane transcriptome.
the transcriptome poses a great challenge when it comes to assembly and annotation.
the differences in transcript abundance and the presence of different isoforms, greatly challenge the assembly of a transcriptome from short-reads ; since the assemblers cannot distinguish between reads originally from different transcripts/isoforms carrying the same exons  <cit> .
to date, most sequencing platforms offer a read-length which is shorter than the typical length of a eukaryotic mrna   <cit> .
the transcriptome sequences obtained from second generation sequencing technology  have been playing an important role in capturing the diversity in the rna populations at a greater sequencing depth  <cit> .
however, a precise prediction and identification of the alternative transcript splicing has not been possible.
algorithms in transcript splice-aware assemblers  have been developed to detect splicing junctions and recover transcript isoforms by using information from short-reads, but these have not always been confirmed.
that is, quite often these approaches overestimate and report spurious computational isoforms rather than picking up only biological ones.
overall, the assemblies from short-read data normally end up identifying more transcripts than expected , which may be attributed not only to the diversity of rnas and diversity of transcript/isoforms in the transcriptome, but also to the limitation in recovering fl transcripts.
most studies use these tools, then filter the transcripts through clustering by retaining the longest sequence in each cluster as representative for analysis, and consider them to be the major isoforms or unigenes  <cit> .
current algorithms such as in bernard et al.
for isoform identification and quantification require longer reads and cannot tackle genes with too many exons.
with the advent of third generation sequencing technology, the cost-per-transcriptome has been reduced, whereas the length of the sequencing reads has been increased significantly.
as of august,  <dig>  the average read length of pacbio single molecule, real-time  sequencing is >10 kb and real length can be up to 60 kb .
this protocol has been applied in some recent studies, for example, detecting  <dig>  alternative splicing events in  <dig>  unique transcripts , covering ~89% of the total sorghum annotated genes  <cit> ; and producing  <dig>  unique transcripts  in the maize transcriptome derived from six different tissues, covering ~70% of the annotated genes  <cit> .
this study represents the first full-length transcriptome reference sequences from sugarcane derived from three different tissues, of different developmental stages, by using the pacbio long-read iso-seq technique.
in addition, rna-seq was used to improve the pacbio transcript isoforms by short-read error correction, and comparison between sugarcane transcripts obtained/assembled from these two different platforms.
annotation of the sugarcane fl transcript isoforms could improve sugarcane genome models, contribute towards understanding of the complexity of the sugarcane genome and serve as reference sequences for differential expression analysis in the future.
RESULTS
sugarcane transcriptome from pacbio isoform sequencing
a pooled sample representing polya rnas from three tissues , of different developmental stages  was sequenced to obtain a wide coverage of the sugarcane transcriptome.
a total of  <dig>  reads of inserts  was generated, with a total of  <dig> , <dig> nucleotides from six smrt cells of non-normalized bins  and normalized bins , including  <dig>   fl non-chimeric rois and  <dig>   non-fl rois.
the length of rois ranged from 300 bp to  <dig>  kb, with an n <dig> of  <dig>  bp.
the length distribution of all roi data is presented in additional file 1: figure s <dig>  a total of  <dig>  high quality sequences and  <dig>  low quality sequences were obtained from quiver polishing, referred to as polished transcript isoforms.
the total unique, non-redundant transcript isoforms included  <dig>  sequences, with the length ranging from 301 bp to  <dig>  bp, n <dig> of  <dig>  bp and n <dig> of  <dig>  bp and  <dig> % gc content.
improving pacbio transcript quality by error correction using rna-seq reads
we followed two error correction pipelines , using three datasets from rna-seq derived from the same experiment;  <dig> , <dig> non-normalized reads,  <dig> , <dig> reads bbnorm-normalized and  <dig> , <dig> trinity-normalized reads.
overall, the error correction led to improvement in transcript prediction, more transcripts covered the full-length of known proteins, longer open reading frames , better completeness results in cegma/busco assessments, and higher alignment rate of transcript isoforms to the sorghum genome.
this only resulted in a slight change in the total number of transcripts isoforms after removing all exact 100% identical sequences.
the lordec error correction outperformed proovread with our transcript data.
the best corrected set was derived from lordec using the trinity-normalized reads, which resulted in  <dig> % of the total transcripts with orfs passing the evigene score in transcript prediction, while that of the non-corrected transcripts was only  <dig> % .
there were  <dig>  orfs  detected by transdecoder in this corrected pacbio transcript isoform set, compared to  <dig>  orfs  for the non-corrected dataset.
the retained orfs with pfam and viridiplantae protein hits from the corrected dataset had an n <dig> of  <dig>  bp, while that of the non-corrected was 708 bp.
the cegma and busco alignments showed that the corrected pacbio dataset had a higher completeness level than the non-corrected .
about  <dig> % of the corrected transcripts aligned to the sorghum genome, while the alignment of the non-corrected transcript isoforms was  <dig> %.
we chose the pacbio dataset corrected by lordec, using reads normalized by trinity , for downstream analysis.
this final pacbio transcript isoforms set had  <dig>  sequences after removal of strictly identical sequences; of total length ~193 mb, with individual transcript length ranging from  <dig> to  <dig>  bp, n <dig> of  <dig> , n <dig> of  <dig> , and  <dig> % gc content.table  <dig> summary of correction of pacbio transcript isoform data using illumina short-reads
de novo assembly of the sugarcane transcriptome from short-reads
the assembly of the sugarcane transcriptome from illumina rna-seq short-reads was carried out to provide a comparative reference for the transcript isoform sequences obtained from pacbio iso-seq, since rna-seq has been utilized widely in construction of transcriptomes.
of  <dig>  million total raw reads generated ,  <dig> , <dig>  reads survived after trimming, having a quality score cutoff of q <dig>  trinity normalization retained only 6%  at a maximum coverage of  <dig>  while retaining 15%  and 21%  at maximum coverage of  <dig> and  <dig> , respectively.
a qc report of all read datasets was generated by fastqc v <dig> .
<dig>  <cit>  .
the use of bbnorm was selected due to the fact that trinity normalization heavily reduced the reads compared to bbnorm at the same maximum coverage cutoff.
about 37%  of the total reads remained at maximum coverage  <dig>  by bbnorm package for de novo assembly.
after assembly and individual clustering, four initial assemblies were obtained from trinity, clc-gwb, velvet/oases, and soapdenovo-trans, respectively.
we observed varied contig number, n <dig>  cumulative length and length distribution in each of the assemblies.
the total number of contigs from the trinity-assembly was  <dig>  , while that from the clc-gwb assembly was  <dig>  , velvet/oases assembly gave  <dig>   and soaptrans assembly  <dig>  .
table  <dig> summarizes the de novo assembly results in this study, including all major statistics and qc, more details see additional file 1: figure s <dig>  the final total number of contigs clustered by cd-hit-est at 95% identity and after retaining transcripts with length from 300 bp to 10 kb, was  <dig> , of ~967 mb, having an n <dig> of  <dig>  bp, n <dig> of 745 bp and  <dig> % gc.
this clustered assembly was referred to as the de novo transcript contigs.table  <dig> comparison of de novo assemblies used in this study
when compared to pacbio transcript isoforms  by blastn ,  <dig> % of the de novo transcript contigs  exhibited similarity to  <dig> % of the pacbio transcript isoforms .
there were  <dig> % of de novo transcript contigs  and  <dig> % of pacbio transcript isoforms  that were unique to each of the datasets.fig.
<dig> comparison between the sugarcane pacbio transcript isoforms and de novo transcript contigs
analysis of reads mapping back to transcripts
in mapping back the rna-seq read data to transcripts from both pacbio and de novo transcripts, we observed  <dig> % of reads mapped back to pacbio transcript isoforms, while  <dig> % mapped to the de novo transcript contigs.
figure  <dig> shows the average coverage plotted against the transcript length for both assemblies.fig.
<dig> average coverage of sugarcane de novo contigs and pacbio isoforms obtained from read mapping.
a, coverage of de novo transcript contigs.
b, coverage of pacbio transcript isoforms
transcriptome completeness analysis
in both cegma and busco alignments , the pacbio dataset showed a lower completeness level than the de novo dataset.
the pacbio transcripts had  <dig> % cegma alignment , and the de novo assembly had  <dig> % cegma alignment .
there were no missing cegma in the de novo transcript contigs, and there was  <dig> % missing cegmas  in the pacbio dataset.
similarly, in busco alignment to  <dig> conserved proteins, the pacbio transcript isoforms had lower completeness than that of de novo assembly, by having  <dig> % completeness compared to  <dig> % .
the de novo assembly had a higher level of duplication in the busco alignment, suggesting that the assembly contains duplicate contigs of different lengths  assembled by different assemblers and retained after clustering.
using the same cegma and busco protein alignments, we assessed the unigene dataset from  <cit>  and the sogi database to determine the consistency of the methods.
the unigenes had  <dig> %  cegma completeness,  <dig> %  busco completeness; and the sogi dataset had  <dig> %  cegma completeness and  <dig> %  busco alignment completeness.
the sogi dataset had the largest proportion of partial/fragmented  and missing proteins , since this dataset contained gene indices and ests .table  <dig> transcriptome coverage analysis based upon cegma and busco alignment
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
counting the full-length and nearly full-length transcripts
the number of transcripts appearing to be full-length with at least 90% and 70% alignment coverage of the viridiplantae uniprot proteins was estimated and compared between pacbio transcript isoforms and de novo transcript contigs .
the pacbio dataset had  <dig>  transcripts appearing to be full-length , and  <dig>  transcripts , and that of the de novo transcript contigs were  <dig>  and  <dig>  at 90% and 70%, respectively.
analysis of the matched proteins at the cutoff of 70% alignment coverage from both assemblies indicated that only  <dig>  proteins  were common between these two assemblies, with  <dig>   unique to pacbio dataset, and  <dig>   unique to the de novo transcript contigs.fig.
<dig> full-length analysis between sugarcane pacbio transcript isoforms and de novo transcript contigs.
a, counts of proteins covered by transcripts at different thresholds.
b, comparison between the protein hits from pacbio and de novo transcripts which covered at least 70% viridiplantae protein length
as these results considered each protein from the uniprot database as only one count, regardless of the presence of different isoforms that carry the same protein sequence part  getting matched many times.
a modified approach to counting the full-length transcripts for isoforms was applied, in which we took all counts of isoforms that hit the same protein into consideration and estimated the number of full-length proteins covered.
using this strategy, we found in the pacbio data,  <dig>  transcripts that covered ≥90% of viridiplantae proteins and  <dig>  transcripts that covered ≥70% of viridiplantae proteins.
in the de novo transcript contigs, it was  <dig>  and  <dig>  protein hits covered by ≥90% and ≥70%, respectively.
de novo transcript contigs had more proteins covered at lower percentage due to the greater duplication retained in the assembly, and inclusion of more partial gene content.
when protein hits of ≥90% alignment coverage from the two results were compared, the unique protein hits of pacbio was  <dig>  and de novo was  <dig> , which were the same as in the first approach.
investigation of  <dig> full-length genes from sugarcane and other grass species showed that pacbio dataset resulted in a better performance in terms of recovering the full-length sequence of these specific genes .
at an e-value ≤1e- <dig>  there were more genes covered by transcripts at 90%  and 70%  in pacbio transcript isoforms than that in de novo transcript contigs , unigene set  and sogi database , respectively.
the lower full-length gene count in unigenes could be due to only main isoforms being retained for this dataset, while lower full-length gene count in sogi database could be due to the fraction of ests in it.table  <dig> alignment and full-length assessment of a selected gene set
evidence of alternative splicing in the sugarcane transcriptome from pacbio iso-seq
using the results from pacbio transcript isoforms mapped against the sorghum genome , our in-house sugarcane whole genome assembly from sugarcane cultivar q <dig>  and particular contigs that spanned through the sucrose phosphate synthase a and cellulase  <dig> genes , we were able to visualize alternative splicing in the pacbio transcript isoforms.
using the transcriptome analysis pipeline for isoform sequencing  described in  <cit> , we detected amongst those transcript isoforms aligned against the sorghum genome,  <dig>  alternative splicing events, including  <dig>   intron retention,  <dig>  skipped exon,  <dig>   alternative 5’splice sites and  <dig>   alternative 3’splice sites.
an estimation of exons per transcripts amongst the transcript isoforms aligned against the sorghum genome is also included .fig.
<dig> evidence of different transcript isoforms of sugarcane transcriptome present in the pacbio transcript dataset.
a, isoforms aligned against the sorghum chromosome  <dig>  b, isoforms aligned to contigs of our in-house sugarcane whole genome de novo assembly.
c, different transcript isoforms aligned to sucrose phosphate synthase gene and cellulase  <dig> gene contigs.
d, average exons per transcript estimated based on the transcript isoforms aligned against sorghum genome
prediction of potential coding regions, main and alternate transcripts analysis
we analyzed the candidate coding regions in the transcript sequences by retaining only open reading frames  that exhibited homology with the pfam protein domain database or the uniprot viridiplantae known protein database, which were more likely to be biologically real.
there were  <dig>  non-overlapped orfs detected in the pacbio transcript isoform sequences, belonged to  <dig>  orf-containing transcript isoforms .
only  <dig>  isoforms did not contain orfs and these were used for characterization of non-coding rnas, in the next section.
of the total orf-containing transcripts,  <dig>  matched the viridiplantae proteins , while  <dig>  matched the pfam protein domains.
the total number of transcript sequences retained in combination with the transdecoder frame-score was  <dig>  , with lengths ranging from 300 bp to  <dig>  bp, and n <dig> of  <dig>  bp.
there were more sequences retained in the de novo assembly, since it started with more data, however transcript contig length was shorter than that of the pacbio dataset.
a total  <dig>  orfs were found in  <dig>  orf-containing contigs .
when combined with the homology search results,  <dig>  transcript contigs matched the viridiplantae protein database, and  <dig>  contigs matched the pfam protein domains.
the final number of contigs retained by transdecoder was  <dig> , accounting for  <dig> % of the total de novo assembly.
this final contig set had lengths ranging from 300 bp to  <dig>  bp and n <dig> of 738 bp.
when predicting the potential coding genes using the evigene pipeline, the total number of predicted transcripts in the pacbio transcript isoforms was  <dig>  , while the dropset had  <dig>  transcripts .
the average length of the largest  <dig>  proteins from the dataset was reported to be  <dig>  aa, and all candidate transcripts  had an n <dig> of  <dig>  bp and cds length ranged from 186 bp to  <dig>  bp.
as for the prediction by transdecoder, the de novo contig set had more orfs, and therefore, more predicted transcripts both main and alternate.
there were  <dig>  predicted transcripts  with an n <dig> of  <dig>  compared to pacbio transcript isoforms, the de novo predicted transcripts had much shorter length distribution and average length of the largest  <dig>  proteins .
using the same prediction approach on the unigenes and sogi, we found  <dig>  predicted main transcripts in the unigenes , and  <dig>  predicted transcripts  in sogi.
the average length of the largest  <dig>  proteins for the unigenes was  <dig> aa, while that of sogi was  <dig> aa.
all the results of orf and transcript prediction are presented in table  <dig> and fig.
<dig> table  <dig> open reading frame and transcript prediction analysis of sugarcane transcriptome sequence data
retained transcriptsa
ave length 1 k proteinsb
atranscripts with pfam and viridiplantae hits.
baverage length  of the largest  <dig>  proteins
fig.
<dig> analysis of orfs and transcript prediction of sugarcane transcriptome.
a, length distribution of orf-containing transcripts resulted from transdecoder and evigen.
b, length distribution of predicted transcripts by evigene in pacbio data.
c, length distribution of predicted transcripts by evigene in de novo contig data
analysis of candidate non-coding rnas
the proportion of candidate non-coding transcripts  was different between the pacbio dataset  and de novo dataset .
due to a large number of non-coding contigs in the de novo contig dataset, which were likely from different non-coding rna classes, such as transfer rnas , small rnas, micrornas  and ribosomal rnas , and also de novo assembly artifacts; only the candidate long non-coding transcripts from the pacbio  were used for further characterization.
the pacbio non-coding transcript set had a length ranging from  <dig> to  <dig>  bp.
when compared against the ncbi nucleotide nr database, it was found that  <dig> transcripts matched sequences from bacterial, fungal and insect sources.
the remaining  <dig>  sequences included  <dig>  sequences matching the ncbi nr nucleotide database, and  <dig>  transcript isoforms that did not match any entries in nr database.
a total of  <dig>  sequences exhibited similarity to protein-encoding sequences , and  <dig>  sequences matched ncbi non-coding entries belonging to genomic sequences of the grass family.
these  <dig>  sequences contained  <dig> transcript sequences matching the predicted non-coding rnas of zea mays and setaria italica.
the final retained candidate long non-coding rnas for this dataset were  <dig>  transcript sequences , including  <dig>  non-orf, non-blast hit transcript sequences .
repeat content analysis
the repeat masking against the customized repeat library for viridiplantae showed that the total number of interspersed repeats within the pacbio transcript isoform data was  <dig>  ; including  <dig> % retroelements,  <dig> % dna transposons and  <dig> % unclassified repeat elements.
the retroelements included short interspersed nuclear elements - sines , long interspersed nuclear elements - lines  and long terminal repeat  elements .
amongst all repeat classes, the ltr gypsy/dirs <dig> was the most abundant which made up to  <dig> %, following by ltr ty1/copia , lines l1/cin <dig>  and dna transposon tourist/harbinger .
in the de novo dataset, there were  <dig>  interspersed repeats identified , including  <dig> % retroelements ,  <dig> % dna transposons and  <dig> % unclassified repeats).
gypsy/dirs <dig> , tourist/harbinger , ty1/copia  and lines l1/cin <dig>  were the dominant repeats in the de novo dataset.
all details are presented in additional file 1: table s <dig>
a total of  <dig>  ssrs were discovered in  <dig>  pacbio transcript isoforms , while a total of  <dig>  ssrs were identified in  <dig>  de novo contigs .
in both cases, the most abundant motifs detected were tri-nucleotide , followed by di-nucleotide motifs .
the ssrs detected in both pacbio transcript isoforms and de novo transcript contigs are presented in additional file 1: table s <dig>
transcript annotation
using  <dig>  pacbio transcript isoforms  , we found that there were  <dig> additional sequences from other sources, such as bacteria, fungi and insects, present in the samples.
this made up a total of  <dig> cross-contaminated sequences  in the original dataset and these were subsequently removed prior to functional annotation.
of  <dig>  remaining sequences,  <dig>   matched the ncbi nr nucleotide database, and  <dig>  sequences that did not return any matches while containing an orf, which could potentially be novel transcripts in the sugarcane transcriptome.
when compared against the viridiplantae protein database,  <dig>  transcript isoforms  showed similarity against  <dig>  plant protein sequences.
there were  <dig>  transcript isoforms  of pacbio transcript isoforms matching  <dig>  sorghum proteins,  <dig>   matching  <dig>  sucest entries  transcript isoforms remained).
the comparison between pacbio transcript isoforms and de novo transcript contigs in table  <dig>  showed that de novo contig dataset matched more viridiplantae proteins , sorghum proteins  and sucest entries .table  <dig> annotation of sugarcane transcriptome
viridiplantae proteinsa
aat an e-value ≤1e-5
there were  <dig> pacbio transcript isoforms matching the sugarcane chloroplast genome, and  <dig> matching the sorghum and maize mitochondrial genomes, while of the de novo transcript contigs,  <dig> matched the chloroplast genome and  <dig> matched the mitochondrial genome .
even though chloroplast and mitochondrial reads were removed prior to de novo assembly to reduce the read abundancy, there could be some chloroplast and mitochondrial reads still remaining in the rna-seq data for which a stringent setting and the mitochondrial genomes from closely related species were used for mapping.
using the plant transcription factor  database, we identified a total of  <dig>  tfs in the pacbio transcript isoforms, including  <dig>  similar to those in sorghum,  <dig> similar to maize tfs,  <dig> similar to rice tfs and  <dig> similar to tfs in other plant species .
there were  <dig> additional tfs identified using the grassius sugarcane tfs, and all  <dig>  identified tfs were distributed on  <dig>  tf-encoding transcript isoforms, which accounted for ~ <dig> % of the total sequences.
these tfs were from  <dig> annotated tf families, and their distribution is presented in the additional file 1: figure s <dig>  in the de novo transcript contigs, we identified  <dig>  tfs from  <dig> tf families, belonging to  <dig>  tf-encoding transcript contigs.
two families  were not found in the de novo transcript contigs, compared to the pacbio dataset.table  <dig> sugarcane transcription factors analysis of pacbio transcript isoforms
using all orf-containing transcripts in functional annotation, a total of  <dig>  go terms were assigned to  <dig>  pacbio transcript isoforms .
these go terms were classified into three main classes, cellular component, molecular function and biological process.
among the cellular component category, the highest proportion of transcript isoforms was involved in cell and cell part , organelle  and macromolecular complex .
in molecular function, binding was dominant , followed by catalytic activity , transporter , structural molecule activity  and transcription regulator activity .
in biological process, the most transcript isoforms were assigned to metabolic process , cellular process , localization , biological regulation , pigmentation , response to stimulus  and cellular component organization .
a comparison of enriched go terms between the pacbio transcript isoforms and de novo transcript contigs  is presented in fig.
<dig> fig.
<dig> gene ontology enrichment analysis of sugarcane transcript sequences.
for de novo transcript contigs, only go terms represented for  <dig>  sequences were used
kegg metabolic pathway analysis provided additional possible functional information showing the pathways that the transcript isoforms take part in, since one gene could be assigned to more than one go term in the gene ontology annotation.
the results are presented in fig.
<dig>  expressing the percentage of transcripts involved in the pathways.
a total of  <dig>  pacbio transcript isoforms  matched  <dig>  kegg pathway annotations , while  <dig>  de novo transcript contigs  matched  <dig>  ko annotations.
the largest functional pathway was metabolic pathway, representing  <dig> % and  <dig> % for pacbio and de novo transcript datasets, respectively; followed by biosynthesis of secondary metabolites , biosynthesis of antibiotics , ribosome , splicesome , biosynthesis of amino acids  and carbon metabolism .
additional file 1: figure s <dig> shows some important pathways for sugarcane including purine metabolism, starch and sucrose metabolism, phenylpropanoid biosynthesis  and carbon fixation pathway.fig.
<dig> kegg metabolic pathway classification of sugarcane pacbio transcript isoforms and de novo transcript contigs
comparative analysis with closely related species
it was found that,  <dig> % of total pacbio transcript isoforms  and 41% of de novo contigs were aligned to the sorghum genome.
when considering only retained orfs from transdecoder in both datasets,  <dig> % of pacbio orfs and 70% of de novo orfs mapped to the sorghum genome.
there were  <dig> % of evigene predicted pacbio transcripts and only 37% of the de novo predicted transcripts that aligned to the sorghum genome.
details are provided in additional file 1: table s <dig> fig.
<dig> pacbio transcript isoforms aligned against the sorghum chromosomes.
purple blocks represent for the transcript isoforms distribution along the sorghum chromosomes
discussions
in the eukaryotic cell, about 95% of genes undergo rna transcript splicing where most introns are removed and exons are retained, resulting in multiple alternative transcripts  of the gene  <cit> .
isoforms of each gene can be formed by cis-splicing or trans-splicing, where different exons are combined together to create an mrna molecule.
in general, cis-splicing involves processing a single molecule  <cit> , whereas in trans-splicing, many pre-mrnas are processed and their exons joined and ligated  <cit> .
most nuclear gene-related splicing in plants have been found to involve cis-splicing, in different modes such as intron retention, alternative splice or exon skipping/inclusion.
for instance, in the model plant, arabidopsis thaliana, it was found that about 61% of multi-exonic genes displayed alternative splicing, including different modes, ~40% intron retention, ~ <dig> % alternative 3’ splice site, ~8% exon skipping/inclusion and ~ <dig> % alternative 5’ splice site  <cit> .
similar proportions were reported in the transcriptomes of sorghum  <cit>  and maize  <cit>  with intron retention being the most abundant splicing mode, accounting for about 40%.
trans-splicing has been observed mostly in plant organellar genomes, such as in mitochondria  but was recently also found in maize nuclear genes  <cit> .
therefore, it is estimated that there are more transcripts than genes in a given genome, for example, in arabidopsis cells, there are on average  <dig>  transcripts from about  <dig>  genes  <cit> .
different spliced isoforms can be translated into different proteins and could be present in the sample at different levels of expression, at different developmental stages  <cit> .
within the total transcript population, about 20% is comprised of high abundance transcripts of a few genes , about 40-60% is from the intermediate abundance transcripts , while 20–40% is from the rare transcripts  <cit> .
we generated an initial collection of  <dig>  unique sugarcane transcript isoforms in our experiment.
this transcriptome dataset provides direct evidence of alternative splicing of transcripts for each of the genes in the sugarcane genome with higher confidence, compared to alternative splicing events reported in the de novo assembly from short-reads.
even though pacbio offers longer reads than other current platforms , it has a higher error rate  <cit> .
in iso-seq, the error rate is expected to be lower since the reads are a consensus from multiple sequencing passes of the circular cdna in the smrt cell .
however, due to relatively low supporting reads for the low quality reads in the quiver self-correction pipeline, it was observed that the transcript isoforms produced a large proportion of fragmented orfs, and the transcript prediction resulted in a low number of transcripts.
in a study in  <cit>  on sorghum transcripts, it was estimated that the error rate in pacbio iso-seq was  <dig> %, including  <dig> % mismatches,  <dig> % insertion  and  <dig> % deletion .
in this study, after a further error correction using the rna-seq reads obtained from the same experiment, the pacbio transcripts isoforms generated longer orfs, better prediction results, better performance in completeness assessments and more reads aligned to the closely related sorghum genome.
it is important to note that the rna-seq reads were generated from only internodal samples, while pacbio data included internode, leaf and root samples.
therefore, it is expected that there were low quality transcripts originally from leaf and root tissues, and also rare transcripts resulting from the normalization, left un-corrected in this second error correction.
using the pacbio iso-seq to capture full-length transcripts without assembly overcomes the difficulty posed by the short-read data.
the comparative analysis with the de novo assembled contigs from illumina rna-seq reads, allowed us to evaluate the benefits of each of the assemblies in constructing the sugarcane transcriptome.
the short-reads from the illumina platform have been used widely for rna-seq differential gene expression analysis  <cit>  since it provides sufficient depth and a lower error rate compared to reads generated from pacbio.
however, due to the complexity of the alternative spicing mechanism of eukaryotic cells, recovering full-length transcripts has been a challenge for most of the assemblers using short reads, such as trinity, soapdenovo-trans or velvet/oases.
many more transcripts were generated from the de novo assembly in this study compared to the pacbio transcript isoforms, as well as the expected number of transcripts.
this was in agreement with most de novo assembly studies, such as in  <cit> ,  <cit>  and  <cit> .
it was found that de novo assembly from combined multiple settings/assemblers showed that this analysis represented well the sample from which it was derived, with  <dig> % of reads mapping back to transcripts, compared to  <dig> % to pacbio transcript isoforms.
the average proportion mapping to the reference transcriptome found in the literature is around 70-90%, i.e.
in  <cit> , since there is a proportion of reads from the lowly expressed transcripts  that are not assembled into contigs.
the higher read mapping rate in our de novo assembly could also be attributed to our library preparation, in which 150 bp paired reads were generated from a library of average 200 bp fragments, creating overlapped reads easy to assemble .
it could also be due to the great depth of reads used for assembly .
even though the pacbio data included the same internodal rnas as de novo dataset, it has been through different library preparation, where the cdnas were produced from only polya rnas .
the comparison of transcripts between the two assemblies suggests that the common transcripts between the two assemblies would mostly be polya rnas; while the transcripts unique to pacbio  could be rare transcripts that come from the normalization process and wider tissue coverage, and those unique to de novo assembly  could be attributed to other types of rnas in the samples.
it has been proposed that de novo transcript contigs could be a good resource for studying the diversity of non-coding rnas  <cit> .
the higher cegma/busco completeness alignment of de novo transcript contigs  compared to pacbio transcript isoforms  suggests that this dataset contained more expected core conserved genes, and indirectly indicates that more genes were captured.
this result was also consistent with the blast search, in which more viridiplantae proteins matched the de novo transcript contigs  than that of the pacbio transcript isoforms  at a length coverage cutoff of 70%.
the de novo transcript contigs incorporated more gene content  compared to the pacbio transcript isoforms.
this could be due to the high coverage of rna-seq sequencing and the use of multiple settings/assemblers in our de novo assembly, while the sequencing depth in pacbio iso-seq was still relatively low.
the pacbio iso-seq, on the other hand, was shown to be better in recovering full-length transcript isoforms , to include more coding transcripts , and to have a much longer orfs.
even though the number of predicted orfs in the pacbio dataset was less than that in the de novo dataset, the orfs had a longer n <dig>  while de novo orfs’ n <dig> was only 738 bp.
similarly, there were less predicted transcripts in the pacbio dataset than that in the de novo contig set.
the number of predicted main transcripts  in the pacbio dataset was  <dig> , which was approximately 71% of the total predicted genes in sugarcane .
combined with the cegma/busco alignment, this suggests that a proportion of the genes were missing in the pacbio data.
many of these could be genes expressed in different tissues or developmental stages to those sampled here.
in the de novo transcript contigs,  <dig>  predicted main transcripts were obtained, which was about ~162% of the predicted sugarcane genes.
it could be that in de novo assembly, not all transcripts were recovered in full-length.
there could be genes that were represented by several different contigs resulting in a total predicted transcript number greater than the true number expected.
this result compares with the unigenes dataset  <cit> , which had  <dig> – <dig> % cegma/busco completeness, but resulted in only  <dig>  main predicted transcripts.
the unigene dataset originally had  <dig>  contigs, assembled from ~ <dig> million of 72 bp paired-end reads, of which only unigenes  were retained for further analysis.
the high cegma/busco alignment could be due to good representation of contigs for the samples, while a low number of predicted transcripts could have resulted from only the main isoforms being retained in the final contig set.
analysis of the predicted transcripts in the sogi dataset concluded that was in agreement with the estimation reported in  <cit>  and indicates that this dataset represents ~90% of predicted sugarcane genes.
of a total  <dig>  predicted transcripts by evigene,  <dig>  main transcripts were identified, which was equivalent to  <dig> % of the total sugarcane predicted genes, and close to the figure above.
this dataset had a low cegma/busco alignment, ranging from  <dig>  to  <dig> %.
it could be that the cegma/busco alignment required 70% alignment to the conserved proteins, whereas sogi database contained a proportion of short est sequences  making the alignment length less than the threshold used by the programs, and therefore resulted in low completeness.
amongst all dataset, pacbio transcript isoforms had the longest average length of the largest  <dig>  proteins .
even though the de novo assembly was shown to have better completion, suggesting more gene content included, it represented fragmented sequences.
the number of alternate transcripts reported for pacbio data was  <dig> , and for the de novo dataset was  <dig> , despite many more de novo input transcript contigs in this dataset.
the pacbio predicted transcripts could be used to improve the length of sugarcane predicted gene models, and sugarcane protein sequences, which are covered by the sogi database but are not full-length sequences.
long non-coding rnas are rnas longer than 300 bp that do not encode proteins  and potentially play important roles in gene regulation of eukaryotic cells  <cit> .
their numbers, characteristics and genetic patterns in the genome still remain unclear.
the prediction and annotation of long non-coding rnas is normally challenging since unlike the coding rnas, they are not orthologous and lack homology between closely related species.
therefore, the information from one species is not useful in non-coding rna prediction for other species  <cit> .
more often, long non-coding rnas of a given genome are predicted by subjecting the un-known non-orf-containing rnas to a model, which is built on a set of high confidence non-coding rnas and a partition of coding transcripts of that genome  <cit> .
in this study, de novo assembly included more non-coding rnas  compared to the pacbio dataset .
we identified  <dig>  transcript sequences  as putative long-noncoding rnas in the pacbio transcript isoforms.
this was done by comparing the non-coding transcripts against the available protein and genomic databases.
CONCLUSIONS
the transcript data generated in this study probably accounts for about 71% of the total predicted genes in the sugarcane genome.
the pacbio iso-seq analysis recovered more full-length transcripts, with a longer n <dig>  more orfs and predicted transcripts and higher average length of the largest  <dig>  proteins, compared to that of the de novo contigs from rna-seq.
analysis of the gene content in the two assemblies suggests that rna-seq covered more gene content, and more rna classes, probably as a result of the greater sequencing depth.
the majority of transcript isoforms captured in pacbio iso-seq were protein-coding sequences , whereas only  <dig> % of the total rna-seq de novo contigs contained orfs.
about  <dig> % and  <dig> % of pacbio and de novo transcripts matched the viridiplantae protein database, respectively.
the use of normalization and the inclusion of more tissue  in the iso-seq library preparation may have contributed to the recovery of the unique fraction  attributed to rare and tissue-specific transcripts that were not covered in the rna-seq.
comparative analysis with the sorghum genome indicated a high content of orthologous genes between the two genomes.
the total set of  <dig>  predicted transcripts in the study could be used to improve sugarcane gene models and sugarcane proteins in the sugarcane databases like sogi that lack full-length sequences.
this dataset will serve as reference sequences representing full-length transcript isoforms that are expressed in sugarcane leaf, internode and root tissues, and facilitate differential expression analysis allowing exploration of different isoforms of genes to be studied.
this reference database is termed as the sugit database .

 understanding of the sugarcane transcriptome is limited due to the complexity in gene copy number, repetitive content, and heterozygosity in the genome. it is not clear how many transcript isoforms result from the alternative splicing in this potentially very complex transcriptome. sugarcane is a polyploid hybrid between saccharum officinarum and s. spontaneum, and each sugarcane hybrid has its own unique chromosome set, containing up to copies of each gene and a total ~ predicted genes. therefore, it is expected that sugarcane transcripts represent transcription of genes/homoelogues that are not only unique to the progenitor genomes but also transcription from alternate splicing, which we collectively refer to as transcript isoforms. sugarcane genes such as those in the sucrose phosphate synthase gene families, invertase genes and sucrose synthase family have been shown to be comprised of many isoforms. most of the sugarcane studies, including transcriptome studies found in the literature, i.e. in and, have been based on sorghum genomic/transcript sequences which have the highest gene synteny and orthologous alignment with the sugarcane genome; sugarcane expressed sequence tags; saccharum officinarum gene indices - sogi v representing ~90% of the estimated genes in s. officinarum; and other resources reviewed in. studies based on these databases have provided useful information on the sugarcane transcriptome, while a whole genome sequence is not yet available. however, it is thought that there are still many sugarcane genes missing in these databases and in addition, the full-length sequences of distinct transcript isoforms are not included. use of these transcript databases for rna-seq analysis leading to the identification differentially expressed genes does not provide information on the corresponding isoform/s or the homoelogue/s contributing to the differential expression. there is a need to construct fl transcript sequences including such isoforms to facilitate analysis of isoform differential expression, and also to extend our understanding of the sugarcane transcriptome. the transcriptome poses a great challenge when it comes to assembly and annotation. the differences in transcript abundance and the presence of different isoforms, greatly challenge the assembly of a transcriptome from short-reads; since the assemblers cannot distinguish between reads originally from different transcripts/isoforms carrying the same exons. to date, most sequencing platforms offer a read-length which is shorter than the typical length of a eukaryotic mrna. the transcriptome sequences obtained from second generation sequencing technology have been playing an important role in capturing the diversity in the rna populations at a greater sequencing depth. however, a precise prediction and identification of the alternative transcript splicing has not been possible. algorithms in transcript splice-aware assemblers have been developed to detect splicing junctions and recover transcript isoforms by using information from short-reads, but these have not always been confirmed. that is, quite often these approaches overestimate and report spurious computational isoforms rather than picking up only biological ones. overall, the assemblies from short-read data normally end up identifying more transcripts than expected, which may be attributed not only to the diversity of rnas and diversity of transcript/isoforms in the transcriptome, but also to the limitation in recovering fl transcripts. most studies use these tools, then filter the transcripts through clustering by retaining the longest sequence in each cluster as representative for analysis, and consider them to be the major isoforms or unigenes. current algorithms such as in bernard et al. for isoform identification and quantification require longer reads and cannot tackle genes with too many exons. with the advent of third generation sequencing technology, the cost-per-transcriptome has been reduced, whereas the length of the sequencing reads has been increased significantly. as of august, the average read length of pacbio single molecule, real-time sequencing is >10 kb and real length can be up to 60 kb. this protocol has been applied in some recent studies, for example, detecting alternative splicing events in unique transcripts, covering ~89% of the total sorghum annotated genes; and producing unique transcripts in the maize transcriptome derived from six different tissues, covering ~70% of the annotated genes. this study represents the first full-length transcriptome reference sequences from sugarcane derived from three different tissues, of different developmental stages, by using the pacbio long-read iso-seq technique. in addition, rna-seq was used to improve the pacbio transcript isoforms by short-read error correction, and comparison between sugarcane transcripts obtained/assembled from these two different platforms. annotation of the sugarcane fl transcript isoforms could improve sugarcane genome models, contribute towards understanding of the complexity of the sugarcane genome and serve as reference sequences for differential expression analysis in the future.  sugarcane transcriptome from pacbio isoform sequencing a pooled sample representing polya rnas from three tissues, of different developmental stages was sequenced to obtain a wide coverage of the sugarcane transcriptome. a total of reads of inserts was generated, with a total of, nucleotides from six smrt cells of non-normalized bins and normalized bins, including fl non-chimeric rois and non-fl rois. the length of rois ranged from 300 bp to  kb, with an n of  bp. the length distribution of all roi data is presented in additional file 1: figure s a total of high quality sequences and low quality sequences were obtained from quiver polishing, referred to as polished transcript isoforms. the total unique, non-redundant transcript isoforms included sequences, with the length ranging from 301 bp to  bp, n of  bp and n of  bp and % gc content. improving pacbio transcript quality by error correction using rna-seq reads we followed two error correction pipelines, using three datasets from rna-seq derived from the same experiment;, non-normalized reads,, reads bbnorm-normalized and, trinity-normalized reads. overall, the error correction led to improvement in transcript prediction, more transcripts covered the full-length of known proteins, longer open reading frames, better completeness results in cegma/busco assessments, and higher alignment rate of transcript isoforms to the sorghum genome. this only resulted in a slight change in the total number of transcripts isoforms after removing all exact 100% identical sequences. the lordec error correction outperformed proovread with our transcript data. the best corrected set was derived from lordec using the trinity-normalized reads, which resulted in % of the total transcripts with orfs passing the evigene score in transcript prediction, while that of the non-corrected transcripts was only %. there were orfs detected by transdecoder in this corrected pacbio transcript isoform set, compared to orfs for the non-corrected dataset. the retained orfs with pfam and viridiplantae protein hits from the corrected dataset had an n of  bp, while that of the non-corrected was 708 bp. the cegma and busco alignments showed that the corrected pacbio dataset had a higher completeness level than the non-corrected. about % of the corrected transcripts aligned to the sorghum genome, while the alignment of the non-corrected transcript isoforms was %. we chose the pacbio dataset corrected by lordec, using reads normalized by trinity, for downstream analysis. this final pacbio transcript isoforms set had sequences after removal of strictly identical sequences; of total length ~193 mb, with individual transcript length ranging from to  bp, n of, n of, and % gc content.table summary of correction of pacbio transcript isoform data using illumina short-reads de novo assembly of the sugarcane transcriptome from short-reads the assembly of the sugarcane transcriptome from illumina rna-seq short-reads was carried out to provide a comparative reference for the transcript isoform sequences obtained from pacbio iso-seq, since rna-seq has been utilized widely in construction of transcriptomes. of million total raw reads generated,, reads survived after trimming, having a quality score cutoff of q trinity normalization retained only 6% at a maximum coverage of while retaining 15% and 21% at maximum coverage of and, respectively. a qc report of all read datasets was generated by fastqc v.. the use of bbnorm was selected due to the fact that trinity normalization heavily reduced the reads compared to bbnorm at the same maximum coverage cutoff. about 37% of the total reads remained at maximum coverage by bbnorm package for de novo assembly. after assembly and individual clustering, four initial assemblies were obtained from trinity, clc-gwb, velvet/oases, and soapdenovo-trans, respectively. we observed varied contig number, n cumulative length and length distribution in each of the assemblies. the total number of contigs from the trinity-assembly was, while that from the clc-gwb assembly was, velvet/oases assembly gave and soaptrans assembly. table  summarizes the de novo assembly results in this study, including all major statistics and qc, more details see additional file 1: figure s the final total number of contigs clustered by cd-hit-est at 95% identity and after retaining transcripts with length from 300 bp to 10 kb, was, of ~967 mb, having an n of  bp, n of 745 bp and % gc. this clustered assembly was referred to as the de novo transcript contigs.table comparison of de novo assemblies used in this study when compared to pacbio transcript isoforms by blastn, % of the de novo transcript contigs exhibited similarity to % of the pacbio transcript isoforms. there were % of de novo transcript contigs and % of pacbio transcript isoforms that were unique to each of the datasets.fig. comparison between the sugarcane pacbio transcript isoforms and de novo transcript contigs analysis of reads mapping back to transcripts in mapping back the rna-seq read data to transcripts from both pacbio and de novo transcripts, we observed % of reads mapped back to pacbio transcript isoforms, while % mapped to the de novo transcript contigs. figure  shows the average coverage plotted against the transcript length for both assemblies.fig. average coverage of sugarcane de novo contigs and pacbio isoforms obtained from read mapping. a, coverage of de novo transcript contigs. b, coverage of pacbio transcript isoforms transcriptome completeness analysis in both cegma and busco alignments, the pacbio dataset showed a lower completeness level than the de novo dataset. the pacbio transcripts had % cegma alignment, and the de novo assembly had % cegma alignment. there were no missing cegma in the de novo transcript contigs, and there was % missing cegmas in the pacbio dataset. similarly, in busco alignment to conserved proteins, the pacbio transcript isoforms had lower completeness than that of de novo assembly, by having % completeness compared to %. the de novo assembly had a higher level of duplication in the busco alignment, suggesting that the assembly contains duplicate contigs of different lengths assembled by different assemblers and retained after clustering. using the same cegma and busco protein alignments, we assessed the unigene dataset from and the sogi database to determine the consistency of the methods. the unigenes had % cegma completeness, % busco completeness; and the sogi dataset had % cegma completeness and % busco alignment completeness. the sogi dataset had the largest proportion of partial/fragmented and missing proteins, since this dataset contained gene indices and ests.table transcriptome coverage analysis based upon cegma and busco alignment counting the full-length and nearly full-length transcripts the number of transcripts appearing to be full-length with at least 90% and 70% alignment coverage of the viridiplantae uniprot proteins was estimated and compared between pacbio transcript isoforms and de novo transcript contigs. the pacbio dataset had transcripts appearing to be full-length, and transcripts, and that of the de novo transcript contigs were and at 90% and 70%, respectively. analysis of the matched proteins at the cutoff of 70% alignment coverage from both assemblies indicated that only proteins were common between these two assemblies, with unique to pacbio dataset, and unique to the de novo transcript contigs.fig. full-length analysis between sugarcane pacbio transcript isoforms and de novo transcript contigs. a, counts of proteins covered by transcripts at different thresholds. b, comparison between the protein hits from pacbio and de novo transcripts which covered at least 70% viridiplantae protein length as these results considered each protein from the uniprot database as only one count, regardless of the presence of different isoforms that carry the same protein sequence part getting matched many times. a modified approach to counting the full-length transcripts for isoforms was applied, in which we took all counts of isoforms that hit the same protein into consideration and estimated the number of full-length proteins covered. using this strategy, we found in the pacbio data, transcripts that covered ≥90% of viridiplantae proteins and transcripts that covered ≥70% of viridiplantae proteins. in the de novo transcript contigs, it was and protein hits covered by ≥90% and ≥70%, respectively. de novo transcript contigs had more proteins covered at lower percentage due to the greater duplication retained in the assembly, and inclusion of more partial gene content. when protein hits of ≥90% alignment coverage from the two results were compared, the unique protein hits of pacbio was and de novo was, which were the same as in the first approach. investigation of full-length genes from sugarcane and other grass species showed that pacbio dataset resulted in a better performance in terms of recovering the full-length sequence of these specific genes. at an e-value ≤1e- there were more genes covered by transcripts at 90% and 70% in pacbio transcript isoforms than that in de novo transcript contigs, unigene set and sogi database, respectively. the lower full-length gene count in unigenes could be due to only main isoforms being retained for this dataset, while lower full-length gene count in sogi database could be due to the fraction of ests in it.table alignment and full-length assessment of a selected gene set evidence of alternative splicing in the sugarcane transcriptome from pacbio iso-seq using the results from pacbio transcript isoforms mapped against the sorghum genome, our in-house sugarcane whole genome assembly from sugarcane cultivar q and particular contigs that spanned through the sucrose phosphate synthase a and cellulase genes, we were able to visualize alternative splicing in the pacbio transcript isoforms. using the transcriptome analysis pipeline for isoform sequencing described in, we detected amongst those transcript isoforms aligned against the sorghum genome, alternative splicing events, including intron retention, skipped exon, alternative 5’splice sites and alternative 3’splice sites. an estimation of exons per transcripts amongst the transcript isoforms aligned against the sorghum genome is also included.fig. evidence of different transcript isoforms of sugarcane transcriptome present in the pacbio transcript dataset. a, isoforms aligned against the sorghum chromosome b, isoforms aligned to contigs of our in-house sugarcane whole genome de novo assembly. c, different transcript isoforms aligned to sucrose phosphate synthase gene and cellulase gene contigs. d, average exons per transcript estimated based on the transcript isoforms aligned against sorghum genome prediction of potential coding regions, main and alternate transcripts analysis we analyzed the candidate coding regions in the transcript sequences by retaining only open reading frames that exhibited homology with the pfam protein domain database or the uniprot viridiplantae known protein database, which were more likely to be biologically real. there were non-overlapped orfs detected in the pacbio transcript isoform sequences, belonged to orf-containing transcript isoforms. only isoforms did not contain orfs and these were used for characterization of non-coding rnas, in the next section. of the total orf-containing transcripts, matched the viridiplantae proteins, while matched the pfam protein domains. the total number of transcript sequences retained in combination with the transdecoder frame-score was, with lengths ranging from 300 bp to  bp, and n of  bp. there were more sequences retained in the de novo assembly, since it started with more data, however transcript contig length was shorter than that of the pacbio dataset. a total orfs were found in orf-containing contigs. when combined with the homology search results, transcript contigs matched the viridiplantae protein database, and contigs matched the pfam protein domains. the final number of contigs retained by transdecoder was, accounting for % of the total de novo assembly. this final contig set had lengths ranging from 300 bp to  bp and n of 738 bp. when predicting the potential coding genes using the evigene pipeline, the total number of predicted transcripts in the pacbio transcript isoforms was, while the dropset had transcripts. the average length of the largest proteins from the dataset was reported to be aa, and all candidate transcripts had an n of  bp and cds length ranged from 186 bp to  bp. as for the prediction by transdecoder, the de novo contig set had more orfs, and therefore, more predicted transcripts both main and alternate. there were predicted transcripts with an n of compared to pacbio transcript isoforms, the de novo predicted transcripts had much shorter length distribution and average length of the largest proteins. using the same prediction approach on the unigenes and sogi, we found predicted main transcripts in the unigenes, and predicted transcripts in sogi. the average length of the largest proteins for the unigenes was aa, while that of sogi was aa. all the results of orf and transcript prediction are presented in table  and fig.  table open reading frame and transcript prediction analysis of sugarcane transcriptome sequence data retained transcriptsa ave length 1 k proteinsb atranscripts with pfam and viridiplantae hits. baverage length of the largest proteins fig. analysis of orfs and transcript prediction of sugarcane transcriptome. a, length distribution of orf-containing transcripts resulted from transdecoder and evigen. b, length distribution of predicted transcripts by evigene in pacbio data. c, length distribution of predicted transcripts by evigene in de novo contig data analysis of candidate non-coding rnas the proportion of candidate non-coding transcripts was different between the pacbio dataset and de novo dataset. due to a large number of non-coding contigs in the de novo contig dataset, which were likely from different non-coding rna classes, such as transfer rnas, small rnas, micrornas and ribosomal rnas, and also de novo assembly artifacts; only the candidate long non-coding transcripts from the pacbio were used for further characterization. the pacbio non-coding transcript set had a length ranging from to  bp. when compared against the ncbi nucleotide nr database, it was found that transcripts matched sequences from bacterial, fungal and insect sources. the remaining sequences included sequences matching the ncbi nr nucleotide database, and transcript isoforms that did not match any entries in nr database. a total of sequences exhibited similarity to protein-encoding sequences, and sequences matched ncbi non-coding entries belonging to genomic sequences of the grass family. these sequences contained transcript sequences matching the predicted non-coding rnas of zea mays and setaria italica. the final retained candidate long non-coding rnas for this dataset were transcript sequences, including non-orf, non-blast hit transcript sequences. repeat content analysis the repeat masking against the customized repeat library for viridiplantae showed that the total number of interspersed repeats within the pacbio transcript isoform data was; including % retroelements, % dna transposons and % unclassified repeat elements. the retroelements included short interspersed nuclear elements - sines, long interspersed nuclear elements - lines and long terminal repeat elements. amongst all repeat classes, the ltr gypsy/dirs was the most abundant which made up to %, following by ltr ty1/copia, lines l1/cin and dna transposon tourist/harbinger. in the de novo dataset, there were interspersed repeats identified, including % retroelements, % dna transposons and % unclassified repeats). gypsy/dirs, tourist/harbinger, ty1/copia and lines l1/cin were the dominant repeats in the de novo dataset. all details are presented in additional file 1: table s a total of ssrs were discovered in pacbio transcript isoforms, while a total of ssrs were identified in de novo contigs. in both cases, the most abundant motifs detected were tri-nucleotide, followed by di-nucleotide motifs. the ssrs detected in both pacbio transcript isoforms and de novo transcript contigs are presented in additional file 1: table s transcript annotation using pacbio transcript isoforms, we found that there were additional sequences from other sources, such as bacteria, fungi and insects, present in the samples. this made up a total of cross-contaminated sequences in the original dataset and these were subsequently removed prior to functional annotation. of remaining sequences, matched the ncbi nr nucleotide database, and sequences that did not return any matches while containing an orf, which could potentially be novel transcripts in the sugarcane transcriptome. when compared against the viridiplantae protein database, transcript isoforms showed similarity against plant protein sequences. there were transcript isoforms of pacbio transcript isoforms matching sorghum proteins, matching sucest entries transcript isoforms remained). the comparison between pacbio transcript isoforms and de novo transcript contigs in table  showed that de novo contig dataset matched more viridiplantae proteins, sorghum proteins and sucest entries.table annotation of sugarcane transcriptome viridiplantae proteinsa aat an e-value ≤1e-5 there were pacbio transcript isoforms matching the sugarcane chloroplast genome, and matching the sorghum and maize mitochondrial genomes, while of the de novo transcript contigs, matched the chloroplast genome and matched the mitochondrial genome. even though chloroplast and mitochondrial reads were removed prior to de novo assembly to reduce the read abundancy, there could be some chloroplast and mitochondrial reads still remaining in the rna-seq data for which a stringent setting and the mitochondrial genomes from closely related species were used for mapping. using the plant transcription factor database, we identified a total of tfs in the pacbio transcript isoforms, including similar to those in sorghum, similar to maize tfs, similar to rice tfs and similar to tfs in other plant species. there were additional tfs identified using the grassius sugarcane tfs, and all identified tfs were distributed on tf-encoding transcript isoforms, which accounted for ~ % of the total sequences. these tfs were from annotated tf families, and their distribution is presented in the additional file 1: figure s in the de novo transcript contigs, we identified tfs from tf families, belonging to tf-encoding transcript contigs. two families were not found in the de novo transcript contigs, compared to the pacbio dataset.table sugarcane transcription factors analysis of pacbio transcript isoforms using all orf-containing transcripts in functional annotation, a total of go terms were assigned to pacbio transcript isoforms. these go terms were classified into three main classes, cellular component, molecular function and biological process. among the cellular component category, the highest proportion of transcript isoforms was involved in cell and cell part, organelle and macromolecular complex. in molecular function, binding was dominant, followed by catalytic activity, transporter, structural molecule activity and transcription regulator activity. in biological process, the most transcript isoforms were assigned to metabolic process, cellular process, localization, biological regulation, pigmentation, response to stimulus and cellular component organization. a comparison of enriched go terms between the pacbio transcript isoforms and de novo transcript contigs is presented in fig.  fig. gene ontology enrichment analysis of sugarcane transcript sequences. for de novo transcript contigs, only go terms represented for sequences were used kegg metabolic pathway analysis provided additional possible functional information showing the pathways that the transcript isoforms take part in, since one gene could be assigned to more than one go term in the gene ontology annotation. the results are presented in fig.  expressing the percentage of transcripts involved in the pathways. a total of pacbio transcript isoforms matched kegg pathway annotations, while de novo transcript contigs matched ko annotations. the largest functional pathway was metabolic pathway, representing % and % for pacbio and de novo transcript datasets, respectively; followed by biosynthesis of secondary metabolites, biosynthesis of antibiotics, ribosome, splicesome, biosynthesis of amino acids and carbon metabolism. additional file 1: figure s shows some important pathways for sugarcane including purine metabolism, starch and sucrose metabolism, phenylpropanoid biosynthesis and carbon fixation pathway.fig. kegg metabolic pathway classification of sugarcane pacbio transcript isoforms and de novo transcript contigs comparative analysis with closely related species it was found that, % of total pacbio transcript isoforms and 41% of de novo contigs were aligned to the sorghum genome. when considering only retained orfs from transdecoder in both datasets, % of pacbio orfs and 70% of de novo orfs mapped to the sorghum genome. there were % of evigene predicted pacbio transcripts and only 37% of the de novo predicted transcripts that aligned to the sorghum genome. details are provided in additional file 1: table s fig. pacbio transcript isoforms aligned against the sorghum chromosomes. purple blocks represent for the transcript isoforms distribution along the sorghum chromosomes discussions in the eukaryotic cell, about 95% of genes undergo rna transcript splicing where most introns are removed and exons are retained, resulting in multiple alternative transcripts of the gene. isoforms of each gene can be formed by cis-splicing or trans-splicing, where different exons are combined together to create an mrna molecule. in general, cis-splicing involves processing a single molecule, whereas in trans-splicing, many pre-mrnas are processed and their exons joined and ligated. most nuclear gene-related splicing in plants have been found to involve cis-splicing, in different modes such as intron retention, alternative splice or exon skipping/inclusion. for instance, in the model plant, arabidopsis thaliana, it was found that about 61% of multi-exonic genes displayed alternative splicing, including different modes, ~40% intron retention, ~ % alternative 3’ splice site, ~8% exon skipping/inclusion and ~ % alternative 5’ splice site. similar proportions were reported in the transcriptomes of sorghum and maize with intron retention being the most abundant splicing mode, accounting for about 40%. trans-splicing has been observed mostly in plant organellar genomes, such as in mitochondria but was recently also found in maize nuclear genes. therefore, it is estimated that there are more transcripts than genes in a given genome, for example, in arabidopsis cells, there are on average transcripts from about genes. different spliced isoforms can be translated into different proteins and could be present in the sample at different levels of expression, at different developmental stages. within the total transcript population, about 20% is comprised of high abundance transcripts of a few genes, about 40-60% is from the intermediate abundance transcripts, while 20–40% is from the rare transcripts. we generated an initial collection of unique sugarcane transcript isoforms in our experiment. this transcriptome dataset provides direct evidence of alternative splicing of transcripts for each of the genes in the sugarcane genome with higher confidence, compared to alternative splicing events reported in the de novo assembly from short-reads. even though pacbio offers longer reads than other current platforms, it has a higher error rate. in iso-seq, the error rate is expected to be lower since the reads are a consensus from multiple sequencing passes of the circular cdna in the smrt cell. however, due to relatively low supporting reads for the low quality reads in the quiver self-correction pipeline, it was observed that the transcript isoforms produced a large proportion of fragmented orfs, and the transcript prediction resulted in a low number of transcripts. in a study in on sorghum transcripts, it was estimated that the error rate in pacbio iso-seq was %, including % mismatches, % insertion and % deletion. in this study, after a further error correction using the rna-seq reads obtained from the same experiment, the pacbio transcripts isoforms generated longer orfs, better prediction results, better performance in completeness assessments and more reads aligned to the closely related sorghum genome. it is important to note that the rna-seq reads were generated from only internodal samples, while pacbio data included internode, leaf and root samples. therefore, it is expected that there were low quality transcripts originally from leaf and root tissues, and also rare transcripts resulting from the normalization, left un-corrected in this second error correction. using the pacbio iso-seq to capture full-length transcripts without assembly overcomes the difficulty posed by the short-read data. the comparative analysis with the de novo assembled contigs from illumina rna-seq reads, allowed us to evaluate the benefits of each of the assemblies in constructing the sugarcane transcriptome. the short-reads from the illumina platform have been used widely for rna-seq differential gene expression analysis since it provides sufficient depth and a lower error rate compared to reads generated from pacbio. however, due to the complexity of the alternative spicing mechanism of eukaryotic cells, recovering full-length transcripts has been a challenge for most of the assemblers using short reads, such as trinity, soapdenovo-trans or velvet/oases. many more transcripts were generated from the de novo assembly in this study compared to the pacbio transcript isoforms, as well as the expected number of transcripts. this was in agreement with most de novo assembly studies, such as in, and. it was found that de novo assembly from combined multiple settings/assemblers showed that this analysis represented well the sample from which it was derived, with % of reads mapping back to transcripts, compared to % to pacbio transcript isoforms. the average proportion mapping to the reference transcriptome found in the literature is around 70-90%, i.e. in, since there is a proportion of reads from the lowly expressed transcripts that are not assembled into contigs. the higher read mapping rate in our de novo assembly could also be attributed to our library preparation, in which 150 bp paired reads were generated from a library of average 200 bp fragments, creating overlapped reads easy to assemble. it could also be due to the great depth of reads used for assembly. even though the pacbio data included the same internodal rnas as de novo dataset, it has been through different library preparation, where the cdnas were produced from only polya rnas. the comparison of transcripts between the two assemblies suggests that the common transcripts between the two assemblies would mostly be polya rnas; while the transcripts unique to pacbio could be rare transcripts that come from the normalization process and wider tissue coverage, and those unique to de novo assembly could be attributed to other types of rnas in the samples. it has been proposed that de novo transcript contigs could be a good resource for studying the diversity of non-coding rnas. the higher cegma/busco completeness alignment of de novo transcript contigs compared to pacbio transcript isoforms suggests that this dataset contained more expected core conserved genes, and indirectly indicates that more genes were captured. this result was also consistent with the blast search, in which more viridiplantae proteins matched the de novo transcript contigs than that of the pacbio transcript isoforms at a length coverage cutoff of 70%. the de novo transcript contigs incorporated more gene content compared to the pacbio transcript isoforms. this could be due to the high coverage of rna-seq sequencing and the use of multiple settings/assemblers in our de novo assembly, while the sequencing depth in pacbio iso-seq was still relatively low. the pacbio iso-seq, on the other hand, was shown to be better in recovering full-length transcript isoforms, to include more coding transcripts, and to have a much longer orfs. even though the number of predicted orfs in the pacbio dataset was less than that in the de novo dataset, the orfs had a longer n while de novo orfs’ n was only 738 bp. similarly, there were less predicted transcripts in the pacbio dataset than that in the de novo contig set. the number of predicted main transcripts in the pacbio dataset was, which was approximately 71% of the total predicted genes in sugarcane. combined with the cegma/busco alignment, this suggests that a proportion of the genes were missing in the pacbio data. many of these could be genes expressed in different tissues or developmental stages to those sampled here. in the de novo transcript contigs, predicted main transcripts were obtained, which was about ~162% of the predicted sugarcane genes. it could be that in de novo assembly, not all transcripts were recovered in full-length. there could be genes that were represented by several different contigs resulting in a total predicted transcript number greater than the true number expected. this result compares with the unigenes dataset, which had – % cegma/busco completeness, but resulted in only main predicted transcripts. the unigene dataset originally had contigs, assembled from ~ million of 72 bp paired-end reads, of which only unigenes were retained for further analysis. the high cegma/busco alignment could be due to good representation of contigs for the samples, while a low number of predicted transcripts could have resulted from only the main isoforms being retained in the final contig set. analysis of the predicted transcripts in the sogi dataset concluded that was in agreement with the estimation reported in and indicates that this dataset represents ~90% of predicted sugarcane genes. of a total predicted transcripts by evigene, main transcripts were identified, which was equivalent to % of the total sugarcane predicted genes, and close to the figure above. this dataset had a low cegma/busco alignment, ranging from to %. it could be that the cegma/busco alignment required 70% alignment to the conserved proteins, whereas sogi database contained a proportion of short est sequences making the alignment length less than the threshold used by the programs, and therefore resulted in low completeness. amongst all dataset, pacbio transcript isoforms had the longest average length of the largest proteins. even though the de novo assembly was shown to have better completion, suggesting more gene content included, it represented fragmented sequences. the number of alternate transcripts reported for pacbio data was, and for the de novo dataset was, despite many more de novo input transcript contigs in this dataset. the pacbio predicted transcripts could be used to improve the length of sugarcane predicted gene models, and sugarcane protein sequences, which are covered by the sogi database but are not full-length sequences. long non-coding rnas are rnas longer than 300 bp that do not encode proteins and potentially play important roles in gene regulation of eukaryotic cells. their numbers, characteristics and genetic patterns in the genome still remain unclear. the prediction and annotation of long non-coding rnas is normally challenging since unlike the coding rnas, they are not orthologous and lack homology between closely related species. therefore, the information from one species is not useful in non-coding rna prediction for other species. more often, long non-coding rnas of a given genome are predicted by subjecting the un-known non-orf-containing rnas to a model, which is built on a set of high confidence non-coding rnas and a partition of coding transcripts of that genome. in this study, de novo assembly included more non-coding rnas compared to the pacbio dataset. we identified transcript sequences as putative long-noncoding rnas in the pacbio transcript isoforms. this was done by comparing the non-coding transcripts against the available protein and genomic databases.  the transcript data generated in this study probably accounts for about 71% of the total predicted genes in the sugarcane genome. the pacbio iso-seq analysis recovered more full-length transcripts, with a longer n more orfs and predicted transcripts and higher average length of the largest proteins, compared to that of the de novo contigs from rna-seq. analysis of the gene content in the two assemblies suggests that rna-seq covered more gene content, and more rna classes, probably as a result of the greater sequencing depth. the majority of transcript isoforms captured in pacbio iso-seq were protein-coding sequences, whereas only % of the total rna-seq de novo contigs contained orfs. about % and % of pacbio and de novo transcripts matched the viridiplantae protein database, respectively. the use of normalization and the inclusion of more tissue in the iso-seq library preparation may have contributed to the recovery of the unique fraction attributed to rare and tissue-specific transcripts that were not covered in the rna-seq. comparative analysis with the sorghum genome indicated a high content of orthologous genes between the two genomes. the total set of predicted transcripts in the study could be used to improve sugarcane gene models and sugarcane proteins in the sugarcane databases like sogi that lack full-length sequences. this dataset will serve as reference sequences representing full-length transcript isoforms that are expressed in sugarcane leaf, internode and root tissues, and facilitate differential expression analysis allowing exploration of different isoforms of genes to be studied. this reference database is termed as the sugit database. 
the sugarcane transcriptome was sequenced using pacbio isoform sequencing of a pooled rna sample derived from leaf, internode and root tissues, of different developmental stages, from varieties, to explore the potential for capturing full-length transcript isoforms. only 62% of pacbio transcript isoforms matched 67% of de novo contigs, while the non-matched proportions were attributed to the inclusion of leaf/root tissues and the normalization in pacbio, and the representation of more gene content and rna classes in the de novo assembly, respectively. sugarcanepolyploid transcriptometranscriptome assemblyde novo assemblyisoform sequencinghybrid assemblysugit databasehttp://dx.doi.org/ /501100001794university of queenslandhttp://dx.doi.org/ /501100003550queensland governmenthttp://dx.doi.org/ /100008709sugar research australiaissue-copyright-statement© the author 2017 the transcriptome dataset should contribute to improved sugarcane gene models and sugarcane protein predictions; and will serve as a reference database for analysis of transcript expression in sugarcane. comparison with de novo contigs from illumina rna-sequencing of the internode samples from the same experiment and public databases showed that the iso-seq method recovered more full-length transcript isoforms, had a higher n and average length of largest proteins; whereas a greater representation of the gene content and rna diversity was captured in rna-seq. about 69% of pacbio transcript isoforms and 41% of de novo contigs aligned with the sorghum genome, indicating the high conservation of orthologs in the genic regions of the two genomes. most of the sugarcane transcriptomic studies have been based on saccharum officinarum gene indices, expressed sequence tags and de novo assembled transcript contigs from short-reads; hence knowledge of the sugarcane transcriptome is limited in relation to transcript length and number of transcript isoforms. the majority of this dataset matched the plant protein database, while just over 2% was novel transcripts, and over 2% was putative long non-coding rnas. a total of unique transcript isoforms were obtained, representing about 71% of the total number of predicted sugarcane genes. about 56% and 23% of total sequences were annotated against the gene ontology and kegg pathway databases, respectively.
 despite the economic importance of sugarcane in sugar and bioenergy production, there is not yet a reference genome available. most of the sugarcane transcriptomic studies have been based on saccharum officinarum gene indices, expressed sequence tags and de novo assembled transcript contigs from short-reads; hence knowledge of the sugarcane transcriptome is limited in relation to transcript length and number of transcript isoforms.  the sugarcane transcriptome was sequenced using pacbio isoform sequencing of a pooled rna sample derived from leaf, internode and root tissues, of different developmental stages, from varieties, to explore the potential for capturing full-length transcript isoforms. a total of unique transcript isoforms were obtained, representing about 71% of the total number of predicted sugarcane genes. the majority of this dataset matched the plant protein database, while just over 2% was novel transcripts, and over 2% was putative long non-coding rnas. about 56% and 23% of total sequences were annotated against the gene ontology and kegg pathway databases, respectively. comparison with de novo contigs from illumina rna-sequencing of the internode samples from the same experiment and public databases showed that the iso-seq method recovered more full-length transcript isoforms, had a higher n and average length of largest proteins; whereas a greater representation of the gene content and rna diversity was captured in rna-seq. only 62% of pacbio transcript isoforms matched 67% of de novo contigs, while the non-matched proportions were attributed to the inclusion of leaf/root tissues and the normalization in pacbio, and the representation of more gene content and rna classes in the de novo assembly, respectively. about 69% of pacbio transcript isoforms and 41% of de novo contigs aligned with the sorghum genome, indicating the high conservation of orthologs in the genic regions of the two genomes.  the transcriptome dataset should contribute to improved sugarcane gene models and sugarcane protein predictions; and will serve as a reference database for analysis of transcript expression in sugarcane. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords sugarcanepolyploid transcriptometranscriptome assemblyde novo assemblyisoform sequencinghybrid assemblysugit databasehttp://dx.doi.org/ /501100001794university of queenslandhttp://dx.doi.org/ /501100003550queensland governmenthttp://dx.doi.org/ /100008709sugar research australiaissue-copyright-statement© the author 2017 