<commit_before># RUscripts

To get started either download this repository using the ZIP options or (preferably)

git clone https://github.com/mattloose/RUscripts



Important Note:

Running read until will influence the behaviour of your flow cell and change the output of your sequencing experiment. You are strongly advised to run simulations of read until prior to running on a live flow cell. The code as presented here is a demonstration of read until and one implementation. Users run this code entirely at their own risk.

#Getting help with read until.

We are happy to help where we can. Please feel free to contact via twitter ({at}mattloose) or email <EMAIL> . We will try to update this read me file whenever significant changes are made to the API or minKNOW versions that might affect read until. For issues with the API itself or the ws_event_sampler program, users will most likely be redirected to ONT.

#Getting started

Scripts for implementing read until and other examples.

To use read until at this time you have to use the command line. To access the command line on windows, click on the start menu and type 'cmd' in the search box. You should see a program called 'cmd'.

Once at the command line, navigate to the folder containing this Readme file.

All scripts are python 2.7.

We assume that the default version of python in your path is the ONT Anaconda installation. To test this, type

    C:\path\to\RUscripts>python

And you should see:

    Python 2.7.5 |Anaconda 1.8.0 (64-bit)| (default, Jul  1 2013, 12:37:52) [MSC v.1500 64 bit (AMD64)] on win32
    Type "help", "copyright", "credits" or "license" for more information.
    >>>

Note that this shows Anaconda 1.8.0

To exit this window type

    >>> exit()

Please note these scripts depend upon APIs which are not provided here. The Read Until API is available from Oxford Nanopore on request to MAP participants.

# minKNOW version compatability statement

These scripts and files are compatible with versions of minKNOW predating the v 0.51.1.62 (released February 22nd 2016). Read Until in this latest version of minKNOW is currently unstable. ONT are fixing this at this time. We provide scripts that compensate for these errors in a simulation, but these will not work for a live read until run at this time.

Currently supported versions of minKNOW are v 0.48.1.3 to 0.51.1.51.

We highlight where the version of minKNOW used might impact results in the documentation below.

*** You must have minKNOW installed on the windows machine you are running this code on. The version of Anaconda it comes with contains some of the required packages. These scripts will also run on Linux or OSX but at this time we provide no explicit documentation to support this.

# Python modules required for running read until using source code.

h5py

configargparse

biopython

mlpy # Must be version 3.5.0 or above - see below

watchdog

sklearn #Available as default with Anaconda


# To install these modules on the ONT Anaconda type:

    C:\path\to\RUscripts>easy_install thrift ws4py h5py configargparse biopython watchdog

NOTE: You may need to add C:\Anaconda\Scripts to your path to access easy_install and pip.

# To install mlpy

Installing mlpy with easy_install results in a version less than that required. To overcome this on windows systems visit http://www.lfd.uci.edu/~gohlke/pythonlibs/#mlpy to obtain an mlpy wheel compiled for the appropriate version of windows and python.

On our system this file is called:

mlpy-3.5.0-cp27-none-win_amd64.whl


To install this you need to first install pip.

    C:\path\to\RUscripts>easy_install install pip

Then to install your downloaded wheel:

    C:\path\to\RUscripts>pip install \path\to\mlpy-3.5.0-cp27-non-win_amd64.whl

Now you should have all the prerequisites required to run these scripts on windows using the native python.


# General comments

All scripts here adhere to the following standard arguments:

-h Will print a help statement and exit.
-ver Will print version information and exit.
-v Will print additional debug information and more verbose output to screen while the code is executing.

Additional options are described below for each script.

# Files and Folders in this repository

ampbalance.py - an illustrative script for selecting a fixed number of reads from individual amplicons using squiggle matching alone.

ampbalancetest (folder - contains 110 example fast5 files)

ampliconSPLIT.py - an illustrative script that will split reads into folders based on squiggle matching to individual amplicons.

exampleread\llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5

getmodels.py - A script for extracting model files to generate a squiggle reference for a given sequence.

J02459.fasta - A reference fasta file for bacteriophage lambda. This is the reference used for all our work.

lambda_amplicons.txt - A file that defines the amplicons used in this study. The file is in the format seq_id:start-stop - so here the seq_id is J02456 (as in the file J02459.fasta) and the start and stop coordinates match our amplicons eg. J02456:52-1980

LICENSE - MIT License

README.md - this readme file.

ReadUntil\aReadUntil.py - A script for amplicon balancing via the read until API.

ReadUntil\gReadUntil.py - A script for selecting a specific region of the genome via the read until API

ReadUntil\lambda_amplicons.fasta - This is a reference file containing copies of individual amplicons for bacteriophage lambda. It is used in testing the behaviour of aReadUntil.py - described below.

ReadUntil\ruutils.py - This file contains a number of useful functions for various scripts. It will never be run itself but is called by many of the other scripts.

test_gReadUntil.py - A script for testing the functionality of gReadUntil.py in the absence of the API itself.

RUtestset (folder - contains 55 example fast5 files)


# Running Offline Read Until Scripts

# getmodels.py

In order to squiggle match between an incoming read and a reference sequence it is necessary to generate a reference in squiggle space. To do this you require a model file to enable the conversion of kmers in base space to a current value. Currently the appropriate model files are embedded within reads as they are returned by the basecaller.

It is important to have a model file for the chemistry and pore type that you are using.

For read until you only require a template model file. If you wish to run ampbalance (which identifies potential 2D reads from a dataset you also require the complement model file).

It is not necessary to generate a model file each time you carry out a new read until experiment. Instead, new files should be obtained whenever a chemistry or pore change occurs. We suggest that it is good practice to extract a model file from the most recent read data available to you when preparing for a read until run as these models have been silently updated in the past.

Ideally the read from which to extract a model should come from the pass folder of a base called run. To extract both template and complement models, this read must be from a 2D base called run. getmodels.py will output whichever models it finds within the file.

To print the getmodels.py help statement at the prompt type:

    C:\path\to\RUscripts>python getmodels.py -h

which will output:

    usage: getmodels.py [-h] [-r READ] [-v] [-ver]

    getmodels: A program to extract model files from Nanopore reads. Ideally you
    should provide a read that has given a 2D read - i.e comes from the pass
    folder on a 2D basecalling run. For Read Until purposes, extraction of just a
    template model is sufficient. Model files will be appended with a number to
    identify the kmer length of the model.

    optional arguments:
    -h, --help            show this help message and exit
    -r READ, --read READ  Provide a read file to extract the current model from.
    -v, --verbose-true    Print detailed messages while processing files to aid
                        in debugging.
    -ver, --version       Print version information and exit.

To run getmodels.py at the prompt ($) type:

    $python getmodels.py -r /path/to/your/read.fast5

We provide a suitable read in the folder exampleread

So:

    C:\path\to\RUscripts>python getmodels.py -r exampleread/llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5

will output:

    Processing exampleread/llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5
    Looking for template model.
    ***: Model Found
    template model file write completed.
    Kmer length is: 6
    File format is:
    Kmer	Mean	Standard Dev
    Looking for complement model.
    ***: Model Found
    complement model file write completed.
    Kmer length is: 6
    File format is:
    Kmer	Mean	Standard Dev

This will have generated two files - one for the template and one for the complement model.

They will be named:

[readtype]_[model_name]_[kmerlength].model

Where [readtype] is either template or complement, [model_name] is the pore/chemistry type (in this case r7.3_e6_70bps_6mer) and [kmerlen] is the locally calculated length of the kmers used in the model.

These model files are used by all other scripts here. Note that the model files will be overwritten by getmodels.py if it is run again on another read with the same pore/chemistry type and kmer length. This is not a problem.

*** We recommend that you update the model regularly.

# ampbalance.py

Ampbalance is a script which will map entire 2D reads to a reference sequence, identify those likely to give a good 2D read by determining if the template and complement reads overlap with one another appropriately and subsequently copy a defined number of those reads up to a specific depth into a folder for further analysis. This script was originally written to facilitate the sequencing and analysis of Ebola virus and to simplify the transfer of read data prior to base calling where network bandwidth was limited (see <NAME>., <NAME>., <NAME>., <NAME>. & <NAME>. Real-time, portable genome sequencing for Ebola Surveillance. Nature, doi:10.1038/nature16996 (2016).)

Unlike read until applications this script uses the entire read.

Here the script is implemented on reads that have been basecalled to facilitate testing and comparisons.

To print the ampbalance.py help statement at the prompt ($) type:

    C:\path\to\RUscripts>python ampbalance.py -h

which will output:

    usage: ampbalance.py [-h] -fasta FASTA -ids IDS -w WATCHDIR -o TARGETPATH -d
                         DEPTH -procs PROCS [-cautious] [-l LENGTH] -t TEMP_MODEL
                         -c COMP_MODEL [-v] [-ver]

    ampbalance: A program designed to balance amplicons from a specific reference
    sequence post sequencing on ONT minIONs but prebasecalling. Developed by Matt
    Loose @mattloose or <EMAIL> for help! Args that start with
    '--' (eg. --reference_fasta_file) can also be set in a config file
    (/path/to/your/script/amp.config or ) by using .ini or
    .yaml-style syntax (eg. reference_fasta_file=value). If an arg is specified in
    more than one place, then command-line values override config file values
    which override defaults.

    optional arguments:
      -h, --help            show this help message and exit
      -fasta FASTA, --reference_fasta_file FASTA
                            The fasta format file for the reference sequence for
                            your organism.
      -ids IDS, --reference_amplicon_positions IDS
                            A file containing a list of amplicon positions defined
                            for the reference sequence. 1 amplicon per line in the
                            format fasta_sequence_name:start-stop e.g
                            J02459:27-1938
      -w WATCHDIR, --watch-dir WATCHDIR
                            The path to the folder containing the downloads
                            directory with fast5 reads to analyse - e.g.
                            C:\data\minion\downloads (for windows).
      -o TARGETPATH, --output-dir TARGETPATH
                            The path to the destination folder for the
                            preprocessed reads
      -d DEPTH, --depth DEPTH
                            The desired coverage depth for each amplicon. Note
                            this is unlikely to be achieved for each amplicon and
                            should probably be an overestimate of the minimum
                            coverage required.
      -procs PROCS, --proc_num PROCS
                            The number of processors to run this on.
      -cautious, --cautious
                            DTW of long reads on low memory systems can cause
                            unexpected crashes. This option will prevent automatic
                            skipping on any reads over 10,000 events. You can
                            optionally increase this length with the -l parameter.
                            USE WITH CAUTION AS THIS MAY CAUSE A SYSTEM TO CRASH.
      -l LENGTH, --length LENGTH
                            A limit on the length of read that ampbalance will
                            attempt to align using DTW - Long reads can cause
                            problems on low memory systems
      -t TEMP_MODEL, --template_model TEMP_MODEL
                            The appropriate template model file to use
      -c COMP_MODEL, --complement_model COMP_MODEL
                            The appropriate complement model file to use
      -v, --verbose-true    Print detailed messages while processing files.
      -ver, --version       show program's version number and exit

This script is designed to match amplicon sequences of known and approximately uniform length to a reference sequence. The reference sequence should be a single sequence.

We provide an example set of reads using lambda, found in the ampbalancetest folder. A typical usage command would be:

    C:\path\to\RUscripts>python ampbalance.py -fasta J02459.fasta -ids lambda_amplicons.txt -w ampbalancetest -o ampbalanceoutput -d 5 -procs 4 -t template_r7.3_e6_70bps_6mer_6.model -c complement_r7.3_e6_70bps_6mer_6.model -l 3000

This will output the following:

    Reading amplicons
    ******AMP DICTIONARY*******
    <type 'dict'>
    {1: 52, 2: 2065, 3: 4070, 4: 6059, 5: 8012, 6: 10008, 7: 12006, 8: 14011, 9: 16076, 10: 18022, 11: 20053}
    {'DO': 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'HF': 0, 'NH': 0}
    Now we are going to try and open the raw reads and do the same as we have done above...
    {'DO': 0, 1: 0, 2: 0, 3: 1, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'TF': 109, 'HF': 4, 'NH': 0}
    {'DO': 0, 1: 0, 2: 0, 3: 1, 4: 0, 5: 0, 6: 0, 7: 1, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'TF': 108, 'HF': 5, 'NH': 0}
    .
    . (many more lines here)
    .
    {'DO': 0, 1: 10, 2: 10, 3: 10, 4: 10, 5: 9, 6: 10, 7: 10, 8: 10, 9: 10, 10: 10, 11: 10, 'BF': 0, 'TF': 1, 'HF': 110, 'NH': 0}
    {'DO': 0, 1: 10, 2: 10, 3: 10, 4: 10, 5: 10, 6: 10, 7: 10, 8: 10, 9: 10, 10: 10, 11: 10, 'BF': 0, 'TF': 0, 'HF': 110, 'NH': 0}
    Amplicon Read Counts
    Amplicon Number: 1 Reads: 10
    Amplicon Number: 2 Reads: 10
    Amplicon Number: 3 Reads: 10
    Amplicon Number: 4 Reads: 10
    Amplicon Number: 5 Reads: 10
    Amplicon Number: 6 Reads: 10
    Amplicon Number: 7 Reads: 10
    Amplicon Number: 8 Reads: 10
    Amplicon Number: 9 Reads: 10
    Amplicon Number: 10 Reads: 10
    Amplicon Number: 11 Reads: 10
    Copying Amplicon Data
    Amplicon Number 1
    Amplicon Number 2
    Amplicon Number 3
    Amplicon Number 4
    Amplicon Number 5
    Amplicon Number 6
    Amplicon Number 7
    Amplicon Number 8
    Amplicon Number 9
    Amplicon Number 10
    Amplicon Number 11

Key:
DO: Template and Complement Strands don't overlap on amplicons as expected. Read Discarded. Unlikely to generate good 2D.
1...11: Amplicon numbers followed by current coverage counts satisfying highest stringency tests.
BF: Bad Files - these files contain more events than the -length threshold. They are unlikely to be of interest in an amplicon run and are skipped unless the -cautious flag is set.
TF: Total Files - this counts down the number of files left to process.
HF: Hairpin Found - this is the number of reads containing a hairpin - these are tested.
NH: No Hairpin - these reads are ignored as no hairpin has been found.

This program will output files to the specified output directory. Files which are already basecalled will be written to a subfolder called "Downloads". Raw files will just be written to the specified output directory.

Note that the example read set here has been preselected to contain 10 reads mapping to each amplicon and all the reads are 2D - thus the HF value is 110 and the DO value is 0.

It is also important to note that this script uses the hairpin_found flag in the read files. Since minKNOW version 0.51.1.62 (released February 22nd 2016) the writing of this flag to read files has been problematic. Thus reads generated with this version of minKNOW may appear to have far fewer 2D reads than they really do.


# ampliconSPLIT.py

ampliconSPLIT.py simulates read until on either raw or basecalled reads. It ignores the first 50 events of a read and matches the subsequent 250 events against a reference squiggle. Read until only ever processes template data and so only a template model is required.

This script will process reads into subfolders in the targetpath corresponding to each amplicon described in the -ids file - e.g /targetpath/{amplicon_number}. If the read is basecalled, it will be placed in a downloads subfolder - e.g /targetpath/{amplicon_number}/downloads .

To print the ampliconSPLIT.py help statement at the prompt ($) type:

    C:\path\to\RUscripts>python ampliconSPLIT.py -h

which will output:

    usage: ampliconSPLIT.py [-h] -fasta FASTA -ids IDS -w WATCHDIR -o TARGETPATH
                        -d DEPTH -procs PROCS -t TEMP_MODEL [-v] [-ver]

    ampliconSPLIT: A program designed to identify and group individual amplicons
    from minION reads prior to base calling. The depth setting limits the number
    of reads copied to each sub folder. Developed by <NAME> @mattloose or
    <EMAIL> for help! Args that start with '--' (eg.
    --reference_fasta_file) can also be set in a config file
    (/Users/mattloose/fixes/RUscripts/UPDATES/amp.config or ) by using .ini or
    .yaml-style syntax (eg. reference_fasta_file=value). If an arg is specified in
    more than one place, then command-line values override config file values
    which override defaults.

    optional arguments:
    -h, --help            show this help message and exit
    -fasta FASTA, --reference_fasta_file FASTA
                        The fasta format file for the reference sequence for
                        your organism.
    -ids IDS, --reference_amplicon_positions IDS
                        A file containing a list of amplicon positions defined
                        for the reference sequence. 1 amplicon per line in the
                        format fasta_sequence_name:start-stop e.g
                        EM_079517:27-1938
    -w WATCHDIR, --watch-dir WATCHDIR
                        The path to the folder containing the downloads
                        directory with fast5 reads to analyse - e.g.
                        C:\data\minion\downloads (for windows).
    -o TARGETPATH, --output-dir TARGETPATH
                        The path to the destination folder for the
                        preprocessed reads
    -d DEPTH, --depth DEPTH
                        The desired coverage depth for each amplicon. Note
                        this is unlikely to be achieved for each amplicon and
                        should probably be an overestimate of the minimum
                        coverage required.
    -procs PROCS, --proc_num PROCS
                        The number of processors to run this on.
    -t TEMP_MODEL, --template_model TEMP_MODEL
                        The appropriate template model file to use. This file
                        can be generated uing the getmodels.py script.
    -v, --verbose-true    Print detailed messages while processing files.
    -ver, --version       show program's version number and exit

An example run command using test data:

    C:\path\to\RUscripts>python ampliconSPLIT.py -fasta J02459.fasta -ids lambda_amplicons.txt -w RUtestset/ -o test -d 10 -procs 8 -t template_r7.3_e6_70bps_6mer_6.model

Which will output:

    processing the reference fasta.
    ID J02459
    length 48502
    FORWARD STRAND
    REVERSE STRAND
    Groking amplicons
    ******AMP DICTIONARY*******
    <type 'dict'>
    {1: 52, 2: 2065, 3: 4070, 4: 6059, 5: 8012, 6: 10008, 7: 12006, 8: 14011, 9: 16076, 10: 18022, 11: 20053}
    {'DO': 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'HF': 0, 'NH': 0}
    We want to build a custom reference that is smaller than the original reference.
    First get a list of all the positions we will need to search
    Generating a custom fasta
    J02459
    processing the custom fasta
    Attempting to match reads and split into folders based on 250 events, excluding the first 50.
    Amplicon Read Counts
        Amplicon Number: 1 Reads: 4
        Amplicon Number: 2 Reads: 5
        Amplicon Number: 3 Reads: 4
        Amplicon Number: 4 Reads: 5
        Amplicon Number: 5 Reads: 6
        Amplicon Number: 6 Reads: 6
        Amplicon Number: 7 Reads: 5
        Amplicon Number: 8 Reads: 5
        Amplicon Number: 9 Reads: 5
        Amplicon Number: 10 Reads: 5
        Amplicon Number: 11 Reads: 5
    Copying Amplicon Data
    Amplicon Number 1
    Amplicon Number 2
    Amplicon Number 3
    Amplicon Number 4
    Amplicon Number 5
    Amplicon Number 6
    Amplicon Number 7
    Amplicon Number 8
    Amplicon Number 9
    Amplicon Number 10
    Amplicon Number 11

The depth parameter (-d) sets the number of reads that will be copied. Setting this to a value greater than the number of reads analysed will sort all reads in the dataset.


# test_gReadUntil.py

*** Note this script is located with the ReadUntil folder of this repository

This script runs entirely independently of the read until API and allows for simulation of selective sequencing of a genome. You can provide a pool of reads and the script will use the first 250 events to map a read and copy those mapping within the desired area to a specified folder.

Help is available by typing:

    C:\path\to\RUscripts\ReadUntil\python test_gReadUntil.py -h

which will output:

    usage: test_gReadUntil.py [-h] -fasta FASTA -targets [TARGETS [TARGETS ...]]
                          -procs PROCS -m TEMP_MODEL [-log LOGFILE] -w
                          WATCHDIR [-o OUTPUT_FOLDER] [-v] [-ver]

    real_read_until: A program providing read until with the Oxford Nanopore
    minION device. This program will ultimately be driven by minoTour to enable
    selective remote sequencing. This program is heavily based on original code
    generously provided by Oxford Nanopore Technologies.

    optional arguments:
    -h, --help            show this help message and exit
    -fasta FASTA, --reference_fasta_file FASTA
                        The fasta format file describing the reference
                        sequence for your organism.
    -targets [TARGETS [TARGETS ...]]
                        Positional IDs to enrich for in the form seqid:start-
                        stop . Can be space seperated eg: J02459:10000-15000
                        J02459:35000-40000
    -procs PROCS, --proc_num PROCS
                        The number of processors to run this on.
    -m TEMP_MODEL, --model TEMP_MODEL
                        The appropriate template model file to use
    -log LOGFILE, --log-file LOGFILE
                        The name of the log file that data will be written to
                        regarding the decision made by this program to process
                        read until.
    -w WATCHDIR, --watch-dir WATCHDIR
                        The path to the folder containing the downloads
                        directory with fast5 reads to analyse - e.g.
                        C:\data\minion\downloads (for windows).
    -o OUTPUT_FOLDER, --output OUTPUT_FOLDER
                        Path to a folder to symbolically place reads
                        representing match and not match.
    -v, --verbose-true    Print detailed messages while processing files.
    -ver, --version       show program's version number and exit

A typical command line to select reads mapping from 10-15kb in the lambda genome would be:

    C:\path\to\RUscripts\ReadUntil\python test_gReadUntil.py -fasta ../J02459.fasta -targets J02459:10000-15000 -procs 4  -m ../template_r7.3_e6_70bps_6mer_6.model -w ../RUtestset/ -o test

This would give the following output:

    ***********************************************************************************************
    **** This code will open a collection of reads and simulate read until on them. It will    ****
    **** copy reads into a secondary folder for subsequent processing by another analysis      ****
    **** package.                                                                              ****
    ***********************************************************************************************
    processing the reference fasta.
    J02459:10000-15000
    ID J02459
    length 48502
    1 J02459:10000-15000 10000 15000 J02459
    We want to extract this chunk J02459_1
    ID J02459_1
    length 5000
    FORWARD STRAND
    REVERSE STRAND
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch10_file49_strand.fast5 No Match
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch10_file69_strand.fast5 Sequence Found
    .
    .
    .
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch8_file58_strand.fast5 Sequence Found
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch9_file100_strand.fast5 No Match
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch8_file11_strand.fast5 No Match

Sequence found indicates that the read is derived from the desired region. No Match indicates that the read is from another region.


# Real Read Until

Important Note:

Running read until will influence the behaviour of your flow cell and change the output of your sequencing experiment. You are strongly advised to run simulations of read until prior to running on a live flow cell. The code as presented here is a demonstration of read until and one implementation. Users run this code entirely at their own risk.


We provide the two scripts used to generate the data presented in "Real time selective sequencing using nanopore technology".

These scripts absolutely require the read until API from Oxford Nanopore.

# Read Until API

This version of Read Until requires that the user runs a command window as Administrator.

The read until API changes depending on the version of minKNOW in use. The read until API for minKNOW versions 0.48.1.3 - 0.51.1.51 is available on request from Oxford Nanopore and the code provided here is compatible with Git Commit

    50ee65c1404870a2e04cb6dc211dcdeeddd3fd4f

and minKNOW versions <= 0.51.1.51

Note that this is not the latest commit. The latest version as of 30th March 2016 does not function as anticipated. It will work for simulation runs with ws_event_sample, but minKNOW itself does not function correctly with the API. This is actively being resolved by ONT. Thus, testing and simulation is possible with our code and later versions of minKNOW and ws_event_sampler but these scripts WILL NOT WORK AT THIS TIME with minKNOW > 0.51.1.51

# Installing the API

Copy the following files and folders (and there contents) into the RUscripts\ReadUntil folder:

    read_until.py
	event_sampler.thrift
    event_sampler/* (note - you need to copy the folder and its contents)
    example.py

# Running ws_event_sampler

Read until works by streaming data from minKNOW during a run to an external analysis pipeline and receiving instructions to reject reads as required. The program streaming the data is called ws_event_sampler.exe and this is found in the folder

    C:\grouper\binaries\bin

which is available after installing minKNOW.

ws_event_sampler can run either stand alone in simulation mode, or alongside minKNOW as it is running. For simulations, we suggest you start minKNOW and have it open (but not actually sequencing). Then start ws_event_sampler in simluation mode. In this way you will see the output from our scripts in the minKNOW messages window.

Running ws_event_sampler is done as follows:

    C:\grouper\binaries>bin\ws_event_sampler.exe

*** NOTE: The command itself is bin\ws_event_sampler.exe so the command is executed from the folder binaries, not from within the folder bin. Executing the file from within the bin folder will result in the program crashing with a message about being unable to open a config file.

*** NOTE: ws_event_sampler provides no visible indication that it is running other than blocking the command line. If it fails, it will terminate - sometimes with an error message, but often not. Support for ws_event_sampler is available direct from ONT. We hope we have provided sufficient guidance for its use in the context of our scripts here.

ws_event_sampler can be run in simulation mode to facilitate the development of analysis scripts. When running in simulation mode, ws_event_sampler provides reads driven from a model file stored locally. This same model file must be used to analyse data when in simulation mode. When streaming live data direct from minKnow it is essential to use the correct model for the chemistry and pore type in use.

The read until API provides comprehensive instructions on how to run ws_event_sampler. Below we provide worked examples for our scripts.

To test that ws_event_sampler is functioning correctly, do the following:

1) Open a command window as administrator - to do this, right click on the cmd icon and select "Run as Administrator"

2) In this window navigate to C:\grouper\binaries

3) To see help for ws_event_sampler type:

    C:\grouper\binaries>bin\ws_event_sampler.exe -h

which will output:

    Usage: ws_event_sampler [Options]

    Common options:
      -h [ --help ]                         Produce help message
      -c [ --config ] arg (=./conf/global.conf)
                                            Application config file.

    Specific options:
      -p [ --port ] arg                     Port to listen for requests. If this is
                                            not specified it is taken from the
                                            ws_event_sampler_port field in the
                                            config file.
      -s [ --sim ]                          Start in standalone simulation mode.
      --sim-channels arg                    Number of channels simulation uses.
                                            E.g. use a low number when debugging to
                                            make the pushed data manageable. This
                                            overrides the channel_count setting in
                                            the conf file but only during
                                            simulation.
      --sim-model arg                       Path to model file to be used by
                                            simulation when generating events.
      --sim-fasta arg                       Path to fasta file from which to draw
                                            fragments to be used by simulation.
      --sim-delay arg (=0)                  Seconds delay before starting
                                            simulation. By default the simulation
                                            starts as soon as this server starts
                                            up.
      --sim-fragment-length arg             Random lengths are bounded, approximate
                                            normal distributions, picked by rolling
                                            lots of dice. You need to specify die
                                            sides, die count and a base amount. For
                                            example, the default setting of '5 250
                                            500' will generate reads with mean
                                            length 1000 and bounds [500,1500].
      --sim-log arg                         Filename to use for csv-format log of
                                            the simulation. Will overwrite if
                                            already exists.

To set up a simulation on a specific port type:

    C:\grouper\binaries>bin\ws_event_sampler.exe -s --sim-channels 100 --sim-fasta C:\path\to\RUscripts\J02459.fasta --sim-log log.txt

This will establish the data stream on port 12345, simulating 100 channels with a read distribution as described derived from the lambda sequence we have provided here. ws_event_sampler in simulation mode will write out a log file to log.txt enabling tracking of events.

In the example.py script, edit the line which states:

    host="ws://localhost:<port>"

to

    host="ws://localhost:12345"

Note - you can specify a different port value if you wish.

Then from another command window execute the examply.py script (note this is supplied by the ONT API)

    C:\path\to\RUscripts\ReadUntil>python example.py

This should output something like (depending on the version of the API in use):

    Client connection started. Beginning unblock loop...
    channel_name= 344  read_number= 12  events_in_sample= 200  first event: start= 1104931.0  mean= 71.1864189873
    channel_name= 345  read_number= 12  events_in_sample= 200  first event: start= 1111035.0  mean= 68.1394555562
    channel_name= 346  read_number= 12  events_in_sample= 200  first event: start= 1104312.0  mean= 67.1975140682
    channel_name= 347  read_number= 15  events_in_sample= 200  first event: start= 1203518.0  mean= 72.436242832
    channel_name= 340  read_number= 12  events_in_sample= 200  first event: start= 1104326.0  mean= 60.9232571465
    channel_name= 341  read_number= 12  events_in_sample= 200  first event: start= 1102520.0  mean= 69.496870232
    channel_name= 342  read_number= 12  events_in_sample= 200  first event: start= 1097286.0  mean= 65.2731245584
    channel_name= 343  read_number= 12  events_in_sample= 200  first event: start= 1102224.0  mean= 65.2680067692
    channel_name= 348  read_number= 12  events_in_sample= 200  first event: start= 1119138.0  mean= 77.3645669056
    channel_name= 349  read_number= 12  events_in_sample= 200  first event: start= 1092464.0  mean= 72.4693304753
    channel_name= 298  read_number= 12  events_in_sample= 200  first event: start= 1109404.0  mean= 66.3816659282
    channel_name= 299  read_number= 12  events_in_sample= 200  first event: start= 1103390.0  mean= 66.1452746381

If you see:

    Client connection started. Beginning unblock loop...
    Unhandled exception in thread started by <bound method ReadUntil.data_received_own_thread of <read_until.ReadUntil object at 0x00000000023AA
    FD0>><commit_msg>Changing suggested port to 9200 for compatability<commit_after># RUscripts

To get started either download this repository using the ZIP options or (preferably)

git clone https://github.com/mattloose/RUscripts



Important Note:

Running read until will influence the behaviour of your flow cell and change the output of your sequencing experiment. You are strongly advised to run simulations of read until prior to running on a live flow cell. The code as presented here is a demonstration of read until and one implementation. Users run this code entirely at their own risk.

#Getting help with read until.

We are happy to help where we can. Please feel free to contact via twitter ({at}mattloose) or email <EMAIL> . We will try to update this read me file whenever significant changes are made to the API or minKNOW versions that might affect read until. For issues with the API itself or the ws_event_sampler program, users will most likely be redirected to ONT.

#Getting started

Scripts for implementing read until and other examples.

To use read until at this time you have to use the command line. To access the command line on windows, click on the start menu and type 'cmd' in the search box. You should see a program called 'cmd'.

Once at the command line, navigate to the folder containing this Readme file.

All scripts are python 2.7.

We assume that the default version of python in your path is the ONT Anaconda installation. To test this, type

    C:\path\to\RUscripts>python

And you should see:

    Python 2.7.5 |Anaconda 1.8.0 (64-bit)| (default, Jul  1 2013, 12:37:52) [MSC v.1500 64 bit (AMD64)] on win32
    Type "help", "copyright", "credits" or "license" for more information.
    >>>

Note that this shows Anaconda 1.8.0

To exit this window type

    >>> exit()

Please note these scripts depend upon APIs which are not provided here. The Read Until API is available from Oxford Nanopore on request to MAP participants.

# minKNOW version compatability statement

These scripts and files are compatible with versions of minKNOW predating the v 0.51.1.62 (released February 22nd 2016). Read Until in this latest version of minKNOW is currently unstable. ONT are fixing this at this time. We provide scripts that compensate for these errors in a simulation, but these will not work for a live read until run at this time.

Currently supported versions of minKNOW are v 0.48.1.3 to 0.51.1.51.

We highlight where the version of minKNOW used might impact results in the documentation below.

*** You must have minKNOW installed on the windows machine you are running this code on. The version of Anaconda it comes with contains some of the required packages. These scripts will also run on Linux or OSX but at this time we provide no explicit documentation to support this.

# Python modules required for running read until using source code.

h5py

configargparse

biopython

mlpy # Must be version 3.5.0 or above - see below

watchdog

sklearn #Available as default with Anaconda


# To install these modules on the ONT Anaconda type:

    C:\path\to\RUscripts>easy_install thrift ws4py h5py configargparse biopython watchdog

NOTE: You may need to add C:\Anaconda\Scripts to your path to access easy_install and pip.

# To install mlpy

Installing mlpy with easy_install results in a version less than that required. To overcome this on windows systems visit http://www.lfd.uci.edu/~gohlke/pythonlibs/#mlpy to obtain an mlpy wheel compiled for the appropriate version of windows and python.

On our system this file is called:

mlpy-3.5.0-cp27-none-win_amd64.whl


To install this you need to first install pip.

    C:\path\to\RUscripts>easy_install install pip

Then to install your downloaded wheel:

    C:\path\to\RUscripts>pip install \path\to\mlpy-3.5.0-cp27-non-win_amd64.whl

Now you should have all the prerequisites required to run these scripts on windows using the native python.


# General comments

All scripts here adhere to the following standard arguments:

-h Will print a help statement and exit.
-ver Will print version information and exit.
-v Will print additional debug information and more verbose output to screen while the code is executing.

Additional options are described below for each script.

# Files and Folders in this repository

ampbalance.py - an illustrative script for selecting a fixed number of reads from individual amplicons using squiggle matching alone.

ampbalancetest (folder - contains 110 example fast5 files)

ampliconSPLIT.py - an illustrative script that will split reads into folders based on squiggle matching to individual amplicons.

exampleread\llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5

getmodels.py - A script for extracting model files to generate a squiggle reference for a given sequence.

J02459.fasta - A reference fasta file for bacteriophage lambda. This is the reference used for all our work.

lambda_amplicons.txt - A file that defines the amplicons used in this study. The file is in the format seq_id:start-stop - so here the seq_id is J02456 (as in the file J02459.fasta) and the start and stop coordinates match our amplicons eg. J02456:52-1980

LICENSE - MIT License

README.md - this readme file.

ReadUntil\aReadUntil.py - A script for amplicon balancing via the read until API.

ReadUntil\gReadUntil.py - A script for selecting a specific region of the genome via the read until API

ReadUntil\lambda_amplicons.fasta - This is a reference file containing copies of individual amplicons for bacteriophage lambda. It is used in testing the behaviour of aReadUntil.py - described below.

ReadUntil\ruutils.py - This file contains a number of useful functions for various scripts. It will never be run itself but is called by many of the other scripts.

test_gReadUntil.py - A script for testing the functionality of gReadUntil.py in the absence of the API itself.

RUtestset (folder - contains 55 example fast5 files)


# Running Offline Read Until Scripts

# getmodels.py

In order to squiggle match between an incoming read and a reference sequence it is necessary to generate a reference in squiggle space. To do this you require a model file to enable the conversion of kmers in base space to a current value. Currently the appropriate model files are embedded within reads as they are returned by the basecaller.

It is important to have a model file for the chemistry and pore type that you are using.

For read until you only require a template model file. If you wish to run ampbalance (which identifies potential 2D reads from a dataset you also require the complement model file).

It is not necessary to generate a model file each time you carry out a new read until experiment. Instead, new files should be obtained whenever a chemistry or pore change occurs. We suggest that it is good practice to extract a model file from the most recent read data available to you when preparing for a read until run as these models have been silently updated in the past.

Ideally the read from which to extract a model should come from the pass folder of a base called run. To extract both template and complement models, this read must be from a 2D base called run. getmodels.py will output whichever models it finds within the file.

To print the getmodels.py help statement at the prompt type:

    C:\path\to\RUscripts>python getmodels.py -h

which will output:

    usage: getmodels.py [-h] [-r READ] [-v] [-ver]

    getmodels: A program to extract model files from Nanopore reads. Ideally you
    should provide a read that has given a 2D read - i.e comes from the pass
    folder on a 2D basecalling run. For Read Until purposes, extraction of just a
    template model is sufficient. Model files will be appended with a number to
    identify the kmer length of the model.

    optional arguments:
    -h, --help            show this help message and exit
    -r READ, --read READ  Provide a read file to extract the current model from.
    -v, --verbose-true    Print detailed messages while processing files to aid
                        in debugging.
    -ver, --version       Print version information and exit.

To run getmodels.py at the prompt ($) type:

    $python getmodels.py -r /path/to/your/read.fast5

We provide a suitable read in the folder exampleread

So:

    C:\path\to\RUscripts>python getmodels.py -r exampleread/llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5

will output:

    Processing exampleread/llssbzms2p35x_lambda11ladderup_1208_1_ch11_file27_strand.fast5
    Looking for template model.
    ***: Model Found
    template model file write completed.
    Kmer length is: 6
    File format is:
    Kmer	Mean	Standard Dev
    Looking for complement model.
    ***: Model Found
    complement model file write completed.
    Kmer length is: 6
    File format is:
    Kmer	Mean	Standard Dev

This will have generated two files - one for the template and one for the complement model.

They will be named:

[readtype]_[model_name]_[kmerlength].model

Where [readtype] is either template or complement, [model_name] is the pore/chemistry type (in this case r7.3_e6_70bps_6mer) and [kmerlen] is the locally calculated length of the kmers used in the model.

These model files are used by all other scripts here. Note that the model files will be overwritten by getmodels.py if it is run again on another read with the same pore/chemistry type and kmer length. This is not a problem.

*** We recommend that you update the model regularly.

# ampbalance.py

Ampbalance is a script which will map entire 2D reads to a reference sequence, identify those likely to give a good 2D read by determining if the template and complement reads overlap with one another appropriately and subsequently copy a defined number of those reads up to a specific depth into a folder for further analysis. This script was originally written to facilitate the sequencing and analysis of Ebola virus and to simplify the transfer of read data prior to base calling where network bandwidth was limited (see <NAME>., <NAME>., <NAME>., <NAME>. & Severi, E. Real-time, portable genome sequencing for Ebola Surveillance. Nature, doi:10.1038/nature16996 (2016).)

Unlike read until applications this script uses the entire read.

Here the script is implemented on reads that have been basecalled to facilitate testing and comparisons.

To print the ampbalance.py help statement at the prompt ($) type:

    C:\path\to\RUscripts>python ampbalance.py -h

which will output:

    usage: ampbalance.py [-h] -fasta FASTA -ids IDS -w WATCHDIR -o TARGETPATH -d
                         DEPTH -procs PROCS [-cautious] [-l LENGTH] -t TEMP_MODEL
                         -c COMP_MODEL [-v] [-ver]

    ampbalance: A program designed to balance amplicons from a specific reference
    sequence post sequencing on ONT minIONs but prebasecalling. Developed by <NAME> @mattloose or <EMAIL> for help! Args that start with
    '--' (eg. --reference_fasta_file) can also be set in a config file
    (/path/to/your/script/amp.config or ) by using .ini or
    .yaml-style syntax (eg. reference_fasta_file=value). If an arg is specified in
    more than one place, then command-line values override config file values
    which override defaults.

    optional arguments:
      -h, --help            show this help message and exit
      -fasta FASTA, --reference_fasta_file FASTA
                            The fasta format file for the reference sequence for
                            your organism.
      -ids IDS, --reference_amplicon_positions IDS
                            A file containing a list of amplicon positions defined
                            for the reference sequence. 1 amplicon per line in the
                            format fasta_sequence_name:start-stop e.g
                            J02459:27-1938
      -w WATCHDIR, --watch-dir WATCHDIR
                            The path to the folder containing the downloads
                            directory with fast5 reads to analyse - e.g.
                            C:\data\minion\downloads (for windows).
      -o TARGETPATH, --output-dir TARGETPATH
                            The path to the destination folder for the
                            preprocessed reads
      -d DEPTH, --depth DEPTH
                            The desired coverage depth for each amplicon. Note
                            this is unlikely to be achieved for each amplicon and
                            should probably be an overestimate of the minimum
                            coverage required.
      -procs PROCS, --proc_num PROCS
                            The number of processors to run this on.
      -cautious, --cautious
                            DTW of long reads on low memory systems can cause
                            unexpected crashes. This option will prevent automatic
                            skipping on any reads over 10,000 events. You can
                            optionally increase this length with the -l parameter.
                            USE WITH CAUTION AS THIS MAY CAUSE A SYSTEM TO CRASH.
      -l LENGTH, --length LENGTH
                            A limit on the length of read that ampbalance will
                            attempt to align using DTW - Long reads can cause
                            problems on low memory systems
      -t TEMP_MODEL, --template_model TEMP_MODEL
                            The appropriate template model file to use
      -c COMP_MODEL, --complement_model COMP_MODEL
                            The appropriate complement model file to use
      -v, --verbose-true    Print detailed messages while processing files.
      -ver, --version       show program's version number and exit

This script is designed to match amplicon sequences of known and approximately uniform length to a reference sequence. The reference sequence should be a single sequence.

We provide an example set of reads using lambda, found in the ampbalancetest folder. A typical usage command would be:

    C:\path\to\RUscripts>python ampbalance.py -fasta J02459.fasta -ids lambda_amplicons.txt -w ampbalancetest -o ampbalanceoutput -d 5 -procs 4 -t template_r7.3_e6_70bps_6mer_6.model -c complement_r7.3_e6_70bps_6mer_6.model -l 3000

This will output the following:

    Reading amplicons
    ******AMP DICTIONARY*******
    <type 'dict'>
    {1: 52, 2: 2065, 3: 4070, 4: 6059, 5: 8012, 6: 10008, 7: 12006, 8: 14011, 9: 16076, 10: 18022, 11: 20053}
    {'DO': 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'HF': 0, 'NH': 0}
    Now we are going to try and open the raw reads and do the same as we have done above...
    {'DO': 0, 1: 0, 2: 0, 3: 1, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'TF': 109, 'HF': 4, 'NH': 0}
    {'DO': 0, 1: 0, 2: 0, 3: 1, 4: 0, 5: 0, 6: 0, 7: 1, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'TF': 108, 'HF': 5, 'NH': 0}
    .
    . (many more lines here)
    .
    {'DO': 0, 1: 10, 2: 10, 3: 10, 4: 10, 5: 9, 6: 10, 7: 10, 8: 10, 9: 10, 10: 10, 11: 10, 'BF': 0, 'TF': 1, 'HF': 110, 'NH': 0}
    {'DO': 0, 1: 10, 2: 10, 3: 10, 4: 10, 5: 10, 6: 10, 7: 10, 8: 10, 9: 10, 10: 10, 11: 10, 'BF': 0, 'TF': 0, 'HF': 110, 'NH': 0}
    Amplicon Read Counts
    Amplicon Number: 1 Reads: 10
    Amplicon Number: 2 Reads: 10
    Amplicon Number: 3 Reads: 10
    Amplicon Number: 4 Reads: 10
    Amplicon Number: 5 Reads: 10
    Amplicon Number: 6 Reads: 10
    Amplicon Number: 7 Reads: 10
    Amplicon Number: 8 Reads: 10
    Amplicon Number: 9 Reads: 10
    Amplicon Number: 10 Reads: 10
    Amplicon Number: 11 Reads: 10
    Copying Amplicon Data
    Amplicon Number 1
    Amplicon Number 2
    Amplicon Number 3
    Amplicon Number 4
    Amplicon Number 5
    Amplicon Number 6
    Amplicon Number 7
    Amplicon Number 8
    Amplicon Number 9
    Amplicon Number 10
    Amplicon Number 11

Key:
DO: Template and Complement Strands don't overlap on amplicons as expected. Read Discarded. Unlikely to generate good 2D.
1...11: Amplicon numbers followed by current coverage counts satisfying highest stringency tests.
BF: Bad Files - these files contain more events than the -length threshold. They are unlikely to be of interest in an amplicon run and are skipped unless the -cautious flag is set.
TF: Total Files - this counts down the number of files left to process.
HF: Hairpin Found - this is the number of reads containing a hairpin - these are tested.
NH: No Hairpin - these reads are ignored as no hairpin has been found.

This program will output files to the specified output directory. Files which are already basecalled will be written to a subfolder called "Downloads". Raw files will just be written to the specified output directory.

Note that the example read set here has been preselected to contain 10 reads mapping to each amplicon and all the reads are 2D - thus the HF value is 110 and the DO value is 0.

It is also important to note that this script uses the hairpin_found flag in the read files. Since minKNOW version 0.51.1.62 (released February 22nd 2016) the writing of this flag to read files has been problematic. Thus reads generated with this version of minKNOW may appear to have far fewer 2D reads than they really do.


# ampliconSPLIT.py

ampliconSPLIT.py simulates read until on either raw or basecalled reads. It ignores the first 50 events of a read and matches the subsequent 250 events against a reference squiggle. Read until only ever processes template data and so only a template model is required.

This script will process reads into subfolders in the targetpath corresponding to each amplicon described in the -ids file - e.g /targetpath/{amplicon_number}. If the read is basecalled, it will be placed in a downloads subfolder - e.g /targetpath/{amplicon_number}/downloads .

To print the ampliconSPLIT.py help statement at the prompt ($) type:

    C:\path\to\RUscripts>python ampliconSPLIT.py -h

which will output:

    usage: ampliconSPLIT.py [-h] -fasta FASTA -ids IDS -w WATCHDIR -o TARGETPATH
                        -d DEPTH -procs PROCS -t TEMP_MODEL [-v] [-ver]

    ampliconSPLIT: A program designed to identify and group individual amplicons
    from minION reads prior to base calling. The depth setting limits the number
    of reads copied to each sub folder. Developed by <NAME> @mattloose or
    <EMAIL> for help! Args that start with '--' (eg.
    --reference_fasta_file) can also be set in a config file
    (/Users/mattloose/fixes/RUscripts/UPDATES/amp.config or ) by using .ini or
    .yaml-style syntax (eg. reference_fasta_file=value). If an arg is specified in
    more than one place, then command-line values override config file values
    which override defaults.

    optional arguments:
    -h, --help            show this help message and exit
    -fasta FASTA, --reference_fasta_file FASTA
                        The fasta format file for the reference sequence for
                        your organism.
    -ids IDS, --reference_amplicon_positions IDS
                        A file containing a list of amplicon positions defined
                        for the reference sequence. 1 amplicon per line in the
                        format fasta_sequence_name:start-stop e.g
                        EM_079517:27-1938
    -w WATCHDIR, --watch-dir WATCHDIR
                        The path to the folder containing the downloads
                        directory with fast5 reads to analyse - e.g.
                        C:\data\minion\downloads (for windows).
    -o TARGETPATH, --output-dir TARGETPATH
                        The path to the destination folder for the
                        preprocessed reads
    -d DEPTH, --depth DEPTH
                        The desired coverage depth for each amplicon. Note
                        this is unlikely to be achieved for each amplicon and
                        should probably be an overestimate of the minimum
                        coverage required.
    -procs PROCS, --proc_num PROCS
                        The number of processors to run this on.
    -t TEMP_MODEL, --template_model TEMP_MODEL
                        The appropriate template model file to use. This file
                        can be generated uing the getmodels.py script.
    -v, --verbose-true    Print detailed messages while processing files.
    -ver, --version       show program's version number and exit

An example run command using test data:

    C:\path\to\RUscripts>python ampliconSPLIT.py -fasta J02459.fasta -ids lambda_amplicons.txt -w RUtestset/ -o test -d 10 -procs 8 -t template_r7.3_e6_70bps_6mer_6.model

Which will output:

    processing the reference fasta.
    ID J02459
    length 48502
    FORWARD STRAND
    REVERSE STRAND
    Groking amplicons
    ******AMP DICTIONARY*******
    <type 'dict'>
    {1: 52, 2: 2065, 3: 4070, 4: 6059, 5: 8012, 6: 10008, 7: 12006, 8: 14011, 9: 16076, 10: 18022, 11: 20053}
    {'DO': 0, 1: 0, 2: 0, 3: 0, 4: 0, 5: 0, 6: 0, 7: 0, 8: 0, 9: 0, 10: 0, 11: 0, 'BF': 0, 'HF': 0, 'NH': 0}
    We want to build a custom reference that is smaller than the original reference.
    First get a list of all the positions we will need to search
    Generating a custom fasta
    J02459
    processing the custom fasta
    Attempting to match reads and split into folders based on 250 events, excluding the first 50.
    Amplicon Read Counts
        Amplicon Number: 1 Reads: 4
        Amplicon Number: 2 Reads: 5
        Amplicon Number: 3 Reads: 4
        Amplicon Number: 4 Reads: 5
        Amplicon Number: 5 Reads: 6
        Amplicon Number: 6 Reads: 6
        Amplicon Number: 7 Reads: 5
        Amplicon Number: 8 Reads: 5
        Amplicon Number: 9 Reads: 5
        Amplicon Number: 10 Reads: 5
        Amplicon Number: 11 Reads: 5
    Copying Amplicon Data
    Amplicon Number 1
    Amplicon Number 2
    Amplicon Number 3
    Amplicon Number 4
    Amplicon Number 5
    Amplicon Number 6
    Amplicon Number 7
    Amplicon Number 8
    Amplicon Number 9
    Amplicon Number 10
    Amplicon Number 11

The depth parameter (-d) sets the number of reads that will be copied. Setting this to a value greater than the number of reads analysed will sort all reads in the dataset.


# test_gReadUntil.py

*** Note this script is located with the ReadUntil folder of this repository

This script runs entirely independently of the read until API and allows for simulation of selective sequencing of a genome. You can provide a pool of reads and the script will use the first 250 events to map a read and copy those mapping within the desired area to a specified folder.

Help is available by typing:

    C:\path\to\RUscripts\ReadUntil\python test_gReadUntil.py -h

which will output:

    usage: test_gReadUntil.py [-h] -fasta FASTA -targets [TARGETS [TARGETS ...]]
                          -procs PROCS -m TEMP_MODEL [-log LOGFILE] -w
                          WATCHDIR [-o OUTPUT_FOLDER] [-v] [-ver]

    real_read_until: A program providing read until with the Oxford Nanopore
    minION device. This program will ultimately be driven by minoTour to enable
    selective remote sequencing. This program is heavily based on original code
    generously provided by Oxford Nanopore Technologies.

    optional arguments:
    -h, --help            show this help message and exit
    -fasta FASTA, --reference_fasta_file FASTA
                        The fasta format file describing the reference
                        sequence for your organism.
    -targets [TARGETS [TARGETS ...]]
                        Positional IDs to enrich for in the form seqid:start-
                        stop . Can be space seperated eg: J02459:10000-15000
                        J02459:35000-40000
    -procs PROCS, --proc_num PROCS
                        The number of processors to run this on.
    -m TEMP_MODEL, --model TEMP_MODEL
                        The appropriate template model file to use
    -log LOGFILE, --log-file LOGFILE
                        The name of the log file that data will be written to
                        regarding the decision made by this program to process
                        read until.
    -w WATCHDIR, --watch-dir WATCHDIR
                        The path to the folder containing the downloads
                        directory with fast5 reads to analyse - e.g.
                        C:\data\minion\downloads (for windows).
    -o OUTPUT_FOLDER, --output OUTPUT_FOLDER
                        Path to a folder to symbolically place reads
                        representing match and not match.
    -v, --verbose-true    Print detailed messages while processing files.
    -ver, --version       show program's version number and exit

A typical command line to select reads mapping from 10-15kb in the lambda genome would be:

    C:\path\to\RUscripts\ReadUntil\python test_gReadUntil.py -fasta ../J02459.fasta -targets J02459:10000-15000 -procs 4  -m ../template_r7.3_e6_70bps_6mer_6.model -w ../RUtestset/ -o test

This would give the following output:

    ***********************************************************************************************
    **** This code will open a collection of reads and simulate read until on them. It will    ****
    **** copy reads into a secondary folder for subsequent processing by another analysis      ****
    **** package.                                                                              ****
    ***********************************************************************************************
    processing the reference fasta.
    J02459:10000-15000
    ID J02459
    length 48502
    1 J02459:10000-15000 10000 15000 J02459
    We want to extract this chunk J02459_1
    ID J02459_1
    length 5000
    FORWARD STRAND
    REVERSE STRAND
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch10_file49_strand.fast5 No Match
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch10_file69_strand.fast5 Sequence Found
    .
    .
    .
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch8_file58_strand.fast5 Sequence Found
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch9_file100_strand.fast5 No Match
    ../RUtestset/llssbzms2p35x_20151004_readuntiludududududu_RU21_lambdaPCR_2922_1_ch8_file11_strand.fast5 No Match

Sequence found indicates that the read is derived from the desired region. No Match indicates that the read is from another region.


# Real Read Until

Important Note:

Running read until will influence the behaviour of your flow cell and change the output of your sequencing experiment. You are strongly advised to run simulations of read until prior to running on a live flow cell. The code as presented here is a demonstration of read until and one implementation. Users run this code entirely at their own risk.


We provide the two scripts used to generate the data presented in "Real time selective sequencing using nanopore technology".

These scripts absolutely require the read until API from Oxford Nanopore.

# Read Until API

This version of Read Until requires that the user runs a command window as Administrator.

The read until API changes depending on the version of minKNOW in use. The read until API for minKNOW versions 0.48.1.3 - 0.51.1.51 is available on request from Oxford Nanopore and the code provided here is compatible with Git Commit

    <PASSWORD>

and minKNOW versions <= 0.51.1.51

Note that this is not the latest commit. The latest version as of 30th March 2016 does not function as anticipated. It will work for simulation runs with ws_event_sample, but minKNOW itself does not function correctly with the API. This is actively being resolved by ONT. Thus, testing and simulation is possible with our code and later versions of minKNOW and ws_event_sampler but these scripts WILL NOT WORK AT THIS TIME with minKNOW > 0.51.1.51

# Installing the API

Copy the following files and folders (and there contents) into the RUscripts\ReadUntil folder:

    read_until.py
	event_sampler.thrift
    event_sampler/* (note - you need to copy the folder and its contents)
    example.py

# Running ws_event_sampler

Read until works by streaming data from minKNOW during a run to an external analysis pipeline and receiving instructions to reject reads as required. The program streaming the data is called ws_event_sampler.exe and this is found in the folder

    C:\grouper\binaries\bin

which is available after installing minKNOW.

ws_event_sampler can run either stand alone in simulation mode, or alongside minKNOW as it is running. For simulations, we suggest you start minKNOW and have it open (but not actually sequencing). Then start ws_event_sampler in simluation mode. In this way you will see the output from our scripts in the minKNOW messages window.

Running ws_event_sampler is done as follows:

    C:\grouper\binaries>bin\ws_event_sampler.exe

*** NOTE: The command itself is bin\ws_event_sampler.exe so the command is executed from the folder binaries, not from within the folder bin. Executing the file from within the bin folder will result in the program crashing with a message about being unable to open a config file.

*** NOTE: ws_event_sampler provides no visible indication that it is running other than blocking the command line. If it fails, it will terminate - sometimes with an error message, but often not. Support for ws_event_sampler is available direct from ONT. We hope we have provided sufficient guidance for its use in the context of our scripts here.

ws_event_sampler can be run in simulation mode to facilitate the development of analysis scripts. When running in simulation mode, ws_event_sampler provides reads driven from a model file stored locally. This same model file must be used to analyse data when in simulation mode. When streaming live data direct from minKnow it is essential to use the correct model for the chemistry and pore type in use.

The read until API provides comprehensive instructions on how to run ws_event_sampler. Below we provide worked examples for our scripts.

To test that ws_event_sampler is functioning correctly, do the following:

1) Open a command window as administrator - to do this, right click on the cmd icon and select "Run as Administrator"

2) In this window navigate to C:\grouper\binaries

3) To see help for ws_event_sampler type:

    C:\grouper\binaries>bin\ws_event_sampler.exe -h

which will output:

    Usage: ws_event_sampler [Options]

    Common options:
      -h [ --help ]                         Produce help message
      -c [ --config ] arg (=./conf/global.conf)
                                            Application config file.

    Specific options:
      -p [ --port ] arg                     Port to listen for requests. If this is
                                            not specified it is taken from the
                                            ws_event_sampler_port field in the
                                            config file.
      -s [ --sim ]                          Start in standalone simulation mode.
      --sim-channels arg                    Number of channels simulation uses.
                                            E.g. use a low number when debugging to
                                            make the pushed data manageable. This
                                            overrides the channel_count setting in
                                            the conf file but only during
                                            simulation.
      --sim-model arg                       Path to model file to be used by
                                            simulation when generating events.
      --sim-fasta arg                       Path to fasta file from which to draw
                                            fragments to be used by simulation.
      --sim-delay arg (=0)                  Seconds delay before starting
                                            simulation. By default the simulation
                                            starts as soon as this server starts
                                            up.
      --sim-fragment-length arg             Random lengths are bounded, approximate
                                            normal distributions, picked by rolling
                                            lots of dice. You need to specify die
                                            sides, die count and a base amount. For
                                            example, the default setting of '5 250
                                            500' will generate reads with mean
                                            length 1000 and bounds [500,1500].
      --sim-log arg                         Filename to use for csv-format log of
                                            the simulation. Will overwrite if
                                            already exists.

To set up a simulation on a specific port type:

    C:\grouper\binaries>bin\ws_event_sampler.exe -s --sim-channels 100 --sim-fasta C:\path\to\RUscripts\J02459.fasta --sim-log log.txt

This will establish the data stream on port 12345, simulating 100 channels with a read distribution as described derived from the lambda sequence we have provided here. ws_event_sampler in simulation mode will write out a log file to log.txt enabling tracking of events.

In the example.py script, edit the line which states:

    host="ws://localhost:<port>"

to

    host="ws://localhost:9200"

Note - you can specify a different port value if you wish.

Then from another command window execute the examply.py script (note this is supplied by the ONT API)

    C:\path\to\RUscripts\ReadUntil>python example.py

This should output something like (depending on the version of the API in use):

    Client connection started. Beginning unblock loop...
    channel_name= 344  read_number= 12  events_in_sample= 200  first event: start= 1104931.0  mean= 71.1864189873
    channel_name= 345  read_number= 12  events_in_sample= 200  first event: start= 1111035.0  mean= 68.1394555562
    channel_name= 346  read_number= 12  events_in_sample= 200  first event: start= 1104312.0  mean= 67.1975140682
    channel_name= 347  read_number= 15  events_in_sample= 200  first event: start= 1203518.0  mean= 72.436242832
    channel_name= 340  read_number= 12  events_in_sample= 200  first event: start= 1104326.0  mean= 60.9232571465
    channel_name= 341  read_number= 12  events_in_sample= 200  first event: start= 1102520.0  mean= 69.496870232
    channel_name= 342  read_number= 12  events_in_sample= 200  first event: start= 1097286.0  mean= 65.2731245584
    channel_name= 343  read_number= 12  events_in_sample= 200  first event: start= 1102224.0  mean= 65.2680067692
    channel_name= 348  read_number= 12  events_in_sample= 200  first event: start= 1119138.0  mean= 77.3645669056
    channel_name= 349  read_number= 12  events_in_sample= 200  first event: start= 1092464.0  mean= 72.4693304753
    channel_name= 298  read_number= 12  events_in_sample= 200  first event: start= 1109404.0  mean= 66.3816659282
    channel_name= 299  read_number= 12  events_in_sample= 200  first event: start= 1103390.0  mean= 66.1452746381

If you see:

    Client connection started. Beginning unblock loop...
    Unhandled exception in thread started by <bound method ReadUntil.data_received_own_thread of <read_until.ReadUntil object at 0x00000000023AA
    FD0>>
{'fraction_non_alphanumeric': 0.05657595416631561, 'fraction_numerical': 0.0391917319628163, 'mean_word_length': 2.596535528508661, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 74, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 16, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 122, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1805272', 'n_tokens_mistral': 21426, 'n_tokens_neox': 18518, 'n_words': 9287}