Network architecture.

Deepbinner is implemented using the TensorFlow [17] and Keras [18] code libraries. Its neural network architecture was based on elements developed in the field of image classification: groups of convolutional layers followed by max pooling layers [19]; parallel ‘inception’ modules and low dimension bottlenecks [20]; noise, dropout and batch normalisation layers [21]; and global average pooling [22].

Using these elements, we trialled hundreds of randomised network architectures to search for an effective design. Networks were assessed on their loss (categorical cross-entropy) and classification accuracy on a validation set. To discourage overfitting, we preferred models with fewer parameters and a small ratio of validation set loss to training set loss. The best performing architecture was subsequently refined to produce the final Deepbinner network shown in Fig 1. One notable way that Deepbinner’s architecture differs from image classification networks is the number of filters. Deepbinner uses a constant filter size (48, except for where the parallel module increases the filter count) whereas image classification networks commonly use a smaller number of filters in early layers and a larger number in later layers.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Neural network architecture.

Layers in the network are drawn as coloured blocks and data as groups of vertical lines. Data dimensions are shown for each step of the process as data length × filter count. Gaussian noise and dropout layers are only active during network training, not during classification.

https://doi.org/10.1371/journal.pcbi.1006583.g001

Deepbinner’s input array size of 1024 was chosen for two reasons. First, it is a power of two, allowing the data to halve in length as it progresses through the network’s dimension-reducing layers. Second, an ONT native barcode is 40 bp in length (24 bp for the barcode itself plus 8 bp of flanking sequence on each side) which has a typical raw signal length of 250 to 600 values, making 1024 the smallest power of two which can reliably capture the entire signal.

Network training.

Our training data came from eight R9.4 flow cells (designed for 1D sequencing) and six R9.5 flow cells (designed for 1D² sequencing but also compatible with 1D sequencing), all used with the EXP-NBD103 barcoding kit. We basecalled the reads using Albacore (v2.3.1) and searched for adapters and barcodes in the read sequences using the Edlib library [23]. Both read-start signal (raw signal from the beginning of the read) and read-end signal (raw signal from the end of the read) were extracted from each fast5 file and any open-pore signals (high current values corresponding to the absence of a DNA molecule in the nanopore [5]) were trimmed off. We then conducted semi-global dynamic time warping [24] between read signals and the expected signals for adapters and barcodes to find those elements in the read signals. Instances with a clear barcode signal were used to generate training samples of length 1024 (to match Deepbinner’s input array size).

In addition to training Deepbinner on barcoded signals, we also included a variety of barcode-less signals in the training set that were assigned a corresponding no-barcode class. These included signals from real sequences that lacked barcodes (taken from non-barcode parts of the reads), adapter signals from non-barcoded sequencing runs and multiple types of simulated signals: flat signal, Gaussian noise and Perlin noise (a coherent noise function which generates a smoother signal than Gaussian noise [25]) (S1 Fig). Their presence in the training set ensures that Deepbinner can actively assign reads to a no-barcode class, not just fail to find a strong match. These no-barcode samples were included at a rate to make up approximately of the training samples, resulting in a total of 3 308 789 read-start samples and 1 029 684 read-end samples, each an array of 1024 electrical current values with a corresponding barcode label. There were fewer read-end training samples due to the fact that barcodes more reliably occur at the starts of reads than the ends.

Data augmentation is a method of artificially expanding a training set by duplicating samples with transformations [26]. Deepbinner applies data augmentation during network training by distorting signals in the temporal dimension (S2 Fig). This is carried out by duplicating current values at of the signal’s positions (randomly chosen) and deleting values at of the positions. The result is a signal of equal length but elongated in some places and shortened in others. We did not implement distortion of the signal amplitude because a similar role is carried out by the network’s Gaussian noise layer (Fig 1). For our dataset, we found that an augmentation factor of two (one augmented signal for each unmodified signal) resulted in the best classification accuracy. However, this factor is adjustable, and datasets with fewer training samples may benefit from more data augmentation.

We refined our training datasets by conducting a four-way split, classifying each quarter of the data using a model trained on the other three quarters. Samples which were misclassified (assigned a label which disagreed with the training label) were discarded. This served to remove ambiguous training samples (approximately 0.25%), leaving final training sets of 3 300 075 and 1 027 469 samples for the starts and ends of reads, respectively (totalling 16.5 GB in size).

We produced the final Deepbinner trained models by training the network for 1000 epochs (100000 training samples per epoch) using a data augmentation factor of two and a random 95:5 training:validation split (S3 Fig). Data augmentation was only performed on the training data which, combined with the network’s Gaussian noise and dropout layers (only active during training), explains why validation loss and accuracy were superior to training loss and accuracy. Each model was trained on a single NVIDIA P100 GPU and took approximately 40 hours to complete.

Read classification.

When classifying raw ONT signal, Deepbinner will generate a probability for each possible barcode, along with a no-barcode probability (Fig 1). Deepbinner will assign a barcode to the signal if the highest probability is sufficiently larger (a difference of >0.5) than the second-highest probability. If the best and second-best matches are too close or if the best-match is ‘no barcode’, then Deepbinner assigns the ‘none’ label to the signal.

Deepbinner can classify each read using the read-start signal, the read-end signal or both. Using both start and end is appropriate for library preparations which add barcodes to both sides of a read, such as EXP-NBD103. If both are used, Deepbinner will independently perform classification using the read start and read end. Reads will be binned if a sufficient match was found on either end, but if the two ends match different barcodes, the read will be considered a chimera and put in the ‘none’ bin. These reads can be identified by running Deepbinner in verbose mode which provides detailed information on barcode calls. Deepbinner can optionally require a positive match on both the start and end to bin a read, enabling stringent demultiplexing. Other library preparations, such as the SQK-RBK004 rapid barcoding kit, only add a barcode to the start of the read, in which case Deepbinner’s start-only classification is appropriate.

An ONT barcode signal often appears in the first 1024 values of a read’s raw signal, but this may not be the case for reads which contain a longer-than-normal amount of open-pore signal. Deepbinner therefore examines multiple 1024-length signal windows overlapping by 512 samples. By default, it examines 11 such windows (covering 6144 samples in total), but this can be configured. Deepbinner merges the results across these windows by retaining the maximum probability for each barcode, then renormalising the probabilities.

Study design.

In order to assess the accuracy with which Deepbinner can assign reads to bins based on barcode signals, we aimed to sequence a library of barcoded DNA molecules for which the input source of each molecule could be verified independently of the barcode. To do this, we chose 12 bacterial isolates of different species that had been sequenced via Illumina HiSeq, which produces high accuracy short reads. To identify a region from each that was entirely unique to that genome, we performed a co-assembly of the pooled Illumina reads for all samples with SPAdes (v3.11.1) [27] using a small k-mer (k = 23). The longest contigs from this assembly were matched to their source genome by comparison to individual (non-pooled) genome assemblies, and we choose a single contig per source genome. This produced 12 sequences, one from each genome, each composed entirely of unique 23-mers, making them easily distinguishable from each other.

Amplicon library preparation and MinION sequencing.

The 12 bacterial isolates were grown overnight at 37°C on LB agar plates. Single colonies were then grown overnight at 37°C in Luria broth. Bacterial cell pellets from 1.5 ml of broth culture were generated by centrifugation at 15 000×g for 5 minutes. DNA was extracted from these pellets using Agencourt GenFind v2 (Beckman Coulter) with minor modifications as follows. Cell pellets were resuspended in 400 μl lysis buffer containing 9 μl Proteinase K (96 mg/ml Beckman Coulter) and 1 μl RNase A (100 mg/ml Sigma Aldrich R6513) by gentle tip mixing. Samples were lysed at 37°C for 30 minutes. gDNA was extracted from the lysed samples by completing the remaining steps of the GenFind v2 for 200 μl of blood/serum from the binding step onwards.

We used the Primer3web tool (v4.1.0) to choose 25 bp long-range PCR primers from each of the 12 unique sequences to define amplicons which ranged from 9 to 11 kbp (S4 File). PCR was performed using LongAMP Taq 2X Master Mix (New England Biolabs) with 150–450 ng gDNA as input template, primers at 1 μM, an annealing temperature of 56°C and 35 cycles of amplification. Following the PCR, size selection purification was performed with Agencourt AMPure beads (Beckman Coulter) at 0.6× ratio. The size and specificity of the PCR product was confirmed by capillary electrophoresis (Fragment Analyser AATI). A sequencing library was prepared from the purified amplicons using the Nanopore 1D ligation sequencing kit (SQK-LSK108) with the native barcoding expansion kit (EXP-NBD103) as per the manufacturer’s instructions. The run was performed on a MinION MK1b device using an R9.4 flow cell (FLO-MIN106), MinKNOW v18.03.1 and the NC_48Hr_Sequencing_Run_FLO-MIN106_SQK-LSK108 protocol.

Demultiplexing.

The reads were basecalled with Albacore (v2.3.1) using the following options: barcoding (to enable demultiplexing) and disable_filtering (to include low-quality reads). The resulting FASTQ files were pooled and shuffled, and then given to Porechop (v0.2.3) to be independently demultiplexed. The pre-basecalled fast5 files were demultiplexed with Deepbinner (v0.2.0).

In addition to running the tools with default parameters, we also tested Porechop and Deepbinner with parameters to increase or decrease their demultiplexing stringency. Porechop’s lenient settings reduce the barcode and difference thresholds from their default values (75% and 5%, respectively) to 60% and 1%. This makes Porechop willing to consider very low-quality alignments and assign a barcode even when there is a close second-best match. Porechop’s stringent settings increase the barcode threshold to 85% (so it will only consider high-quality alignments) and require a barcode match on both the start and end of reads. The only difference between Deepbinner’s default and stringent settings is that the latter requires a barcode match on both the start and end of reads.

We assigned ground truth classifications to the ONT reads by aligning their basecalled sequences to the amplicon reference sequences with minimap2 (v2.12) [28] and binning with the assign_reads_to_reference.py script (included with Deepbinner). Reads which failed to meet an alignment threshold (100 bp or 10% of the read length, whichever is smaller) were classified as ‘unknown’. Reads which exceeded an alignment threshold to a secondary amplicon (50 bp or 5% of the read length, whichever is larger) were classified as ‘chimera’. This method is only able to detect cross-bin chimeras, i.e. reads with separate components from two different amplicons. It cannot detect within-bin chimeras, e.g. two copies of the same amplicon concatenated in a single read.

Additional test sets.

We also performed a whole genome sequencing (WGS) run using 12 different bacterial species, as it represents a more realistic sequencing scenario for many users. We followed the same preparation and analysis described above for the amplicon sequencing but with genomic DNA and excluding the PCR-specific steps. The reads were aligned with minimap2 (v2.12) [28] to all twelve complete genome sequences (produced via hybrid Illumina-Nanopore assembly using Unicycler v0.4.6 [4]) and assigned ground truth classifications with the same script and thresholds used for binning the amplicon reads. Due to shared sequences between the genomes, confidently assigning ground truth was more challenging than for the amplicon set, resulting in more reads being classified as ‘unknown’.

Based on examination of basecalled sequences, we estimate approximately 0–2% of reads in a barcoded read set may be genuinely lacking a barcode, i.e. one was not ligated during preparation. These reads can be difficult to distinguish from reads which do have a barcode but where no classification was possible, so to examine the demultiplexing tools’ behaviour on non-barcoded reads, we additionally assessed them using a non-barcoded ONT read set (SQK-LSK108, R9.4 flow cell). This is referred to as the ‘negative control’ set and the ground truth classification for all its reads was ‘none’.

To ensure a fair assessment, no reads from any of the three test sets (amplicon, WGS and negative control) were used in the training of Deepbinner’s models.