Additional file 1 of Reverse engineering environmental metatranscriptomes clarifies best practices for eukaryotic assembly
Creators: name: Krinos, Arianna I.
name: Cohen, Natalie R.
name: Follows, Michael J.
name: Alexander, Harriet
Publisher: name: figshare
Year: 2024
Subjects: subject: Microbiology
subject: FOS: Biological sciences
subject: Genetics

description: Additional file 1: Fig. S1. Visual summary of the effect of clustering with mmseqs2 on recovered enera, KOs (functional annotations) and file size in bytes. While reducing the equence identity threshold for clustering results in up to a 30% average reduction in file size, he number of recovered genera and functional annotations are only modestly impacted, specially at high coverage. An intermediate sequence identity of 0.8 and coverage of .8 would result in a 15-25% average reduction in file size, but leave distinct functional and taxonomic annotations unchanged. Fig. S2. Main text figure facetted by simulated METSP assembly group (two different sets of organisms). Of note is that the eukrhythmic reassembly accurately recapitulates the bimodal distribution in GC-content observed in MMETSP group A’s designer metatranscriptomic sequences. Fig. S3. Protein sequence lengths in the reassemblies as compared to the designer. The 1-to-1 line shows where sequences would fall if the average length of recovered protein sequences via TransDecoder were identical between the designer assemblies and the eukrhythmic-derived reassembled products; the fact that all samples fall in the lower right half of the plot indicates that protein sequences were consistently larger in the designer assemblies as compared to the eukrhythmic reassembled products. Fig. S4. Mean contig length as a function of the number of assemblers that found a sequence that matched the given description. For the 4-assembler cluster, this means that all four assemblers tested identified a sequence that matched the sequence included in the distribution when clustered within eukrhythmic. Panel A corresponds to MMETSP Group A while Panel B corresponds to MMETSP Group B. Welch’s independent T-tests and Kolmogorov-Smirnoff tests for between-distribution goodness of fit computed on these length distributions reveals that the overall distribution of lengths for 1 vs. 2 vs. 3 vs. 4-assembler distributions are all statistically significantly different (p &lt; 1e−6), with greater assemblers within a cluster leading to higheraverage length. Fig. S5. Salmon mapping percentages against the simulated raw reads when computed individually against each of the four assembly tools used by eukrhythmic. rnaSPAdes consistently outperformed the other assemblers with regard to percentage mapping, average length, and number of annotations. Fig. S6. Main text figure, but with A: TPM comparison for the designer assemblies and the reassembled products from euk rhythmic labelled by distinct groupings of simulations relative to their subset of MMETSP organisms; each point is filled according to its “MMETSP group”. Fig. S7. Main text figure with the results split by “MMETSP group” to demonstrate that different taxonomic groupings of organisms included in the simulation does not impact the general trends observed in the results. All three panels are divided by samples using the two “MMETSP groups” of individual organism transcriptomes. Fig. S8. Salmon percentage mapping by MMETSP group when contig products from all assemblers as in eukrhythmic are used (bottom distribution with one or more assembler) as compared to when only contigs agreed upon by multiple assemblers are used. Average percentage mapping decreases progressively with fewer contigs being included as the criteria for inclusion are made stricter. Fig. S9. Proportion of contigs from each clustering subset that had ORFs extracted from the sequence. The vast majority of contigs had a single predicted ORF, but rnaSPAdes alone had the greatest number of nucleotide contigs on which an ORF could not be detected. In practice, these sequences might be assumed to be non-coding. Fig. S10. Proportion of contigs from each clustering subset that had ORFs extracted from the sequence across all contigs identified by the assembler. rnaSPAdes had a higher number of contigs without an identified ORF. Fig. S11. Proportion of contigs from the designer assembly that had ORFs extracted from the sequence. In general, fewer contigs did not have an identified ORF, and a higher number of contigs had multiple predicted ORFs than in the reassembled products from euk rhythmic. Fig. S12. Comparison of the abundance of KO IDs within functional annotations across Narragansett Bay samples and different combinations of metatranscriptome assemblers. A black dotted line indicates a one-to-one relationship, meaning that the abundance of KOs that fall along this line are exactly as abundant using the assembler listed on the x-axis and using the assembler listed on the y-axis. On the top left, Trinity is compared to MEGAHIT, on the top right Trinity is compared to rnaSPAdes, on the bottom left MEGAHIT is compared to rnaSPAdes and on the bottom right MEGAHIT is compared to SPAdes. Each point corresponds to a single KO within a sample. While rnaSPAdes tended to report high abundances of each identified KO relative to the other assemblers, MEGAHIT reported fewer instances of each KO than the other three assemblers in most samples. This may be due to the approaches adopted by the two assemblers. Whereas in a typical assembly, k-mers that appear only once are assumed to be the result of error, these k-mers may represent real and important diversity in a set of low-abundance whole-community sequences [47]. The MEGAHIT assembler is an example of metagenomics-specific software that defines “mercy k-mers” which come into play in between two k-mers within a single read that are sequenced more than once. rnaSPAdes, as an example, does not employ a “mercy” strategy but instead decreases the coverage threshold significantly in comparison to genomic assembly [36]. One is adopting the meta-omic assumption that exceedingly low coverage is plausibly not artificial, while the other is being more generous with respect to coverage and memory use, but intending merely to minimize the influence of sequencing error. The strategy adopted by MEGAHIT may result in incomplete consideration of isoforms, which could have contributed to the relatively low recovery rate of multiple copies of the identified KO groups by the MEGAHIT assemblies of the Narragansett Bay metatranscriptome samples. Fig. S13. Proportion of total normalized TPM value of major taxonomic groups according to the four metatranscriptomic assemblers that were tested. rnaSPAdes had a lower proportion of TPM assigned to Ochrophyta (including diatoms), but further investigation appeared this to be in large part a consequence of the large number of small contigs produced by rnaSPAdes. Fig. S14. Proportion of total normalized TPM as reported by Salmon from the raw reads assigned to each taxonomic category in the assemblies produced by each of the four assemblers (facets). In particular, in sample S2 the taxonomic breakdown of the eukaryotic community differs importantly in the midst of a bloom of the diatom Skeletonema. Fig. S15. Proportion of assigned raw reads as reported by Salmon assigned to each taxonomic category for the four metatranscriptome assemblers, expressed as a proportion of the total. Fig. S16. Proportion of assigned raw reads as reported by Salmon assigned to each taxonomic category for the four metatranscriptome assemblers, expressed as a proportion of the total. Fig. S17. Average contig length by annotation of contigs generated by the four assembly tools across samples. Error bars show the standard error of the mean. Fig. S18. Total number of contigs generated by the four tested assemblers for each of the taxonomic groupings considered. rnaSPAdes tended to produce more contigs than the other four assemblers, but these contigs were often shorter and occasionally led to misleading community composition results. Fig. S19. Analogous stacked abundance plots by taxonomic grouping to the relative abundance plots presented in Alexander et al. [23] for the five in-situ sampling points collected from Narragansett Bay. Fig. S20. Salmon percentage mapping of coding sequences (left) vs. entire contigs (right) for the Tara Oceans samples. Mapping only to coding sequences from the assembly decreased the mean percentage mapped, as reported by Salmon. Fig. S21. Per-sample taxonomic annotations for all euk rhythmic CDSs, including those which were and were not found to have a match to the MATOU database. Fig. S22. Total TPM assigned per sample to sequences that were assigned a EUKulele annotation, but did not have a significant blast match to the MATOU database from Carradec et al. [22]. Fig. S23. Full length distribution of coding sequences recovered by the eukrhythmic assembly and were (left) not found in the MATOU database [22] vs. (right) found in the MATOU database. Fig. S24. Per-sample average length of coding sequences recovered by the eukrhythmic assembly and were (left) not found in the MATOU database [22] vs. (right) found in the MATOU database. Table S1. (Supplement) Effect of clustering designer assembly on assembly size and annotations. Clustering was performed on the original “designer metatranscriptome” set of contigs from the MMETSP references using the mmseqs2 tool (Mirdata et al. 2019). eukrhythmic uses a coverage level of 0.98 and sequence identity of 1 for mmseqs2 clustering. See Supplementary Figure 1 for a graphical summary of the influence of sequence identity and coverage on the size of the recovered assembly and its functional and taxonomic annotations