BACKGROUND
the sequencing and in-silico analysis of messenger rna  is now routinely being applied to complex microbial communities in diverse eco-systems, including, but not limited to: soil
, marine
and intestinal
<cit>  habitats.
the typical goals of metatranscriptomics are to taxonomically classify transcripts, predict their functions and quantify their abundances, and to relate these to environmental data in order to reveal how environmental conditions impact microbial communities in different habitats.
metatranscriptome data sets typically consist of hundreds of thousands of  <dig> sequences, or, more recently tens of millions of illumina sequences per sample.
low taxonomic diversity and/or highly expressed genes can lead to a high degree of data redundancy; that is multiple identical or nearly identical sequence fragments.
in an investigation into the proportion of artificial and natural duplicates in pyrosequenced metatranscriptome data, niu et al.
reported that as much as 60% of all sequences in an early metatranscriptome data set were likely natural duplicates
<cit> .
therefore, some form of data reduction strategy is beneficial before running computationally intensive homology searches.
two approaches that are commonly employed to reduce redundancy in large data sets are  assembly: where sequences are assembled into longer contiguous fragments  and  clustering: sequences are grouped into clusters sharing a defined degree of similarity.
the decisions as to whether to perform data reduction and which method to employ are influenced by several factors:  the availability of reference genomes: if sufficient reference genomes are available for a small number of dominant species then the sequences can be mapped to them and taxonomy and function inferred and the relative abundance of the transcripts calculated.
read length - are the unprocessed reads long enough to return annotations?
current illumina platforms produce shorter reads than  <dig>  and are likely to return a lower hit rate to protein databases compared to longer  <dig> reads
<cit> .
the diversity of the sample: although assembly can produce longer sequences and increase the accuracy of subsequent annotations, the variable coverage of transcripts in metatranscriptomics data sets and the presence of closely related organisms can lead to chimeric contigs.
indeed, for highly diverse metagenomic samples it has been recommended that assembly not be performed at all
<cit> .
the aims of the analysis: if the read length is adequate for annotation and the intention is to count features  then clustering at high identities is a recommended alternative
<cit> .
with the lower coverage but higher read length of  <dig> metatranscriptome data, assembly is relatively uncommon and instead authors tend to either cluster or annotate sequences individually.
clustering is regularly used for detecting and removing sequencing artifacts from  <dig> data
<cit> , grouping rrna data into operational taxonomic units
<cit> , and grouping proteins into families
<cit> .
in addition to the known benefits of a reduction in the size of the data set and therefore computation time, we set out to assess whether, by clustering translated metatranscriptome sequences and transferring protein domain annotation from cluster representatives to cluster members - some of which may only partially cover protein domains used for classification, we can accurately increase the number of classifiable reads.
more specifically, we investigated some popular data reduction tools and assessed their performance on simulated  <dig> and illumina metatranscriptome data in terms of the accuracy of resulting protein annotations.
note that although several approaches have been described to simulate metagenomic data sets
and rna-seq data
<cit> , to date only small scale attempts have been made to simulate metatranscriptome data sets based on a small number of species
<cit> .
RESULTS
simulated  <dig> data
the simulated  <dig> data sets contained  <dig>  sequences each, totalling ~ <dig> megabases of sequence per diversity level.
between  <dig> and 14% of  <dig> sequences from each data set returned matches to pfam-a.
when compared to the theoretical domain content, the correlation coefficients for all read annotation were  <dig> ,  <dig>  and  <dig>  for ld, md and hd respectively .table  <dig>
correlation coefficients between simulated data set annotations and known protein domain content
all
clustered
contigs
debris
assembly
1
clustered assembly
summary of pearson correlation coefficients between processed data sets and the known domain content of sample for low diversity , medium diversity  and high diversity  simulated  <dig> and illumina metatranscriptomes.
1assembly includes annotation from both contigs and debris sequences.
then, taking the parameter set that provided the largest increase in true positives minus false positives, compared to the annotation of all unclustered reads, we found that the best clustering parameters were: ≥ 60% overall similarity and 100% coverage of cluster member sequences for the ld data set; ≥80% similarity and 100% coverage of the cluster members for the md data set; and ≥60% similarity, ≥25% coverage of the cluster representative and between 0-50% minimum coverage of cluster members for the hd data set .
while the best performing clustering parameters produced a net gain  of between  <dig>  and  <dig>  domains , the correlation coefficients were slightly lower than for all read annotation ).
the mira assemblies incorporated ~50% of all sequences into  <dig>  and  <dig>  contigs for the ld and md samples respectively, and ~30% of sequences into  <dig>  contigs for the hd sample.
the average contig lengths were  <dig> ,  <dig>  and  <dig>  base pairs for ld, md and hd, respectively .
the average contig entropy was  <dig> ,  <dig>  and  <dig>  for ld, md and hd respectively  with  <dig> %,  <dig> % and  <dig> % of contigs possessing an entropy of zero.figure  <dig>
contig entropy for assembled simulated metatranscriptomes.
contig entropy plotted against contig length for a) mira assembled simulated  <dig> data sets and b) trinity assembled simulated illumina data sets.
plots represent, from left to right: low diversity , medium diversity  and high diversity  data sets.
for the ld and md data sets, the net gain of true positives  was a ~100% increase, and for the hd data set an increase of ~20% was achieved .
the contigs alone had a weaker correlation with the theoretical domain content than all read or clustered read annotation .
when combined with the debris sequences, the correlation coefficients for all three samples were higher than for all all-read or clustered annotations ).
this could be due to two factors: firstly the low proportion of sequences incorporated into the contigs,  and secondly the assemblies may be biased towards high-abundance transcripts .figure  <dig>
results from pfam-a annotated simulated metatranscriptomes.
percentage of true positives, false positives, true negatives and potential domains  based on a comparison with the known domain content of the data sets for all reads , best clustering , assembly  and clustered assembly .
a) results for simulated  <dig> data sets, from left to right: low, medium and high diversity.
b) results for simulated illumina data sets from left to right: low, medium and high diversity.
correlation between high diversity simulations and known protein domain content.
correlation plots of pfam-a annotations of each processed data set compared to known domain content for a) high diversity  <dig> simulated data set and b) high diversity illumina simulated data set.
top row, left to right: all reads unprocessed; clustered reads; assembly - contigs only.
bottom row, left to right: assembly – debris only; assembly – contigs and debris combined; clustered assembly.
pearson correlation coefficient shown in top left corner.
clustering of the  <dig> assemblies  led to a very slight increase in the detection of true positives  but the overall effect was a very slight reduction in the correlation with the theoretical domain content compared to the unclustered assembly .
simulated illumina data
around 4% of the illumina reads could be annotated with pfam-a domains.
the correlation coefficients for all read annotation with the theoretical domain content were .
the illumina data sets were clustered with the best performing parameter set for the equivalent diversity level identified in the  <dig> simulations described above.
while clustering reduced the data sets by ~40% for ld and md and ~25% for the hd data set the resulting annotations had a weaker correlation to the theoretical domain content of the sample .the trinity assemblies incorporated ~40% of sequences from the ld and md data sets into  <dig>  and  <dig>  contigs respectively with an average length of ~400 nt.
for the hd data set, ~14% of reads from the hd data set into  <dig>  contigs with an average length of 328 nt.
the average contig entropy was  <dig> ,  <dig>  and  <dig>  for ld, md and hd respectively  with  <dig> %,  <dig> % and 92% of contigs possessing an entropy of zero.
the number of domains correctly identified increased by ~ <dig> fold for the ld and md data sets and by ~ <dig> fold for the hd data set compared to individual sequence annotation .
the correlation between the annotation of the contigs alone and the theoretical domain content of the sample were higher than for all read annotation .
again it appears that the contigs capture the majority of the high-abundance transcripts and the unassembled debris capture the lower abundance transcripts , a combination of the two provides a stronger correlation with the known domain content of the samples than either individually .
clustering of the illumina assemblies  produced a net gain of between  <dig>  to  <dig>  extra domains, however this made only a relatively small improvement to the correlations with the known domain content of each sample .
discussion
the simulations show that the diversity of a metatranscriptome sample greatly impact the accuracy of protein domain annotations; with the high diversity simulations producing the weakest correlations with the known domain content of the sample.
with a highly diverse population of organisms and transcripts, the average coverage of each transcript will decrease, thus clustering will result in many small clusters and fewer transcripts will be sequenced to sufficient depth to allow extension into longer contiguous fragments.
however, regardless of the diversity level a better reflection of the domain content of the samples was achieved through applying data reduction techniques.
the largest improvements in the correlation with the known domain content of the samples was achieved through assembly  for the  <dig> simulations and assembly followed by clustering the contigs and debris together for the illumina simulations.
using near default parameters, highly homogeneous  contigs were recreated from both  <dig> and illumina data.
it has been noted previously that assembly of 'omics data is likely to favour highly abundant organisms
<cit> , and it therefore follows that it would also favour highly abundant transcripts.
the results of our simulations suggest that the annotations of contigs alone are insufficient, and we therefore recommend that they should be combined with those of the debris sequences to provide a better reflection of the real domain content of the samples.
overall, the simulated illumina samples produced stronger correlations with the known protein domain content than the dollar cost-equivalent amount of  <dig> sequence data.
while we attempted to perform this analysis as consistently as possible, it was necessary to employ different assembly programs for the  <dig> and illumina data – .
however, the overall pattern of correlations from the different methods is fairly consistent and it seems likely that the stronger correlations of the illumina simulations are due to the greatly increased coverage provided rather than any biases introduced by the methods.
while these simulations have their limitations, the results achieved were consistent with trials on real metatranscriptome data.
we applied the data reduction methods previously employed on simulated data to two real  <dig> metatranscriptome data sets: the mid-bloom, marine metatranscriptome from
<cit> ; and the 110 m marine metatranscriptome from an oxygen minimum zone
<cit> .
although the genuine domain content of a real microbial metatranscriptome is unknown, the results obtained from the gilbert and stewart metatranscriptomes were, in terms of data reduction and annotation rates, consistent with the medium and high diversity  <dig> simulations .
also, a recent study demonstrated that assembly of a simulated low diversity eukaryotic metatranscriptome could recreate a high number of contigs with low chimerism
<cit> .
in the future, these methods could be extended to exploit the increasing availability of microbial genomes and transcriptomes.
for example, in real metatranscriptome data, the most abundant transcripts are often associated with fundamental processes such as biosynthesis
<cit> .
as more microbial transcriptome data become available ), it should be possible to refine these models of transcript abundance to reflect increased levels of transcripts involved in core processes and thereby produce more realistic simulations of metatranscriptome data.
CONCLUSIONS
based on our simulations, it appears that older recommendations to omit the assembly stage when dealing with high-diversity samples do not extend to metatranscriptome data.
our results also show that including unassembled reads in downstream annotation can improve the overall accuracy and we would recommend that they should not be discarded after assembly.
therefore, whether dealing with  <dig> or illumina data, we recommend combining annotations from contigs and unassembled  sequences for  <dig> samples and employing a two-step data reduction of assembly followed by clustering of contigs and debris for illumina.
the high coverage afforded by illumina sequencing has made it an increasingly popular choice for sequencing microbial communities.
as more purpose built de-novo transcript assemblers become available there is a need for a systematic assessment of assembly tools and sequencing protocols for illumina metatranscriptome data.

 the sequencing and in-silico analysis of messenger rna is now routinely being applied to complex microbial communities in diverse eco-systems, including, but not limited to: soil, marine and intestinal habitats. the typical goals of metatranscriptomics are to taxonomically classify transcripts, predict their functions and quantify their abundances, and to relate these to environmental data in order to reveal how environmental conditions impact microbial communities in different habitats. metatranscriptome data sets typically consist of hundreds of thousands of sequences, or, more recently tens of millions of illumina sequences per sample. low taxonomic diversity and/or highly expressed genes can lead to a high degree of data redundancy; that is multiple identical or nearly identical sequence fragments. in an investigation into the proportion of artificial and natural duplicates in pyrosequenced metatranscriptome data, niu et al. reported that as much as 60% of all sequences in an early metatranscriptome data set were likely natural duplicates. therefore, some form of data reduction strategy is beneficial before running computationally intensive homology searches. two approaches that are commonly employed to reduce redundancy in large data sets are assembly: where sequences are assembled into longer contiguous fragments and clustering: sequences are grouped into clusters sharing a defined degree of similarity. the decisions as to whether to perform data reduction and which method to employ are influenced by several factors: the availability of reference genomes: if sufficient reference genomes are available for a small number of dominant species then the sequences can be mapped to them and taxonomy and function inferred and the relative abundance of the transcripts calculated. read length - are the unprocessed reads long enough to return annotations? current illumina platforms produce shorter reads than and are likely to return a lower hit rate to protein databases compared to longer reads. the diversity of the sample: although assembly can produce longer sequences and increase the accuracy of subsequent annotations, the variable coverage of transcripts in metatranscriptomics data sets and the presence of closely related organisms can lead to chimeric contigs. indeed, for highly diverse metagenomic samples it has been recommended that assembly not be performed at all. the aims of the analysis: if the read length is adequate for annotation and the intention is to count features then clustering at high identities is a recommended alternative. with the lower coverage but higher read length of metatranscriptome data, assembly is relatively uncommon and instead authors tend to either cluster or annotate sequences individually. clustering is regularly used for detecting and removing sequencing artifacts from data, grouping rrna data into operational taxonomic units, and grouping proteins into families. in addition to the known benefits of a reduction in the size of the data set and therefore computation time, we set out to assess whether, by clustering translated metatranscriptome sequences and transferring protein domain annotation from cluster representatives to cluster members - some of which may only partially cover protein domains used for classification, we can accurately increase the number of classifiable reads. more specifically, we investigated some popular data reduction tools and assessed their performance on simulated and illumina metatranscriptome data in terms of the accuracy of resulting protein annotations. note that although several approaches have been described to simulate metagenomic data sets and rna-seq data, to date only small scale attempts have been made to simulate metatranscriptome data sets based on a small number of species.  simulated data the simulated data sets contained sequences each, totalling ~ megabases of sequence per diversity level. between and 14% of sequences from each data set returned matches to pfam-a. when compared to the theoretical domain content, the correlation coefficients for all read annotation were, and for ld, md and hd respectively.table correlation coefficients between simulated data set annotations and known protein domain content all clustered contigs debris assembly 1 clustered assembly summary of pearson correlation coefficients between processed data sets and the known domain content of sample for low diversity, medium diversity and high diversity simulated and illumina metatranscriptomes. 1assembly includes annotation from both contigs and debris sequences. then, taking the parameter set that provided the largest increase in true positives minus false positives, compared to the annotation of all unclustered reads, we found that the best clustering parameters were: ≥ 60% overall similarity and 100% coverage of cluster member sequences for the ld data set; ≥80% similarity and 100% coverage of the cluster members for the md data set; and ≥60% similarity, ≥25% coverage of the cluster representative and between 0-50% minimum coverage of cluster members for the hd data set. while the best performing clustering parameters produced a net gain of between and domains, the correlation coefficients were slightly lower than for all read annotation ). the mira assemblies incorporated ~50% of all sequences into and contigs for the ld and md samples respectively, and ~30% of sequences into contigs for the hd sample. the average contig lengths were, and base pairs for ld, md and hd, respectively. the average contig entropy was, and for ld, md and hd respectively with %, % and % of contigs possessing an entropy of zero.figure contig entropy for assembled simulated metatranscriptomes. contig entropy plotted against contig length for a) mira assembled simulated data sets and b) trinity assembled simulated illumina data sets. plots represent, from left to right: low diversity, medium diversity and high diversity data sets. for the ld and md data sets, the net gain of true positives was a ~100% increase, and for the hd data set an increase of ~20% was achieved. the contigs alone had a weaker correlation with the theoretical domain content than all read or clustered read annotation. when combined with the debris sequences, the correlation coefficients for all three samples were higher than for all all-read or clustered annotations ). this could be due to two factors: firstly the low proportion of sequences incorporated into the contigs, and secondly the assemblies may be biased towards high-abundance transcripts.figure results from pfam-a annotated simulated metatranscriptomes. percentage of true positives, false positives, true negatives and potential domains based on a comparison with the known domain content of the data sets for all reads, best clustering, assembly and clustered assembly. a) results for simulated data sets, from left to right: low, medium and high diversity. b) results for simulated illumina data sets from left to right: low, medium and high diversity. correlation between high diversity simulations and known protein domain content. correlation plots of pfam-a annotations of each processed data set compared to known domain content for a) high diversity simulated data set and b) high diversity illumina simulated data set. top row, left to right: all reads unprocessed; clustered reads; assembly - contigs only. bottom row, left to right: assembly – debris only; assembly – contigs and debris combined; clustered assembly. pearson correlation coefficient shown in top left corner. clustering of the assemblies led to a very slight increase in the detection of true positives but the overall effect was a very slight reduction in the correlation with the theoretical domain content compared to the unclustered assembly. simulated illumina data around 4% of the illumina reads could be annotated with pfam-a domains. the correlation coefficients for all read annotation with the theoretical domain content were. the illumina data sets were clustered with the best performing parameter set for the equivalent diversity level identified in the simulations described above. while clustering reduced the data sets by ~40% for ld and md and ~25% for the hd data set the resulting annotations had a weaker correlation to the theoretical domain content of the sample.the trinity assemblies incorporated ~40% of sequences from the ld and md data sets into and contigs respectively with an average length of ~400 nt. for the hd data set, ~14% of reads from the hd data set into contigs with an average length of 328 nt. the average contig entropy was, and for ld, md and hd respectively with %, % and 92% of contigs possessing an entropy of zero. the number of domains correctly identified increased by ~ fold for the ld and md data sets and by ~ fold for the hd data set compared to individual sequence annotation. the correlation between the annotation of the contigs alone and the theoretical domain content of the sample were higher than for all read annotation. again it appears that the contigs capture the majority of the high-abundance transcripts and the unassembled debris capture the lower abundance transcripts, a combination of the two provides a stronger correlation with the known domain content of the samples than either individually. clustering of the illumina assemblies produced a net gain of between to extra domains, however this made only a relatively small improvement to the correlations with the known domain content of each sample. discussion the simulations show that the diversity of a metatranscriptome sample greatly impact the accuracy of protein domain annotations; with the high diversity simulations producing the weakest correlations with the known domain content of the sample. with a highly diverse population of organisms and transcripts, the average coverage of each transcript will decrease, thus clustering will result in many small clusters and fewer transcripts will be sequenced to sufficient depth to allow extension into longer contiguous fragments. however, regardless of the diversity level a better reflection of the domain content of the samples was achieved through applying data reduction techniques. the largest improvements in the correlation with the known domain content of the samples was achieved through assembly for the simulations and assembly followed by clustering the contigs and debris together for the illumina simulations. using near default parameters, highly homogeneous contigs were recreated from both and illumina data. it has been noted previously that assembly of'omics data is likely to favour highly abundant organisms, and it therefore follows that it would also favour highly abundant transcripts. the results of our simulations suggest that the annotations of contigs alone are insufficient, and we therefore recommend that they should be combined with those of the debris sequences to provide a better reflection of the real domain content of the samples. overall, the simulated illumina samples produced stronger correlations with the known protein domain content than the dollar cost-equivalent amount of sequence data. while we attempted to perform this analysis as consistently as possible, it was necessary to employ different assembly programs for the and illumina data –. however, the overall pattern of correlations from the different methods is fairly consistent and it seems likely that the stronger correlations of the illumina simulations are due to the greatly increased coverage provided rather than any biases introduced by the methods. while these simulations have their limitations, the results achieved were consistent with trials on real metatranscriptome data. we applied the data reduction methods previously employed on simulated data to two real metatranscriptome data sets: the mid-bloom, marine metatranscriptome from; and the 110 m marine metatranscriptome from an oxygen minimum zone. although the genuine domain content of a real microbial metatranscriptome is unknown, the results obtained from the gilbert and stewart metatranscriptomes were, in terms of data reduction and annotation rates, consistent with the medium and high diversity simulations. also, a recent study demonstrated that assembly of a simulated low diversity eukaryotic metatranscriptome could recreate a high number of contigs with low chimerism. in the future, these methods could be extended to exploit the increasing availability of microbial genomes and transcriptomes. for example, in real metatranscriptome data, the most abundant transcripts are often associated with fundamental processes such as biosynthesis. as more microbial transcriptome data become available ), it should be possible to refine these models of transcript abundance to reflect increased levels of transcripts involved in core processes and thereby produce more realistic simulations of metatranscriptome data.  based on our simulations, it appears that older recommendations to omit the assembly stage when dealing with high-diversity samples do not extend to metatranscriptome data. our results also show that including unassembled reads in downstream annotation can improve the overall accuracy and we would recommend that they should not be discarded after assembly. therefore, whether dealing with or illumina data, we recommend combining annotations from contigs and unassembled sequences for samples and employing a two-step data reduction of assembly followed by clustering of contigs and debris for illumina. the high coverage afforded by illumina sequencing has made it an increasingly popular choice for sequencing microbial communities. as more purpose built de-novo transcript assemblers become available there is a need for a systematic assessment of assembly tools and sequencing protocols for illumina metatranscriptome data. 
for the illumina simulations, we found that a two-step approach of assembly followed by clustering of contigs and unassembled sequences produced the most accurate reflection of the real protein domain content of the sample. however, it is unclear how such data reduction techniques impact the annotation of metatranscriptome data and thus affect the interpretation of the results. to investigate the effect of such techniques on the annotation of metatranscriptome data we assess two commonly employed methods: clustering and de-novo assembly. to do this, we also developed an approach to simulate and illumina metatranscriptome data sets with varying degrees of taxonomic diversity. based on these data we recommend that assembly be attempted, and that unassembled reads be included in the final annotation for metatranscriptome data, even from highly diverse environments as the resulting annotations should lead to a more accurate reflection of the transcriptional behaviour of the microbial population under investigation. for metagenomic data, it has been observed that the variable coverage and presence of closely related organisms can lead to fragmented assemblies containing chimeric contigs that may reduce the accuracy of downstream analyses and some advocate the use of alternate data reduction techniques. metatranscriptome sequence data can contain highly redundant sequences from diverse populations of microbes and so data reduction techniques are often applied before taxonomic and functional annotation.
 metatranscriptome sequence data can contain highly redundant sequences from diverse populations of microbes and so data reduction techniques are often applied before taxonomic and functional annotation. for metagenomic data, it has been observed that the variable coverage and presence of closely related organisms can lead to fragmented assemblies containing chimeric contigs that may reduce the accuracy of downstream analyses and some advocate the use of alternate data reduction techniques. however, it is unclear how such data reduction techniques impact the annotation of metatranscriptome data and thus affect the interpretation of the results.  to investigate the effect of such techniques on the annotation of metatranscriptome data we assess two commonly employed methods: clustering and de-novo assembly. to do this, we also developed an approach to simulate and illumina metatranscriptome data sets with varying degrees of taxonomic diversity. for the illumina simulations, we found that a two-step approach of assembly followed by clustering of contigs and unassembled sequences produced the most accurate reflection of the real protein domain content of the sample. for the simulations, the combined annotation of contigs and unassembled reads produced the most accurate protein domain annotations.  based on these data we recommend that assembly be attempted, and that unassembled reads be included in the final annotation for metatranscriptome data, even from highly diverse environments as the resulting annotations should lead to a more accurate reflection of the transcriptional behaviour of the microbial population under investigation. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords metatranscriptomicssequence processingdata reductionclusteringassemblyissue-copyright-statement© the author 2014 