BlobTools: Interrogation of genome assemblies

Implementation

BlobTools is written in Python and consists of a main executable that allows the user to interact with the implemented modules (see Table 1). It offers a simple, modular command line interface which can easily be adapted to process multiple datasets simultaneously using GNU parallel (Tange, 2011). Inputs for BlobTools are standard file formats commonly created during the course of genome assembly projects. The primary processing in BlobTools constructs a BlobDB data structure based on user input. From this data structure, BlobTools generates easily interpretable, two-dimensional visualisations ready for publication, in conjunction with tabular output, enabling the user to partition sequences and paired-end (PE) reads contributing to them, for separate downstream processing. We present two recommended workflows, one targeted at de novo genome assembly projects in the absence of a reference genome (Figure 1A) and another for projects where a reference genome is available (Figure 1B).

Figure 1. Two common BlobTools workflows for taxonomic interrogation of paired-end (PE) read datasets.

(A) Workflow A. Targeted at de novo genome assembly projects in the absence of a reference genome. 1: Creation of a BlobDB data structure based on input files. 2: Visualisation of assembly and generation of tabular output. 3: Partitioning of sequence IDs in assembly, based on user-defined parameters informed by the visualisations. 4: Partitioning of PE reads based on sequence IDs. (B) Workflow B. Targeted at projects where a reference genome is available. 1: Reads are mapped against the reference genome. 2: BAM file is processed to generate FASTQ files based on read mapping behaviour. 3: FASTQ file of read pairs where neither read maps to the reference genome (UnUn) are assembled de novo and used in workflow A. 4: partition of read pairs of target taxon recovered from workflow A are assembled together with the other target taxon read pairs from step 2 and used in workflow A.

Taxonomy assignment

Taxonomy assignment in BlobTools is based on user-supplied, tab-separated-value (TSV) files composed of three columns: the input sequence ID, a NCBI TaxID, and a numerical score. We refer to these TSV files as ‘hits’ files below. They can be generated from the output of sequence similarity searches, such as BLAST (Camacho et al., 2009) or Diamond blastx (Buchfink et al., 2015) searches against public or reference databases, or the output of other contaminant identification tools. The BlobTools module taxify allows easy conversion of tabular file formats to BlobTools compatible input, in addition to annotation of similarity search results based on NCBI TaxID mapping files, as available from UniProt and NCBI.

Based on these inputs, BlobTools assigns a single NCBI taxonomy for each sequence in the assembly, based on the highest scoring NCBI TaxID at the following taxonomic ranks: species, genus, family, order, phylum, and superkingdom. Score calculation can be controlled by the user through a minimal score threshold (--min_score) and a minimal difference in scores (--min_diff) between the best and second-best scoring taxonomy. In addition, three non-canonical taxonomic annotations are possible: ‘no-hit’, the suffix ‘-undef’ and ‘unresolved’. Sequences not assigned to any taxonomic group, or not passing the --min_score threshold, are labelled ‘no-hit’. If a NCBI TaxID has no explicit parent at a taxonomic rank, the suffix ‘-undef’ is appended to the next upper taxonomic rank for which one does exist. In cases where the score difference between the best and second-best hits is smaller than --min_diff, sequences are labelled as ‘unresolved’.

Multiple ‘hits’ files can be provided as input. In this case, the behaviour of the taxonomy assignment process can be controlled further through ‘taxrules’. The highest scoring taxonomy can either be inferred across all files (‘bestsum’) or successively (‘bestsumorder’) in the order they were supplied as input, allowing only sequence that received no hits from one file to be considered for taxonomic annotation in the next file, thereby leveraging reliability of scores of different input file sources.

The original blobology pipeline (Kumar et al., 2013) recommended the use of a single, best BLAST hit per sequence for taxonomy assignment. However, taxonomically mis-annotated sequences in databases (derived from inclusion of un-screened genome assemblies) can lead to erroneous taxonomic annotation. BlobTools mitigates this issue by accepting multiple hits per sequence and allocating taxonomy based on the highest sum of scores.

It should be noted that a definitive taxonomic placement for every sequence in the assembly is not required for successful taxonomic partitioning of sequences, since differential coverage and sequence composition profiles between the genomes are often sufficient.

Visualisations

In BlobTools, sequences are depicted as circles in BlobPlots (as opposed to dots in the blobology pipeline), with diameters proportional to sequence length. The scatter-plot is decorated with coverage and GC histograms for each taxonomic group, which are weighted by the total span (cumulative length) of sequences occupying each bin. A legend reflects the taxonomic affiliation of sequences and lists count, total span and N50 by taxonomic group. Taxonomic groups can be plotted at any taxonomic rank and colours are selected dynamically from a colour map. The number of taxonomic groups to be plotted can be controlled (--plotgroups, default is ‘7’) and remaining groups are binned into the category ‘others’. An example is shown in Figure 2A.

Figure 2. Visualisations of the combined assembly of simulated sequencing libraries.

(A) BlobPlot of the assembly. Sequences in the assembly are depicted as circles, with diameter scaled proportional to sequence length and coloured by taxonomic annotation (at the rank of ’order’) based on BLASTn and Diamond blastx similarity search results provided in this order and using taxrule ’bestsumorder’. Circles are positioned on the X-axis based on their GC proportion and on the Y-axis based on the sum of coverage across both library A and library B. (B) ReadCovPlot of library A. (C) ReadCovPlot of library B. In ReadCovPlots, mapped reads are shown by taxonomic group at the rank of ’order’.

The power of differential coverage profiles across different sequencing libraries for partitioning sequences in an assembly prompted the development of CovPlots (Figure 3) (Koutsovoulos et al., 2016), which are analogous to BlobPlots, except that the GC-axis is substituted by the coverage-axis from another sequencing library. CovPlots can be used for the visualisation of patterns of differential coverage signatures between taxonomic groups in the assembly.

Figure 3. CovPlot of the combined assembly of simulated sequencing libraries.

Sequences in the assembly are depicted as circles, with diameter scaled proportional to sequence length and coloured by taxonomic annotation (at the rank of ’order’) based on BLASTn and Diamond blastx similarity search results provided in this order and using taxrule ’bestsumorder’. Circles are positioned on the X-axis based on coverage in library A and on the Y-axis based on coverage in library B. Parameters for partitioning the sequences in the assembly (which were applied to the tabular representation of the BlobDB) are indicated as dotted grey lines and text annotations in the scatter plot.

The modules for generating BlobPlots and CovPlots support additional input parameters controlling visualisation behaviour, including cumulative addition (--cumulative) or generation of separate plots for each taxonomic group (--multiplot), exclusion (--exclude) or relabelling (--relabel) of taxonomic groups, assignment of specific HEX colours to groups (--colour) or labelling sequences based on arbitrary, user defined categories (--catcolour). The latter could be, for instance, binned categories of RNAseq mappings to sequences in the assembly as shown in Koutsovoulos et al. (2016).

ReadCovPlots (Figure 2B and 2C) visualise the proportion of reads of a library that are unmapped or mapped, showing the percentage of mapped reads by taxonomic group, as barcharts. These can be of use for rapid taxonomic screening of multiple sequencing libraries within a single project. The underlying data of ReadCovPlots and additional metrics are written to tabular text files for custom analyses by the user.

Support of multiple coverage libraries

BlobTools supports coverage input (BAM/CAS format) from multiple sequencing libraries. As these data formats contain more information than needed, BlobTools parses coverage information of sequences (normalised base coverage and read coverage) into COV files in TSV format. These files can be generated through the module map2cov prior to construction of a BlobDB.

Within the BlobDB data structure, base and read coverage information is stored for each sequence in the assembly. If more than one coverage file is supplied, BlobTools constructs an additional coverage library (‘cov_sum’) internally, containing the sum of coverages for each sequence across all coverage files. This internal coverage library is considered when extracting views or plotting visualisations.

Operation

System requirements for BlobTools include a UNIX based operating system, Python 2.7, and pip. An installation script is provided, which installs Python dependencies, downloads and processes a copy of the NCBI TaxDump, and downloads and compiles a copy of samtools (Li et al., 2009). Instructions for installation and execution of BlobTools can be found at https://github.com/DRL/blobtools.

Two common BlobTools workflows for taxonomic interrogation of paired-end (PE) read datasets are depicted in the flowchart in Figure 1. Workflow A is targeted at de novo genome assembly projects where there is no preexisting reference genome. Workflow B should be followed where a reference genome is available.

Workflow A (Figure 1A) proceeds through construction a BlobDB data structure based on input files (step A1), visualisation of assembly and generation of tabular output (A2), partitioning of sequence IDs based on user-defined parameters informed by the visualisations (A3) and partitioning of PE reads based on sequence IDs (A4). It should be noted that while the BlobTools module create (step A1) supports multiple mapping formats, it is recommended that these are processed in advance using map2cov. Generation of tabular ‘hits’ files is simplified through the module taxify, which allows annotation of similarity search results based on TaxID mapping files or based on custom user input in tabular format.

BlobTools can process both PE and single-end read files. The module bamfilter in step A4 is only of relevance if PE read data is used, since single end read data can easily be partitioned using GNU grep or other tools. The module bamfilter can be controlled with a list of sequence IDs to include or to exclude. Use of an exclusion list causes all sequence IDs, except those specified, to be included. In both cases it will output up to four interleaved FASTQ files depending on the actual mapping behaviour of the read pairs and whether the parameter --include_unmapped is provided. Possible mapping behaviours of read pairs are: both reads mapping to included sequences (included-included: InIn), one read mapping to an included sequence and the other being unmapped (InUn), and one read mapping to an included sequence and the other mapping to an excluded sequence (ExIn). If the --include_unmapped parameter is specified, the module also writes read pairs where neither read maps to the assembly (UnUn). The latter case can occur if the assembler used for generating the sequences did not make use of all reads in the dataset. The resulting partitioned PE read files can then be assembled separately and the workflow is repeated. Decisions concerning which PE read files to use is left at the discretion of the user. However, as general rule, if target taxa have been sequenced at low coverages it might be preferable to be inclusive (using InIn, InEx, InUn and UnUn FASTQ files for assembly) and risking including non-target reads, than being exclusive (using only InIn and InUn for assembly) and risking losing significant proportions of reads from target genomes.

Workflow B (Figure 1B) should be applied when a reference genome is available. Reads are mapped against the reference genome (B1) and the resulting BAM file is processed with the module bamfilter (B2) using the parameter --include_unmapped and without providing a list of sequences. This will result in three FASTQ files: InIn, InUn and UnUn. Since taxonomic origin of the InIn and InUn reads has been established through the mapping step, only the UnUn reads are assembled de novo (B3) and processed via workflow A. This decreases computational requirements substantially. If workflow A yields a PE read partition of the target organism, which will consist of parts of the organism’s genome not present in the reference, these reads are can be used together with the InIn and InUn reads from step 2 to generate a new assembly (B4), which should be screened again via Workflow A. This iterative procedure can easily be applied to projects studying highly variable species where segmental presence-absence is common and a reference genome is expanded (to form a pangenome) as new samples are sequenced, or holobiomes, where reference genomes of multiple taxa are expanded as new samples are added.

Data

To illustrate workflow A (Figure 1A), we simulated read libraries for the nematode Caenorhabditis elegans contaminated with other organisms (see Table 2). Library A contains C. elegans reads contaminated with reads from Escherichia coli, Homo sapiens chromosome 19 and H. sapiens mitochondrial (mtDNA) genome, mimicking a dataset where the target genome is contaminated with DNA from food (E. coli) and operator (H. sapiens). Library B is composed of C. elegans reads contaminated with Pseudomonas aeruginosa, mimicking a project where the metazoan target species is heavily colonised by a prokaryotic organism.

Taxonomic interrogation and partitioning of read pairs using BlobTools

We assembled both read datasets together and mapped each library individually against the assembly. We supplied the assembly to BlobTools, in addition to coverage information extracted from both BAM files and the results of sequence similarity searches.

To simulate cases where sequences of genomes in the assembly are not part of public sequence databases, we removed all sequences annotated under the taxonomic terms ‘Caenorhabditis elegans’, ‘Hominids’, ‘Escherichia’, ‘Pseudomonas’, and ‘Other sequences’ before conducting sequence similarity searches. The search results provided to BlobTools were BLASTn megablast search against NCBI nt (-outfmt ’6 qseqid staxids bitscore std’ -max-target-seqs 1 -max_hsp 1 -evalue 1e-25) and Diamond blastx searches against UniProt Reference Proteomes (--outfmt 6 --sensitive --max-target-seqs 1 --evalue 1e-25), supplied in this order and using taxrule ‘bestsumorder’.

A BlobPlot (Figure 2A), ReadCovPlots (Figure 2B and C) and a CovPlot (Figure 3) were generated at the taxonomic rank of ’order’. A tabular view of the BlobDB was generated using the module view under the taxrule ’bestsumorder’ and for the taxonomic ranks of ’superkingdom’, ’phylum’, and ’order’. We partitioned sequences based on differential coverage and taxonomy annotation (Figure 3) using the tabular view and the UNIX tools GNU grep, GNU cut, and GNU awk. Subsequently, read pairs were partitioned based on mapping behaviour to these sequence partitions using the module bamfilter and read pairs where both reads mapped to included sequences (i. e. the InIn set) were assembled by taxonomic group.

We then generated BlobPlots for the four assemblies (named ‘rhabditida-BT’, ‘primates-BT’, ‘pseudomonadales-BT’ and ‘enterobacterales-BT’) (Figure 4). Coverage information was based on mapping of both simulated sequencing libraries against all four assemblies and sequences were coloured based on the genome-of-origin of the simulated reads mapping to them.

Figure 4. BlobPlots of assemblies by taxon after read partitioning using BlobTools.

Coverage was obtained by mapping original reads to assemblies. Sequences are taxonomically annotated with ’true’ taxonomy based on origin of simulated reads mapping to them. Sequences labelled as ’no-hit’ did not receive any reads mapped to them. (A) Assembly of partition of Rhabditida reads (’rhabditida-BT’). One P. aeruginosa sequence (span 4,886 nt) remains. (B) Assembly of partition of Primates reads (’primates-BT’). Five E. coli sequences (total span 3,838 nt) remain. (C) Assembly of partition of Pseudomonadales reads (’pseudomonadales-BT’). (D) Assembly of partition of Enterobacterales reads (’enterobacterales-BT’). One sequence of P. aeruginosa (span 254 nt) remains.

Evaluation of results

Cleaned assemblies were evaluated based on the count of simulated reads, by genome-of-origin, mapping to them (Table 3), and based on standard assembly metrics (Table 4).

To account for assembly and mapping biases, the original simulated read sets were also assembled separately by taxon, yielding the assemblies CELEG-SIM (reads simulated from the C. elegans genome), HSAPI-SIM (reads simulated from H. sapiens chromosome 19 and mtDNA), PAERU-SIM (reads simulated from P. aeruginosa genome), and ECOLI-SIM (reads simulated from E. coli genome).

We evaluated the effect of parameters of similarity searches against public databases on taxonomic annotation using BlobTools (see Supplementary File 2). Since exhaustive searches against large databases require time and computing power we focussed on parameters that limit resource usage and control the number of returned results. In both BLASTn and Diamond blastx, the options -max-target-seq and -max-hsps are implemented. The former is an early filter applied during primary search and excludes initial hits from later examination. The latter controls the number of high-scoring pairs (HSPs) reported between a query and a subject in the search. The BLAST specific parameter -culling-limit controls the number of hits that can be allocated to a given region on the query. For this dataset, the best trade-off between false positive and false negative taxonomic annotations was achieved by combining BLAST search (-max-target-seqs 10 -evalue 1e-25) against NCBI nt with Diamond blastx searches (--evalue 1e-25 --max-target-seqs 1) against UniProt Reference Proteomes, in this order, using BlobTools taxrule ’bestsumorder’. However, a much faster search with acceptable outcome was achieved by changing the BLASTn parameters to -max-target-seqs 1 -max_hsps 1.