DATMA: Distributed AuTomatic Metagenomic Assembly and annotation framework

DATMA stages

Workflow design

DATMA is a command line application written in Python and tested in Linux. We provide an installation script in our GitHub to automatically install DATMA source codes and the tools that make up part of the workflow. We tested it on Ubuntu 16.04 and included a user manual for custom compilation and installation of source codes on other Linux distributions. By default, DATMA configures all tools called in the workflow according to the authors recommended parameters, but these values can be modified using a configuration file. In this file, the user specifies the input sequence file, the output directory, the workflow stages, the database directories, the number of threads to use, CLAMEs parameters, etc. The minimum configuration file should contain the input-sequence file, the sequence type (i.e., FASTA, FASTQ, or SFF) and the output directory. We show a complete configuration file in DATMA’s user manual.

Although there are several workflow engines (e.g., Snakemake (Koster & Rahmann, 2012), Nextflow (Di Tommaso et al., 2017), Ibis (Bal et al., 2010), and Swift (Wilde et al., 2011)) we selected COMPSs (Badia et al., 2015) which provided us with the tools that we required: simple python interface and automatic parallel task distribution and synchronization. COMPSs offers a simple programming model, that does not require the use of APIs to modify the original user applications, and enables the execution of the same code on different back-ends. It uses a sequential description of the work, and it identifies and launches asynchronous parallel tasks automatically. A complete description of COMPs and its performance is in Badia et al. (2015).

COMPSs allows DATMA to be executed in single or distributed mode. In single mode, the framework executes all the stages into the same computer. In distributed mode, DATMA uses a master-worker execution strategy, to distribute application tasks across the different computer nodes available. It executes the quality control, 16S rRNA identification, and CLAME binning stages in the master node (these stages can be multi-threaded). Once the bins are generated, DATMA assembles and annotates them using the available nodes. It requires two configuration files (resources.xml and project.xml) within the execution environment. The first file contains the information of the available computing resources, and the second file has information about the computing resources to be used for a specific execution. The user manual has an example of each file.

CAMI dataset

We used the first Critical Assessment of Metagenome Interpretation (CAMI) (Sczyrba et al., 2017) challenges to evaluate DATMA performance. CAMI consortium provides three metagenome datasets at different complexity levels (high, medium, and low complexity). CAMI_low consists of one simulated Illumina HiSeq data, with size 15 Gbp and a total of 40 genomes and 20 circular elements. CAMI_medium includes two samples with a total size of 40 Gbp (132 genomes and 100 circular elements). CAMI_high is a simulated time series benchmark dataset with five samples of size 75 Gbp (596 genomes and 478 circular elements). All the samples in the three datasets are paired-end 150-bp Illumina reads and are available at CAMI web site (Sczyrba et al., 2017). CAMI datasets have been studied by several binning tools (i.e., CONCOCT (Alneberg et al., 2014), MyCC (Lin & Liao, 2016), COCACOLA (Lu et al., 2017), BinSanity (Graham, Heidelberg & Tully, 2017), MaxBin2 (Wu, Simmons & Singer, 2016), and MetaBat2 (Kang et al., 2019)). In their last report MetaBAT2 shows better performance that the other binning tools in all CAMI experiments. Therefore, we use MetaBAT2 to compare our results.

We downloaded the three synthetic datasets from the CAMI website and studied each individually. For the three CAMI experiments, we configured DATMA with default parameters. It starts by removing low-quality reads (quality Q < 30 and length <60 bp). Then, the remaining sequences were merged using the FLASH2 (Magoč & Salzberg, 2011) and forced to combine the reads that were not merged (only for binning purposes) using extra “NNN” because they were too large (see methods section). Then, the 16S rRNA ribosomal sequences were separated using BWA (Li & Durbin, 2009) to map the reads against the Rfam database (Griffiths-Jones et al., 2003). The remaining reads were binned with CLAME using b=70bp and iterating with b=50bp and b=30bp. We set DATMA to report only bins with more than 100.000 reads and selected SPAdes (Nurk et al., 2013) as the assembler tool. We used CheckM results, reported by DATMA, to assess the genome completeness and compare those results against the MetaBAT2 report. For all the experiments, we set the number of threads to twenty and configured DATMA to process blocks of 20 million reads. The complete configuration file, for each dataset, is available in DATMA’s GitHub.

Tables S2, S3 and S3 in Additional File 1, contains the full report generated by CheckM for the three CAMI datasets. Figure 2 summarizes the best bins (completeness >40% and low contamination) from each experiment. In the CAMI-low Complexity dataset, DATMA recovers 27 of 40 genomes present in the sample, with completeness higher than 60%. In the CAMI-medium experiment, it reports 28/132 genomes, most of them with completeness higher than 60%. In the CAMI_high dataset, DATMA reports 33/596 genomes, all of them with completeness higher than 40%. Figure 2 also compares our results against the MetaBAT2 report. It shows that MetaBAT2 has better performance than DATMA in all the experiments; however, it is essential to notice that DATMA uses the raw reads, while MetaBAT2 uses the golden bins (a set of well-defined assembly contigs provides by CAMI). These results suggest that DATMA can recover the predominant genomes from a complex metagenome of a microbial community. However, they also indicate the limitation of our tool to recover species in less abundance into the metagenome (only the most abundant were reported) or separate close-taxonomic genomes (high contamination level is present in some bins).

Brocadia caroliniensis metagenome

We used a metagenome recovered from a full-scale glycerol-fed nitritation-denitritation separate centrate treatment process (NCBI project PRJNA228949). The original paper (Park et al., 2017) reports that 2,448,982 reads were manually analyzed to generate 209 contigs (with size >500 bp) that integrate the draft genome for Brocadia caroliniensis species. We downloaded the raw reads and analyzed them with DATMA. It was executed with default parameters to remove low-quality bases and reads that were too short (Q<30 and length <70 bp). The 1,860,653 leftover reads were aligned against the Rfam database (Griffiths-Jones et al., 2003) to remove 16S rRNA gene sequences. After removing 12,754 reads, DATMA called CLAME (Benavides et al., 2018) with 1,847,899 sequences using b=70bp, as the number of bases alignment parameter. The bins with more than 2000 reads were assembled with SPAdes (Nurk et al., 2013).

Table 1 summarizes the number of bins generated, the assembly metrics, the total ORFs detected, the completeness-contamination of the bins, and the computational time used by DATMA. It also contrasts these results against the report produced by MetaWRAP (Uritskiy, DiRuggiero & Taylor, 2018) and SqueezeMeta (Tamames & Puente-Sánchez, 2019) frameworks. MetaWRAP completeness of the Brocadia genome is higher than the obtained by DATMA; but, DATMA obtains a better N50. SqueezeMeta annotated most reads as Brocadicae family, but it generated a larger number of contigs than the other frameworks. DATMA was the fastest tool.

We configured DATMA to use Quast tool (Gurevich et al., 2013) with a reference. It allowed evaluating the coverage and depth of the bins generated by our pipeline. Table 2 shows that DATMA assembly covers about 97% of the Brocadia genome, and only six contigs did not align with precision to the reference sequence. DATMA covered more of the genome than the other tools, but it presented more unaligned contigs than MetaWRAP. DATMA had similar results than the manual process of the original paper but it was automatic and faster than the other tools.

San Fernando biosolid metagenome

We used DATMA to study the biosolid metagenome produced by the San Fernando wastewater plant located in Medellin-Colombia. Two biosolid samples (each about 0.5 kg) were collected and transferred to the laboratory in refrigeration. The DNA extraction was done using PowerMax^® Soil DNA Isolation Kit supplied by MOBIO Corporation (Diagnostics Products MP Biomedicals, 2019). The samples were then sequenced using ROCHE’s 454 Titanium technology in 3/4 PTP at the Centro Nacional de Secuenciación Genómica-CNSG, Universidad de Antioquia, Medellin, Colombia. A total of 6,206,317 reads were analyzed. A study of the microbial diversity, as well as the methanogenesis pathway of this metagenome, is presented in Bedoya et al. (2019).

DATMA was executed with default parameters to remove low-quality sequences (Q<30 and length <70 bases) and 5,668,260 reads were left. These reads were aligned against the Rfam database (Griffiths-Jones et al., 2003) to identify 16S rRNA ribosomal sequences. A total of 53,557 reads were detected and separated by DATMA. The 5,614,703 leftover sequences were binned with CLAME using default parameters but reporting bins with more than 5000 reads. We selected SPAdes (Nurk et al., 2013) as the assembler tool. We compared DATMA’s results and performance against MetaWRAP (Uritskiy, DiRuggiero & Taylor, 2018) and SqueezeMeta (Tamames & Puente-Sánchez, 2019) frameworks. Moreover, we used MG-RAST server (Wilke et al., 2016) to analyze this metagenome. We set the number of threads to four for all pipelines. We used Squeezeemeta in co-assembly mode using two samples collected in different seasons of the year. However, it generated a No-Consensus output in the merge stage. Therefore, we configured it in sequential mode and executed it with the complete biosolid dataset.

We started rating the species abundance in the metagenome using the 16S rRNA sequences reported by DATMA. Figure 3 shows the taxonomic annotation generated by the RDP Classifier tool (Wang et al., 2007) with the ribosomal reads. It indicates that Bacteria, with 70% of the reads, is the primary domain. Within this domain, Firmicutes, Proteobacteria, Bacteroidetes, and Chloroflexi are the main phyla (with 16%, 14%, 13%, and 3% respectively). Each of them contains several families, except for Chloroflexi, in which case Anaerolineaceae is the dominant family (it corresponds to the 3% of all metagenome).

Table 3 summarizes the following results obtained by DATMA, MetaWRAP (Uritskiy, DiRuggiero & Taylor, 2018) and SqueezeMeta (Tamames & Puente-Sánchez, 2019): number of contigs, their metrics, the ORF detected, the number of genomes reported and the quality according to the CheckM tool (Parks et al., 2015). We have included the report generated by the MG-RAST server (Wilke et al., 2016). MetaWRAP reports eight genomes with completeness higher than 80% and a contamination level less than 7%. However, most of them cannot be assigned with precision into a clade family, and most importantly, no bin was annotated into Anaerolineaceae family (see Fig. S4, in Additional File 3). SqueezeMeta shows that Proteobacteria is the dominant phylum, but no bins belonging to Chlorofexi were reported. MG-RAST indicates that most of the reads are classified as Pseudomonadaceae and Anaerolineaceae families, but because we submitted the data as a private project, no additional information could be collected. We observed that DATMA is the only tool which reports two primary bins annotated in the Anaerolineaceae family. Moreover, it was the fastest tool.

According to MIMAG standards (Bowers et al., 2017) to report a genome, Bin0 with 1292 contigs and Bin1 with 647 contigs have suitable results to propose a draft genome. However, the contamination level indicated in Bin0 is too high to propose a draft genome. We focus our study on Bin 1.

Figure 4 shows DATMA’s annotation report using the Kaiju tool (Menzel, Ng & Krogh, 2016) for Bin 1. It indicates that most of the contigs were annotated into the Chloroflexi phylum and Anaerolineaceae family. Moreover, the relation between the number of ORFs and the genome estimation (1 ORF per Kbp) agrees with the relation reported for similar species from this family (i.e., Pelolinea submarina with 3131 ORFs, 3.5 Mbp and a relation of 0.89 ORFs/Kbp and Leptolinea tardivitalis with 3301 ORFs, 3.69 Mbp and a relation of 0.90 ORFs/Kbps). We mapped the 53,557 reads into the 647 contigs and computed the assembly depth. It presents an average of 18.94 reads per position.

We enriched the contigs with the 16S rRNA sequences removed in the 16S-identification stage. We manually selected the reads annotated within the Anaerolineaceae family and assembled them using the SPAdes (Nurk et al., 2013) to obtain the complete 16S rRNA gene. Analysis of this 16S rRNA ribosomal gene, using BLAST (against the local NT) indicated that the 16S rRNA gene is related to Anaerolinea thermophila UNI-1 DNA species (NCBI Accession AP012029). To improve the taxonomic annotation, we used MEGA 7.0 (Kumar, Stecher & Tamura, 2016) to build a phylogenetic tree using the 16S rRNA sequences and the Ribosomal data project database (Cole et al., 2014). The evolutionary tree, in Fig. 5, was inferred by using the Maximum Likelihood method with the Jukes-Cantor model (Kumar, Stecher & Tamura, 2016) and the process described by Brumm et al. (2015). It indicates that the recovered reads are close to the family Anaerolineaceae and are related to the genus Pelolinea and Leptolinea.

Using the set of standards for the minimum information regarding a metagenome-assembled genome (MIMAG) proposed by Bowers et al. (2017) and the previous results, we can use the contigs from Bin 1 to describe a Low-quality draft genome that belongs to the family Anaerolineaceae, which is closely related to the genus Pelolinea and Leptolinea. We called this draft genome Anaerolineaceae_UdeA_SF1. We submitted the draft genome (assembled contigs and respective reads) to the NCBI database (Bioproject PRJNA529916).

Computational performance

To illustrate the computational performance of DATMA we executed the experiments within two different scenarios: (i) single mode, using only the Master machine, and (ii) distributed mode, using the Master machine with multiple workers like a grid of computers. We simulated the grid of computers using tree servers (Master, Worker1, and Worker2) connected via a secure shell connection. Table S5, in Additional File 1, illustrates the computer specifications of each server. To simulate a more significant number of workers, like a bigger grid of computing, we allow for several tasks to run on the same computer. Applications were configured to use four threads on all the experiments.

Figure 6 and Table S6, in Additional File 1, show the execution time for all the datasets using several scenarios. It shows that computational time decreases as the number of workers increase. Figure 6 also illustrates the memory performance of DATMA. It reports a peak in the binning stage, but it then decreases when DATMA distributes the next tasks into the available computing resources. Because we used blocks of the same size (20 million reads) to bin the raw reads of the CAMI datasets, the peak of memory is the same for all the experiments. It can be reduced by modifying the blocks size parameter, in the configuration file.