Merging multiple isolate data into a new model

[addyblanch]

We are hoping to create a genus specific model due to some difficulties we have had with base calling. We have 4 sets of data off a minION run with corresponding reference genomes. Any info on how to merge them into one model would be great.

Thanks
Adam

[iiSeymour]

Assuming you have already partitioned your data (say into directories {p1, p2, p3, p4}) step one would be to run 4 basecalling jobs with --save-ctc using each reference.

/data/p1 $ bonito basecaller dna_r9.4.1 --save-ctc --reference ref_1.mmi . > p1.sam
/data/p2 $ bonito basecaller dna_r9.4.1 --save-ctc --reference ref_2.mmi . > p2.sam
/data/p2 $ bonito basecaller dna_r9.4.1 --save-ctc --reference ref_3.mmi . > p3.sam
/data/p4 $ bonito basecaller dna_r9.4.1 --save-ctc --reference ref_4.mmi . > p4.sam

Each job will create three .npy files required for training {chunks.npy, reference_lengths.npy, references.npy}. Note - when basecalling with --save-ctc the training data is accumulated in memory before being written out as the .npy files so you'll need a decent amount of RAM (or limit the number of reads with --max-reads). The final step before training a model (either from scratch or fine tuning the pretrained model provided) is to stack the 4 sets of .npy files together (I can provide code for this).

Does that make sense? Happy to clarify further.

[addyblanch]

Yes I've had a go with one set of data through training and re-basecalling so makes sense. Have four isolates to combine hopefully this will improve the assemblies as just one didn't. Is that normal?

Would be great if you could share the code to merge training sets please.
A

[addyblanch]

Any update to this Chris?

[iiSeymour]

Hey @addyblanch

Sorry, here you go - just set OUTDIR and DIRS appropriately.

import os
import numpy as np

OUTDIR = '/tmp'
DIRS = [
    '/data/p1',
    '/data/p2',
    '/data/p3',
    '/data/p4',
]

LENS = []
REFS = []
CHUNKS = []

def rectify(ragged_array, max_len=None):
    lengths = np.array([len(x) for x in ragged_array], dtype=np.uint16)
    padded = np.zeros((len(ragged_array), max_len or np.max(lengths)), dtype=ragged_array[0].dtype)
    for x, y in zip(ragged_array, padded):
        y[:len(x)] = x
    return padded, lengths


for directory in DIRS:
    LENS.append(np.load(os.path.join(directory, 'reference_lengths.npy')))
    REFS.append(np.load(os.path.join(directory, 'references.npy')))
    CHUNKS.append(np.load(os.path.join(directory, 'chunks.npy')))

lens = np.concatenate(LENS)
refs, _ = rectify([x for ref in REFS for x in ref])
chunks = np.concatenate(CHUNKS, axis=0)

np.save(os.path.join(OUTDIR, 'reference_lengths.npy'), lens)
np.save(os.path.join(OUTDIR, 'references.npy'), refs)
np.save(os.path.join(OUTDIR, 'chunks.npy'), chunks)

[addyblanch]

Thats great, really appreciate it. Thank you.

[mbhall88]

@iiSeymour that code snippet was very useful. I wonder if you would consider making it a subcommand of bonito? Something like bonito combine?