BACKGROUND
whole transcriptome analysis using next generation sequencing  or sequencing-by-synthesis  technologies offers the possibility of interrogating genes and their expression en masse without knowledge of their underlying genomes.
transcriptome sequencing is often preferred over genome sequencing because of the reduced size of the sequence target space and the high functional information content  <cit> .
however, sequences generated from ngs platforms are often too short to represent entire protein-coding transcripts, and genomes for reference-guided transcriptome reconstruction are rare.
de bruijn graph assemblers allow de novo assembly of transcripts but represent only approximate computational solutions  <cit> .
the final assembly is one of many possibilities for which there is no universally accepted heuristic verification method; it is often highly redundant and contains mis-assemblies that are difficult to identify  <cit> .
post-assembly processing intended to reduce redundancy typically involves reassembly or clustering of assembled sequences.
this however may lead to propagation of mis-assemblies  <cit>  and assignment of sequences to unrelated gene clusters, resulting in loss of unique annotations  <cit> .
the main objective of transcriptome sbs is to ascribe functional labels to assembled transcribed fragments .
this is usually done via significant sequence similarity  <cit>  or domain signature annotations  <cit> .
similarity-based approaches predominantly rely on transfer of functional labels of the best blast hits to the sequence in question  <cit> .
however, low blast annotation coverage is often observed, in particularly in transcriptomes of non-model organisms  <cit> .
the implementation of significant blast hit as a proxy for functional annotation has further limitations: sequences that produce significant similarity may be functionally unrelated due to divergence  <cit> , low complexity sequences may produce high-scoring hits but have no biological relationships  <cit> , and functional homologs may lack sequence similarity  <cit> .
consequently, a first large-scale assessment of protein function shows that blast alone is often ineffective at predicting functional labels  <cit> .
domain-based annotation methods  appreciate only presence or absence of domains.
given that domains seldom function in isolation  <cit> , a reliable approach should involve a method that recognises the overall domain co-occurrence architecture of the sequences under examination.
prerequisite for domain-based annotation is a reliable protein prediction method that tolerates sequencing errors and frame shifts.
here, we introduce ifrat, which allows for selection and annotation of functionally relevant transfrags  without clustering.
this is achieved through 1) removal of identical subsequences, 2) error tolerant cds prediction, 3) identification of coding potential, and 4) complementation of blast with a multiple domain architecture annotation .
the effectiveness and versatility of this approach is shown on published datasets from non-model organisms.figure  <dig>
a schematic diagram of the ifrat pipeline.
flow diagram to illustrate the method of integrating protein-coding potential and open reading frame prediction to infer bona fide assembly derived-transscripts and multiple domain co-occurrence functional annotation.
methods
availability of supporting data
to establish a robust workflow for prioritizing and selecting functionally relevant  transfrags, we selected the fungal plant pathogen neurospora crassa  <cit>  as a species with a reference genome.
publicly available non-strand specific rna-seq data  from wild type n. crassa 74-or23-1va was obtained from the ncbi sequence read archive .
untranslated 5′ and 3′ regions were procured using ensembl biomart  <cit>  from http://fungi.ensembl.org release- <dig>  the associated genomic, predicted coding sequences and rfam family genes were obtained from the neurospora crassa sequencing project, broad institute of harvard and mit .
we verified the pipeline using recently published transcriptomes of non-model organisms: buckwheat   <cit> ; hydra   <cit> ; fresh water snail   <cit> ; centipede , marine worm , sea cradle , mediterranean sponge , and earthworm   <cit> .
the scripts, assemblies and alignment outputs generated in the ensuing analyses are available on the south africa national bioinformatics institute permanent data archive .
read preprocessing
quality scores of illumina reads generally depreciate towards the 3′-end.
prior to assembly, low quality bases were trimmed from the 3′-end of each sequence if above an error probability of  <dig>   using a custom perl script with snippets from condetri  <cit> .
the quality-based filtering and trimming process ensured that orphan reads whose partner failed the quality threshold, were retained in a separate file and used for de novo transcriptome assembly.
rna-seq assembly
reference-free transcriptome reconstruction was performed separately using either trinity , or velvet  in combination with oases .
trinity implements greedy algorithmic traversal of the k-mer graph prior to building a de bruign graph from clusters of pre-assembled sequences.
as a result, assembled transfrags are represented by actual reads.
oases on the other-hand, interrogates a pre-assembly from velvet to address alternative splicing and coverage variation across transcripts.
trinity was specifically designed for transcriptome assembly using a single, fixed k-mer size .
therefore we tested oases k- <dig> and two variations of multiple k-mer assembly: an additive assembly by pooling  as described by  <cit> , and a merged assembly using the oases-merge pipeline .
only transfrags above  <dig> bp were kept for downstream analysis.
redundancy check
a common attribute of de novo transcriptome assemblies is sequence redundancy.
using in-house perl or python with suffix array scripts, we filtered for 100% identical copies and subsequences  in both, forward and reverse directions.
to compare our filtering approach with a typically applied post-assembly clustering step, we used cd-hit-est  <cit>  with the following parameters: -m  <dig> -t  <dig> -g  <dig> -c  <dig>  - b  <dig> -al  <dig>  -as  <dig>  -n  <dig> -d  <dig> -p  <dig>  and -m  <dig> -t  <dig> -g  <dig> -c  <dig>  -b  <dig> -as  <dig>  -n  <dig> -d  <dig> -p  <dig> .
in addition, we evaluated the redundancy in each assembly using cd-hit-est as describe by  <cit> .
coding potential assessment and conceptual translation
assembled transcripts were evaluated for protein-coding attributes using portrait  <cit> .
we corrected portrait to run angle  <cit>  in  <dig> frames, since the biological orientation of transfrags from non-strand specific libraries cannot be readily ascertained.
the predicted open reading frame  with the highest dynamic programming score was chosen for conceptual translation into protein sequence using the standard codon usage table.
transfrags without an orf were classified as orphan in this study.
we note that they can be evaluated for coding capability through the protein-independent model of portrait.
transfrag annotation
we assigned protein domains to the predicted protein sequences using hmmer version  <dig>   <cit>  with the manually curated protein profile hidden markov models from pfam .
we then applied multipfam2go  <cit>  to explore co-occurrence relationships between the domains of each protein and assigned functional labels  if the underlying domain architectures predicted protein function.
to mimic annotation of non-model organisms, we generated a blast-able database of uniprot knowledgebase  release 2013_ <dig> , excluding n. crassa sequences.
we screened for highly significant blastx hits  using the ncbi blast package  and identified the top hit  using custom perl scripts.
validating bona fide transfrags by mapping to reference genome and predicted cds
the bona fide transfrags were aligned to the reference cds with blat v.  <dig>  <cit>  to assess the integrity of assembly-derived transcripts.
blat alignment in sim <dig> format were generate under intron restriction  with -minscore =  <dig> and post-alignment processing were performed through a series of custom perl scripts.
genome-base clustering was performed to assess gene space coverage by aligning bona fide transfrags to n. crassa reference genome with gmap 2013-09- <dig> v <dig>  <cit> .
the introns for n. crassa were obtained using ensembl api  <cit>  from http://fungi.ensembl.org release- <dig> to compute the maximum total length of intron per gene.
information about intron length statistics in fungi were obtain as described by  <cit>  to parameterize transfrag and cds alignment to the genome: min-intron length =  <dig>  max-intron length =  <dig>  total length =  <dig>  the known gene loci are compared to the loci of aligned transfrags in a pair-wise manner using in-house perl scripts to avoid building cluster chains  <cit> .
transfrags that do not overlap with cds are clustered using bedtools  <cit> .
we aligned sequences belonging to the 5′ and 3′ untranslated regions of predicted genes and the rfam family of predicted genes to the n. crassa genome under absolute condition of no introns using gmap with a threshold of 95% coverage and 95% identity.
the loci of these high-scoring alignments were compared to those of transfrags that did not overlap with cds.
RESULTS
de novo assembly and filtering n. crassa transfrags
the number of reads before and after quality filtering is shown in table  <dig>  of the ~ <dig> million reads that were processed,  <dig> % survived quality trimming and were retained for subsequent de novo assembly.table  <dig>
quality trimming statistics of
n. crassa
rna-seq data
attributes
raw reads
processed reads pairs
processed singletons
a summary on assembly statistics for all four assembly methods is shown in table  <dig>  when comparing the two single k-mer assembly approaches , we see that trinity produced twice as many transfrags as oases- <dig>  but at much shorter transfrags lengths.
these two assemblies had very little redundant transfrags compared to multiple k-mer assemblies.
multiple k-mer assemblies produced a much higher number of transfrags than single k-mer assemblies, but 38% - 56% were redundant.
the median transfrag lengths for these assemblies were seven-fold greater than for the trinity assembly.
to compare our filtering procedure  with a typically applied post-assembly clustering method, we used cd-hit-est and generated a non-redundant assembly at 100% global identity.
at these settings, our filtering method produced comparable results.table  <dig>
attributes of
n. crassa
assemblies produced with different filtering approaches
assembly
№ of tf
№ of unique tf
median unique tf length
% redundant tf perl
% redundant tf cd-hit-est
typically, cd-hit-est is used at settings below 100% identity.
the fraction of redundant transfrags at various identity thresholds for our n. crassa assemblies is shown in figure  <dig>  for the oases-p assembly, at 80% identity nearly 90% of the transfrags are considered redundant by cd-hit-est, ie these transfrags can be incorporated into fewer clusters.
this represents nearly  <dig>  transfrags that are lost for downstream analysis when a representative transfrag is selected for a cluster as compared to clustering at 100% identity.figure  <dig>
an assessment of redundancy in various assemblies using cd-hit-est.
comparing the fraction of redundant transfrags across all assemblies at various identity thresholds  in creating clusters with cd-hit-est.
selecting bona fide transfrag and their functional annotation
each non-redundant assembly was separated into two categories: bona fide  and orphan ; numbers are displayed in table  <dig>  in trinity, the proportion of orphan transfrags was higher  than the proportion of bona fide transfrags.
trinity also produced a considerably higher number of orphan tranfrags than any of the three oases assemblies.
as a result, the number of bona fide transfrags was very similar for the two single k-mer assemblies, and oases-p generated the highest number of bona fide transfrags.table  <dig>
classification and annotation of the non-redundant
n. crassa
transfrags
assembly
№ of utf
№ of orphan utf
№ of
bona fide
utf
№ of orphan utf with blast hit
№ of
bona fide
utf with blast hit
№ of
bona fide
utf with multipfam2go
figure  <dig> shows the distribution of transfrag lengths between bona fide and orphans transfrags.
orphan transfrags were generally much shorter than bona fide transfrags.
for the bona fide transfrags of the three oases assemblies, the median transfrag length  and the distributions are very similar.
we note that the oases assemblies had a considerable number of bona fide transfrags that were substantially longer than  <dig> kb.
the median transfrag length of bona fide transfrags assembled using trinity was  <dig>  kb, and only a few of them were longer than  <dig>  kb.figure  <dig>
distribution of transfrag length for various assemblies.
bona fide transfrags : transfrag with coding potential and predicted cds; orphan : transfrags with no coding potential or coding potential but no predicted cds.
non-redundant assemblies were annotated using blast and multipfam2go .
we note that in all assemblies only a small proportion of orphan transfrags had a blast match.
despite the highest number of orphan transfrags, trinity had the least number of blast hits to transfrags in this category.
in contrast, at least 70% of bona fide transfrags from all assemblies had a blast hit.
this represented over 94% of cumulative blastx retrievable hit of the unfiltered assembly .
this number is higher than the ones typically reported in studies on de novo assembled transcriptomes  <cit> .
in addition, bona fide transfrags were annotated with multipfam2go.
the fraction of transfrags that could be associated with gene ontology terms ranged from 33%-50%, which is also high for domain based annotation  <cit> .
assessing transfrag integrity and gene coverage
to evaluate the number of predicted genes represented by the bona fide transfrags, we aligned them to the predicted coding sequences  as well as to the genome of n. crassa .
between 80% and 90% of the bona fide tranfrags mapped to both datasets at high stringency.
although the numbers of bona fide transfrags between single and multiple k-mer assemblies is very different, the number of identified genes is very similar.
most strikingly, trinity identified the same number of predicted genes and putative unknown n. crassa gene loci as oases-p, independent of the dataset and the alignment thresholds.
as a result, the number of bona fide transfrags per gene is lower in single k-mer versus multiple k-mer assemblies.
orphan transfrags that mapped at the same stringency represented 15-40% of the known gene loci , but ~ 90% were already identified by the longer bona fide categories.
unmapped transfrags mapped to multiple location and some were chimeric.
the number of loci occupied by orphan transfrags ranged from  <dig>  -  <dig> .
a look into the biological relevance of these loci revealed that they represent intronic, 5′ and 3′ untranslated regions and genes of the rfam family .table  <dig>
summary of
bona fide
† and orphan* transfrags integrity and validity
assembly
№ of
bona fide
utf
№ of chimeras in unmapped transfrags ζ
alignment of tf to reference genes
alignment of tf to reference genome
№ of tf cov 50%, id 50%
№ of reference unigenes
№ of tf cov 90%, id 90%
№ of reference unigenes
№ of tf uniquely mapped
№ of
n. crassa
genes identified by tf
trinity†
oases-25†
oases-m†
oases-p†
†bona fide: transfrags with coding potential and predicted cds.
*orphan: transfrags with no coding potential or with coding potential but no predicted cds.
ζ possible chimera with a distinct breakpoint.
selecting bona fide assembly-derived transcripts in other species
we also verified the suitability of the ifrat pipeline for selecting reconstructed transcripts in non-model organisms.
the analysis results for unique transfrags longer than  <dig> bp from each draft assembly are show in table  <dig>  we predicted that up to 70% of the published transcriptome could be categorized as orphan transfrags.
as before, the percentage of orphan transfrags with a blast hit was relatively low.
in contrast, the proportion of bona fide transfrags with significant blast matches was often higher than in the unfiltered draft assemblies.table  <dig>
allocation of blastx hits between
bona fide
and orphan transfrags inferred with ifrat
organism
№ of tf in publication
№ of tf with hit in publication
№ of utf > = 100
№ of orphan utf
№ of orphan utf with blast hit
№ of bona fide utf
№ of bona fide with blast hit
hydra vulgaris
radix balthica
alipes grandidieri
cerebratulus marginatus
chiton olivaceus
crella elegans
hormogaster samnitica
fagopyrum tataricum
discussion
single k-mer assemblies of transcriptomes are considered incomplete because a short k-mer result in a highly diverse but also fragmented and redundant assembly, while a long k-mer provides a more contiguous assembly but misses poorly expressed transcripts  <cit> .
to account for this problem the multiple k-mer transcriptome assembly approach was introduced  <cit> .
however, the number of sequences generated in this way exceeds by far the number of protein coding genes likely to exist in the respective organism  <cit> , making identification of genuine transfrags a major challenge for downstream analysis.
to reduce redundancy, clustering or merging methods are currently being applied  <cit> .
yet, these methods rely on common word heuristics, ignoring the biological nature of assembled transcripts  <cit> .
therefore, reference-free clustering tends to mis-assign transfrags to biologically unrelated clusters  <cit>  which leads to loss of unique functional annotations  <cit>  and creation of chimeric transcripts  <cit> .
here, we propose ifrat, a workflow that allows selection of unique bona fide transfrags  without clustering; and introduce domain co-occurrence analysis as means of tranfsrag assembly verification.
ifrat filters unique transfrags by removing exact duplicates, including identical forward and reverse nucleotide subsequences.
ifrat filtering removes slightly more transfrags than cd-hit-est at 100% identity because this program does not properly process transfrags containing n’s .
our results suggest that single k-mer assemblies may not need this filtering step since the proportion of redundant transfrags in the trinity and oases- <dig> datasets were only about 1%.
in contrast, redundancy filtering is particularly important in multiple k-mer assemblies, considering that nearly half the transfrags in the oases-m and oases-p datasets were exact copies or substrings of other transfrags.
it is unknown at what percent identity clustering results in significant loss of unique functional annotations.
however, as suggested by our analysis, clustering without biological insight should be handled with caution because at 99% identity a significant subset of potentially unique transfrags is removed by cd-hit-est.
after filtering, ifrat classifies the sequences into bona fide and orphan transfrags based on cds prediction and coding potential.
our subsequent blast analysis corroborated this categorization, since 70-80% of bona fide transfrags had significant blast matches while this was true for only 1-6% of orphan tranfrags.
we note that the median length of bona fide transfrags ranged from  <dig> kb  to 2kb , which is consistent with the average coding sequence length in fungi  <cit>  while most of the orphan transfrags were short .
however, our results confirmed previous findings that length is not the only indicator of coding potential  <cit>  and ‘non-blastable’ transfrags  <cit> , since 6-26% of the orphan transfrags with blast matches were less than  <dig> bp long.
all four assembly methods produced high quality datasets, as 76-90% of the transfrags mapped to the genome and the predicted cds of n. crassa at high identity and coverage.
bona fide transfrags represented approximately 73% of the  <dig>  known gene loci in the n. crassa genome.
in addition, they indicated the existence of 715- <dig> unknown potentially coding gene locations.
orphan transfrags also mapped to known gene locations, but most of these locations were represented by longer bona fide transfrags.
these orphan transfrags may represent biologically interesting data, such as truncated assemblies , or immature mrna with intronic regions and long utrs for which coding potential could not be predicted  <cit> .
orphan transfrags that mapped to non-coding regions of the genomes represented ribosomal and non-coding rna  <cit> , and may be of interest.
in any case, it is advisable to verify the correct assembly of orphan transfrags and remove mis-assemblies using a suitable reference dataset, such as a reference genome or est collection.
we integrated multi-domain co-occurrence architecture  <cit>  to complement blast annotation.
this avoids non-specific annotation of promiscuous domains resulting from truncated transfrags.
between 44% and 50% of the bona fide tranfrag peptides from n. crassa were assigned at least one go term.
using ifrat, we improved annotation coverage of published transcriptome datasets from non-model organisms.
the choice of database and to a larger extend the coverage filter threshold accounts for small differences in the number of blast hits between bona fide transfrags and unfiltered assemblies.
we attribute this high annotation coverage to the error tolerant cds prediction  <cit>  and selection of longer proteins with coding potential by ifrat.
ifrat is able to select bona fide transfrags irrespective of the assembler or assembly method used.
profound differences between transcriptome assemblers and assembly methods have been elaborately dealt with elsewhere  <cit> .
we note however that trinity performed very similar to oases-p in identifying cds and known gene loci, requiring substantially less computational resources.
other technical limitations, such as runtime and data-size, may influence the choice of one assembler over the other  <cit> .
since many more transfrags were produced by the multiple k-mer assemblies that identified a comparable number of gene loci, they may be suitable for studies with interest on splice variants.
CONCLUSIONS
we have proposed a framework for post-assembly analysis of transcriptome assembly that is flexible enough to accommodate sequencing error and frame-shifts and that does not rely on the availability of a reference genome.
through this, a catalogue of reliable protein coding transfrags is established that represents a reference transcriptome.
the method described herein is potentially applicable not only to assemblies of transcribed fragments generated with trinity or oases but also to assemblies produced by de bruijn graph assemblers where no reliable sequenced genome is available, as demonstrated with the published datasets.
our framework performs well in segregating functionally relevant transcripts.
we note that the main challenge remains the quality of assembly-derived-transcripts, which is undermined by incorporation of non-coding fragments that reduces the coding potential signal.
one possible avenue for improvement is to screen rna-seq reads for non-coding transcribed fragments prior to de novo assembly.
additional files
additional file 1:
distribution of blastx hits between
bona fide
and orphan transfrags.
the bona fide transfrags are enriched with sequences that have a potential blast hit .
additional file 2:
number of orphan transfrags that overlap with genic features and non-protein coding genes .
competing interests
the authors declare that they have no competing interests.
authors’ contributions
skm conceived, design, performed the experiments and analyzed the data.
skm and uh interpreted the results and drafted the manuscript.
pvh implemented redundancy removal in python with suffix array.
ac critically evaluated the manuscript, provided reagents, materials, participated in design and supervised the implementation.
smk, uh, pvh and ac approved the final manuscript.

 whole transcriptome analysis using next generation sequencing or sequencing-by-synthesis technologies offers the possibility of interrogating genes and their expression en masse without knowledge of their underlying genomes. transcriptome sequencing is often preferred over genome sequencing because of the reduced size of the sequence target space and the high functional information content. however, sequences generated from ngs platforms are often too short to represent entire protein-coding transcripts, and genomes for reference-guided transcriptome reconstruction are rare. de bruijn graph assemblers allow de novo assembly of transcripts but represent only approximate computational solutions. the final assembly is one of many possibilities for which there is no universally accepted heuristic verification method; it is often highly redundant and contains mis-assemblies that are difficult to identify. post-assembly processing intended to reduce redundancy typically involves reassembly or clustering of assembled sequences. this however may lead to propagation of mis-assemblies and assignment of sequences to unrelated gene clusters, resulting in loss of unique annotations. the main objective of transcriptome sbs is to ascribe functional labels to assembled transcribed fragments. this is usually done via significant sequence similarity or domain signature annotations. similarity-based approaches predominantly rely on transfer of functional labels of the best blast hits to the sequence in question. however, low blast annotation coverage is often observed, in particularly in transcriptomes of non-model organisms. the implementation of significant blast hit as a proxy for functional annotation has further limitations: sequences that produce significant similarity may be functionally unrelated due to divergence, low complexity sequences may produce high-scoring hits but have no biological relationships, and functional homologs may lack sequence similarity. consequently, a first large-scale assessment of protein function shows that blast alone is often ineffective at predicting functional labels. domain-based annotation methods appreciate only presence or absence of domains. given that domains seldom function in isolation, a reliable approach should involve a method that recognises the overall domain co-occurrence architecture of the sequences under examination. prerequisite for domain-based annotation is a reliable protein prediction method that tolerates sequencing errors and frame shifts. here, we introduce ifrat, which allows for selection and annotation of functionally relevant transfrags without clustering. this is achieved through 1) removal of identical subsequences, 2) error tolerant cds prediction, 3) identification of coding potential, and 4) complementation of blast with a multiple domain architecture annotation. the effectiveness and versatility of this approach is shown on published datasets from non-model organisms.figure a schematic diagram of the ifrat pipeline. flow diagram to illustrate the method of integrating protein-coding potential and open reading frame prediction to infer bona fide assembly derived-transscripts and multiple domain co-occurrence functional annotation. methods availability of supporting data to establish a robust workflow for prioritizing and selecting functionally relevant transfrags, we selected the fungal plant pathogen neurospora crassa as a species with a reference genome. publicly available non-strand specific rna-seq data from wild type n. crassa 74-or23-1va was obtained from the ncbi sequence read archive. untranslated 5′ and 3′ regions were procured using ensembl biomart from http://fungi.ensembl.org release- the associated genomic, predicted coding sequences and rfam family genes were obtained from the neurospora crassa sequencing project, broad institute of harvard and mit. we verified the pipeline using recently published transcriptomes of non-model organisms: buckwheat; hydra; fresh water snail; centipede, marine worm, sea cradle, mediterranean sponge, and earthworm. the scripts, assemblies and alignment outputs generated in the ensuing analyses are available on the south africa national bioinformatics institute permanent data archive. read preprocessing quality scores of illumina reads generally depreciate towards the 3′-end. prior to assembly, low quality bases were trimmed from the 3′-end of each sequence if above an error probability of using a custom perl script with snippets from condetri. the quality-based filtering and trimming process ensured that orphan reads whose partner failed the quality threshold, were retained in a separate file and used for de novo transcriptome assembly. rna-seq assembly reference-free transcriptome reconstruction was performed separately using either trinity, or velvet in combination with oases. trinity implements greedy algorithmic traversal of the k-mer graph prior to building a de bruign graph from clusters of pre-assembled sequences. as a result, assembled transfrags are represented by actual reads. oases on the other-hand, interrogates a pre-assembly from velvet to address alternative splicing and coverage variation across transcripts. trinity was specifically designed for transcriptome assembly using a single, fixed k-mer size. therefore we tested oases k- and two variations of multiple k-mer assembly: an additive assembly by pooling as described by, and a merged assembly using the oases-merge pipeline. only transfrags above bp were kept for downstream analysis. redundancy check a common attribute of de novo transcriptome assemblies is sequence redundancy. using in-house perl or python with suffix array scripts, we filtered for 100% identical copies and subsequences in both, forward and reverse directions. to compare our filtering approach with a typically applied post-assembly clustering step, we used cd-hit-est with the following parameters: -m -t -g -c - b -al -as -n -d -p and -m -t -g -c -b -as -n -d -p. in addition, we evaluated the redundancy in each assembly using cd-hit-est as describe by. coding potential assessment and conceptual translation assembled transcripts were evaluated for protein-coding attributes using portrait. we corrected portrait to run angle in frames, since the biological orientation of transfrags from non-strand specific libraries cannot be readily ascertained. the predicted open reading frame with the highest dynamic programming score was chosen for conceptual translation into protein sequence using the standard codon usage table. transfrags without an orf were classified as orphan in this study. we note that they can be evaluated for coding capability through the protein-independent model of portrait. transfrag annotation we assigned protein domains to the predicted protein sequences using hmmer version with the manually curated protein profile hidden markov models from pfam. we then applied multipfam2go to explore co-occurrence relationships between the domains of each protein and assigned functional labels if the underlying domain architectures predicted protein function. to mimic annotation of non-model organisms, we generated a blast-able database of uniprot knowledgebase release 2013_, excluding n. crassa sequences. we screened for highly significant blastx hits using the ncbi blast package and identified the top hit using custom perl scripts. validating bona fide transfrags by mapping to reference genome and predicted cds the bona fide transfrags were aligned to the reference cds with blat v. to assess the integrity of assembly-derived transcripts. blat alignment in sim format were generate under intron restriction with -minscore =  and post-alignment processing were performed through a series of custom perl scripts. genome-base clustering was performed to assess gene space coverage by aligning bona fide transfrags to n. crassa reference genome with gmap 2013-09- v. the introns for n. crassa were obtained using ensembl api from http://fungi.ensembl.org release- to compute the maximum total length of intron per gene. information about intron length statistics in fungi were obtain as described by to parameterize transfrag and cds alignment to the genome: min-intron length =  max-intron length =  total length =  the known gene loci are compared to the loci of aligned transfrags in a pair-wise manner using in-house perl scripts to avoid building cluster chains. transfrags that do not overlap with cds are clustered using bedtools. we aligned sequences belonging to the 5′ and 3′ untranslated regions of predicted genes and the rfam family of predicted genes to the n. crassa genome under absolute condition of no introns using gmap with a threshold of 95% coverage and 95% identity. the loci of these high-scoring alignments were compared to those of transfrags that did not overlap with cds.  de novo assembly and filtering n. crassa transfrags the number of reads before and after quality filtering is shown in table  of the ~ million reads that were processed, % survived quality trimming and were retained for subsequent de novo assembly.table quality trimming statistics of n. crassa rna-seq data attributes raw reads processed reads pairs processed singletons a summary on assembly statistics for all four assembly methods is shown in table  when comparing the two single k-mer assembly approaches, we see that trinity produced twice as many transfrags as oases- but at much shorter transfrags lengths. these two assemblies had very little redundant transfrags compared to multiple k-mer assemblies. multiple k-mer assemblies produced a much higher number of transfrags than single k-mer assemblies, but 38% - 56% were redundant. the median transfrag lengths for these assemblies were seven-fold greater than for the trinity assembly. to compare our filtering procedure with a typically applied post-assembly clustering method, we used cd-hit-est and generated a non-redundant assembly at 100% global identity. at these settings, our filtering method produced comparable results.table attributes of n. crassa assemblies produced with different filtering approaches assembly № of tf № of unique tf median unique tf length % redundant tf perl % redundant tf cd-hit-est typically, cd-hit-est is used at settings below 100% identity. the fraction of redundant transfrags at various identity thresholds for our n. crassa assemblies is shown in figure  for the oases-p assembly, at 80% identity nearly 90% of the transfrags are considered redundant by cd-hit-est, ie these transfrags can be incorporated into fewer clusters. this represents nearly transfrags that are lost for downstream analysis when a representative transfrag is selected for a cluster as compared to clustering at 100% identity.figure an assessment of redundancy in various assemblies using cd-hit-est. comparing the fraction of redundant transfrags across all assemblies at various identity thresholds in creating clusters with cd-hit-est. selecting bona fide transfrag and their functional annotation each non-redundant assembly was separated into two categories: bona fide and orphan; numbers are displayed in table  in trinity, the proportion of orphan transfrags was higher than the proportion of bona fide transfrags. trinity also produced a considerably higher number of orphan tranfrags than any of the three oases assemblies. as a result, the number of bona fide transfrags was very similar for the two single k-mer assemblies, and oases-p generated the highest number of bona fide transfrags.table classification and annotation of the non-redundant n. crassa transfrags assembly № of utf № of orphan utf № of bona fide utf № of orphan utf with blast hit № of bona fide utf with blast hit № of bona fide utf with multipfam2go figure  shows the distribution of transfrag lengths between bona fide and orphans transfrags. orphan transfrags were generally much shorter than bona fide transfrags. for the bona fide transfrags of the three oases assemblies, the median transfrag length and the distributions are very similar. we note that the oases assemblies had a considerable number of bona fide transfrags that were substantially longer than kb. the median transfrag length of bona fide transfrags assembled using trinity was kb, and only a few of them were longer than kb.figure distribution of transfrag length for various assemblies. bona fide transfrags: transfrag with coding potential and predicted cds; orphan: transfrags with no coding potential or coding potential but no predicted cds. non-redundant assemblies were annotated using blast and multipfam2go. we note that in all assemblies only a small proportion of orphan transfrags had a blast match. despite the highest number of orphan transfrags, trinity had the least number of blast hits to transfrags in this category. in contrast, at least 70% of bona fide transfrags from all assemblies had a blast hit. this represented over 94% of cumulative blastx retrievable hit of the unfiltered assembly. this number is higher than the ones typically reported in studies on de novo assembled transcriptomes. in addition, bona fide transfrags were annotated with multipfam2go. the fraction of transfrags that could be associated with gene ontology terms ranged from 33%-50%, which is also high for domain based annotation. assessing transfrag integrity and gene coverage to evaluate the number of predicted genes represented by the bona fide transfrags, we aligned them to the predicted coding sequences as well as to the genome of n. crassa. between 80% and 90% of the bona fide tranfrags mapped to both datasets at high stringency. although the numbers of bona fide transfrags between single and multiple k-mer assemblies is very different, the number of identified genes is very similar. most strikingly, trinity identified the same number of predicted genes and putative unknown n. crassa gene loci as oases-p, independent of the dataset and the alignment thresholds. as a result, the number of bona fide transfrags per gene is lower in single k-mer versus multiple k-mer assemblies. orphan transfrags that mapped at the same stringency represented 15-40% of the known gene loci, but ~ 90% were already identified by the longer bona fide categories. unmapped transfrags mapped to multiple location and some were chimeric. the number of loci occupied by orphan transfrags ranged from -. a look into the biological relevance of these loci revealed that they represent intronic, 5′ and 3′ untranslated regions and genes of the rfam family.table summary of bona fide † and orphan* transfrags integrity and validity assembly № of bona fide utf № of chimeras in unmapped transfrags ζ alignment of tf to reference genes alignment of tf to reference genome № of tf cov 50%, id 50% № of reference unigenes № of tf cov 90%, id 90% № of reference unigenes № of tf uniquely mapped № of n. crassa genes identified by tf trinity† oases-25† oases-m† oases-p† †bona fide: transfrags with coding potential and predicted cds. *orphan: transfrags with no coding potential or with coding potential but no predicted cds. ζ possible chimera with a distinct breakpoint. selecting bona fide assembly-derived transcripts in other species we also verified the suitability of the ifrat pipeline for selecting reconstructed transcripts in non-model organisms. the analysis results for unique transfrags longer than bp from each draft assembly are show in table  we predicted that up to 70% of the published transcriptome could be categorized as orphan transfrags. as before, the percentage of orphan transfrags with a blast hit was relatively low. in contrast, the proportion of bona fide transfrags with significant blast matches was often higher than in the unfiltered draft assemblies.table allocation of blastx hits between bona fide and orphan transfrags inferred with ifrat organism № of tf in publication № of tf with hit in publication № of utf > = 100 № of orphan utf № of orphan utf with blast hit № of bona fide utf № of bona fide with blast hit hydra vulgaris radix balthica alipes grandidieri cerebratulus marginatus chiton olivaceus crella elegans hormogaster samnitica fagopyrum tataricum discussion single k-mer assemblies of transcriptomes are considered incomplete because a short k-mer result in a highly diverse but also fragmented and redundant assembly, while a long k-mer provides a more contiguous assembly but misses poorly expressed transcripts. to account for this problem the multiple k-mer transcriptome assembly approach was introduced. however, the number of sequences generated in this way exceeds by far the number of protein coding genes likely to exist in the respective organism, making identification of genuine transfrags a major challenge for downstream analysis. to reduce redundancy, clustering or merging methods are currently being applied. yet, these methods rely on common word heuristics, ignoring the biological nature of assembled transcripts. therefore, reference-free clustering tends to mis-assign transfrags to biologically unrelated clusters which leads to loss of unique functional annotations and creation of chimeric transcripts. here, we propose ifrat, a workflow that allows selection of unique bona fide transfrags without clustering; and introduce domain co-occurrence analysis as means of tranfsrag assembly verification. ifrat filters unique transfrags by removing exact duplicates, including identical forward and reverse nucleotide subsequences. ifrat filtering removes slightly more transfrags than cd-hit-est at 100% identity because this program does not properly process transfrags containing n’s. our results suggest that single k-mer assemblies may not need this filtering step since the proportion of redundant transfrags in the trinity and oases- datasets were only about 1%. in contrast, redundancy filtering is particularly important in multiple k-mer assemblies, considering that nearly half the transfrags in the oases-m and oases-p datasets were exact copies or substrings of other transfrags. it is unknown at what percent identity clustering results in significant loss of unique functional annotations. however, as suggested by our analysis, clustering without biological insight should be handled with caution because at 99% identity a significant subset of potentially unique transfrags is removed by cd-hit-est. after filtering, ifrat classifies the sequences into bona fide and orphan transfrags based on cds prediction and coding potential. our subsequent blast analysis corroborated this categorization, since 70-80% of bona fide transfrags had significant blast matches while this was true for only 1-6% of orphan tranfrags. we note that the median length of bona fide transfrags ranged from kb to 2kb, which is consistent with the average coding sequence length in fungi while most of the orphan transfrags were short. however, our results confirmed previous findings that length is not the only indicator of coding potential and ‘non-blastable’ transfrags, since 6-26% of the orphan transfrags with blast matches were less than bp long. all four assembly methods produced high quality datasets, as 76-90% of the transfrags mapped to the genome and the predicted cds of n. crassa at high identity and coverage. bona fide transfrags represented approximately 73% of the known gene loci in the n. crassa genome. in addition, they indicated the existence of 715- unknown potentially coding gene locations. orphan transfrags also mapped to known gene locations, but most of these locations were represented by longer bona fide transfrags. these orphan transfrags may represent biologically interesting data, such as truncated assemblies, or immature mrna with intronic regions and long utrs for which coding potential could not be predicted. orphan transfrags that mapped to non-coding regions of the genomes represented ribosomal and non-coding rna, and may be of interest. in any case, it is advisable to verify the correct assembly of orphan transfrags and remove mis-assemblies using a suitable reference dataset, such as a reference genome or est collection. we integrated multi-domain co-occurrence architecture to complement blast annotation. this avoids non-specific annotation of promiscuous domains resulting from truncated transfrags. between 44% and 50% of the bona fide tranfrag peptides from n. crassa were assigned at least one go term. using ifrat, we improved annotation coverage of published transcriptome datasets from non-model organisms. the choice of database and to a larger extend the coverage filter threshold accounts for small differences in the number of blast hits between bona fide transfrags and unfiltered assemblies. we attribute this high annotation coverage to the error tolerant cds prediction and selection of longer proteins with coding potential by ifrat. ifrat is able to select bona fide transfrags irrespective of the assembler or assembly method used. profound differences between transcriptome assemblers and assembly methods have been elaborately dealt with elsewhere. we note however that trinity performed very similar to oases-p in identifying cds and known gene loci, requiring substantially less computational resources. other technical limitations, such as runtime and data-size, may influence the choice of one assembler over the other. since many more transfrags were produced by the multiple k-mer assemblies that identified a comparable number of gene loci, they may be suitable for studies with interest on splice variants.  we have proposed a framework for post-assembly analysis of transcriptome assembly that is flexible enough to accommodate sequencing error and frame-shifts and that does not rely on the availability of a reference genome. through this, a catalogue of reliable protein coding transfrags is established that represents a reference transcriptome. the method described herein is potentially applicable not only to assemblies of transcribed fragments generated with trinity or oases but also to assemblies produced by de bruijn graph assemblers where no reliable sequenced genome is available, as demonstrated with the published datasets. our framework performs well in segregating functionally relevant transcripts. we note that the main challenge remains the quality of assembly-derived-transcripts, which is undermined by incorporation of non-coding fragments that reduces the coding potential signal. one possible avenue for improvement is to screen rna-seq reads for non-coding transcribed fragments prior to de novo assembly. additional files additional file 1: distribution of blastx hits between bona fide and orphan transfrags. the bona fide transfrags are enriched with sequences that have a potential blast hit. additional file 2: number of orphan transfrags that overlap with genic features and non-protein coding genes. competing interests the authors declare that they have no competing interests. authors’ contributions skm conceived, design, performed the experiments and analyzed the data. skm and uh interpreted the results and drafted the manuscript. pvh implemented redundancy removal in python with suffix array. ac critically evaluated the manuscript, provided reagents, materials, participated in design and supervised the implementation. smk, uh, pvh and ac approved the final manuscript. 
the median length of bona fide transfrags ranged from kb to 2kb, which is consistent with the average coding sequence length in fungi. de novo transcriptome assembly of short transcribed fragments produced from sequencing-by-synthesis technologies often results in redundant datasets with differing levels of unassembled, partially assembled or mis-assembled transcripts. we demonstrate that independent of the assembler, ifrat selects bona fide transfrags from the transcriptome assembly of a model organism without relying on post-assembly clustering or reassembly. ifrat combines 1) removal of identical subsequences, 2) error tolerant cds prediction, 3) identification of coding potential, and 4) complements blast with a multiple domain architecture annotation that reduces non-specific domain annotation. single k-mer assemblies contained fewer transfrags compared to the multiple k-mer assemblies. the fraction of transfrags that could be associated with gene ontology terms ranged from 33-50%, which is also high for domain based annotation. however, these approaches are mostly based on common word heuristics and often create clusters of biologically unrelated sequences, resulting in loss of unique transfrags annotations and propagation of mis-assemblies. ifrat simplifies post-assembly processing providing a reference transcriptome enriched with functionally relevant assembly-derived transcripts for non-model organism. however, trinity identified a comparable number of predicted coding sequence and gene loci to oases pooled assembly. we showed that unselected transfrags were mostly truncated and represent sequences from intronic, untranslated regions and non-coding gene loci. the robustness of ifrat is inferred on rna-seq data of neurospora crassa assembled using de bruijn graph-based assemblers, in single and multiple k-mer modes. between 4-6% are lost when orphan transfrags are excluded and this represents only a tiny fraction of annotation derived from functional transference by sequence similarity. ifrat selects bona fide transfrags representing over 94% of cumulative blast-derived functional annotations of the unfiltered assemblies.
 de novo transcriptome assembly of short transcribed fragments produced from sequencing-by-synthesis technologies often results in redundant datasets with differing levels of unassembled, partially assembled or mis-assembled transcripts. post-assembly processing intended to reduce redundancy typically involves reassembly or clustering of assembled sequences. however, these approaches are mostly based on common word heuristics and often create clusters of biologically unrelated sequences, resulting in loss of unique transfrags annotations and propagation of mis-assemblies.  here, we propose a structured framework that consists of a few steps in pipeline architecture for inferring functionally relevant assembly-derived transcripts. ifrat combines 1) removal of identical subsequences, 2) error tolerant cds prediction, 3) identification of coding potential, and 4) complements blast with a multiple domain architecture annotation that reduces non-specific domain annotation. we demonstrate that independent of the assembler, ifrat selects bona fide transfrags from the transcriptome assembly of a model organism without relying on post-assembly clustering or reassembly. the robustness of ifrat is inferred on rna-seq data of neurospora crassa assembled using de bruijn graph-based assemblers, in single and multiple k-mer modes. single k-mer assemblies contained fewer transfrags compared to the multiple k-mer assemblies. however, trinity identified a comparable number of predicted coding sequence and gene loci to oases pooled assembly. ifrat selects bona fide transfrags representing over 94% of cumulative blast-derived functional annotations of the unfiltered assemblies. between 4-6% are lost when orphan transfrags are excluded and this represents only a tiny fraction of annotation derived from functional transference by sequence similarity. the median length of bona fide transfrags ranged from kb to 2kb, which is consistent with the average coding sequence length in fungi. the fraction of transfrags that could be associated with gene ontology terms ranged from 33-50%, which is also high for domain based annotation. we showed that unselected transfrags were mostly truncated and represent sequences from intronic, untranslated regions and non-coding gene loci.  ifrat simplifies post-assembly processing providing a reference transcriptome enriched with functionally relevant assembly-derived transcripts for non-model organism. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords transcriptome reconstructiontransfragscoding potentialmultiple domainannotationissue-copyright-statement© the author 2015 