Commit 14a7251d authored by Andreas Tille's avatar Andreas Tille

Imported Upstream version 0.0+20110617

parents
Copyright (c) 2016, Vladimir Boza, Comenius University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Comenius University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL COMENIUS UNIVERSITY BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
### DeepNano: alternative basecaller for MinION reads
DeepNano is alternative basecaller for Oxford Nanopore MinION reads
based on deep recurrent neural networks.
Currently it works with SQK-MAP-006 and SQK-MAP-005 chemistry and as a postprocessor for Metrichor.
Here are our benchmarks for SQK-MAP-006 chemistry, which compare mapping accuracy (we trained on reads which align to one half on the
Ecoli and tested on other half of Ecoli and Klebsiela):
| | Ecoli Metrichor | Ecoli DeepNano | Klebsiella Metrichor | Klebsiella DeepNano |
|------------------|-----------------|----------------|----------------------|---------------------|
| Template reads | 71.3% | 77.9% | 68.1% | 76.3% |
| Complement reads | 71.4% | 76.4% | 69.5% | 75.7% |
| 2D reads | 86.8% | 88.5% | 84.8% | 86.7% |
Links to datasets with reads:
- http://www.ebi.ac.uk/ena/data/view/ERR1147230
- https://www.ebi.ac.uk/ena/data/view/SAMEA3713789
Requirements
================
We use Python 2.7.
Here are versions of Python packages, that we used:
- Cython==0.23.4
- numpy==1.10.2
- h5py==2.5.0
- Theano==0.7.0
- python-dateutil==2.5.0
Basic usage:
================
For SQK-MAP-006 chemistry just use:
`OMP_NUM_THREADS=1 python basecall.py <list of fast5 files>`
or
`OMP_NUM_THREADS=1 python basecall.py --directory <directory with reads>`
It outputs basecalls for template, complement and 2D into file named output.fasta.
For SQK-MAP-005 chemistry use:
`OMP_NUM_THREADS=1 python basecall.py --template_net nets_data/map5temp.npz --complement_net nets_data/map5comp.npz --big_net nets_data/map5-2d.npz <list of fast5 files>`
Advanced arguments:
=================
- `-h` - prints help message
- `--template_net PATH` - path to network which basecalls template (has reasonable default)
- `--complement_net PATH` - path to network which basecalls complement (has reasonable default)
- `--big_net PATH` - path to network which basecalls 2D (has reasonable default)
- `--timing` - if set, display timing information for each read
- `--type template/complement/2d/all` - type of basecalling output (defaults to all)
- `--output FILENAME` - output filename
- `--output_orig` - if set, outputs also Metrichor basecalls
- `--directory DIRECTORY` Directory where read files are stored
Usage without metrichor:
================
Basecaller above used spliting to template and
complement, scaling parameters and alignment from metrichor.
The tools below does not require metrichor at some expense at accuracy.
Currently the accuracy for 1D basecall is almost similar as for the tools above
(except for some reads where we really mess up scaling).
The accuracy for 2D basecall is approximatelly 0.5% lower than for metrichor (whereas
the tool above was 2% better than metrichor).
To use this tool, first compile alignment code:
`g++ -O2 -std=gnu++0x align_2d.cc -o align_2d`
Then run:
`OMP_NUM_THREADS=1 python basecall_no_metrichor.py <list of fast5 files>`
Arguments:
- `-h` - prints help message
- `--template_net PATH` - path to network which basecalls template (has reasonable default)
- `--complement_net PATH` - path to network which basecalls complement (has reasonable default)
- `--big_net PATH` - path to network which basecalls 2D (has reasonable default)
- `--type template/complement/2d/all` - type of basecalling output (defaults to all)
- `--output FILENAME` - output filename
- `--directory DIRECTORY` Directory where read files are stored
If you want to watch a directory for new files, first install:
- watchdog==0.8.3
And then use (the output parameter has no effect, we output separate fasta files for each fast5 file):
`OMP_NUM_THREADS=1 python basecall_no_metrichor.py --watch <directory name>`
#include <cstdio>
#include <vector>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;
struct Prob {
double p[4];
double n;
};
Prob LoadProb(int parity) {
Prob p;
scanf("%lf %lf %lf %lf %lf", &p.p[0], &p.p[1], &p.p[2], &p.p[3], &p.n);
for (int i = 0; i < 4; i++) {
p.p[i] = log(p.p[i]);
}
p.n = log(p.n);
return p;
}
Prob LoadProbC(int parity) {
Prob p;
scanf("%lf %lf %lf %lf %lf", &p.p[3], &p.p[2], &p.p[1], &p.p[0], &p.n);
for (int i = 0; i < 4; i++) {
p.p[i] = log(p.p[i]);
}
p.n = log(p.n);
return p;
}
int main() {
vector<Prob> temp, comp;
int tc, cc;
scanf("%d", &tc);
for (int i = 0; i < tc; i++) {
temp.push_back(LoadProb(i%2));
}
for (int i = 0; i < tc; i+=2) {
temp[i].n -= log(0.9) / 2;
}
scanf("%d", &cc);
for (int i = 0; i < cc; i++) {
comp.push_back(LoadProbC(i%2));
}
for (int i = 0; i < cc; i+=2) {
comp[i].n -= log(0.9) / 2;
}
reverse(comp.begin(), comp.end());
vector<vector<double>> probs(temp.size()+1, vector<double>(comp.size()+1));
// 0 - 3 ACGT
// 4 - Ncomp
// 5 - Ntemp
vector<vector<int>> prevs(temp.size()+1, vector<int>(comp.size()+1));
double minf = -1000000000000.0;
double pen = -0.5;
probs[0][0] = 0;
for (int i = 0; i < probs.size(); i++) {
for (int j = 0; j < probs[i].size(); j++) {
if (i == 0 && j == 0) continue;
probs[i][j] = 0;
// if (i == 0 && (j < 2000 && j < probs[i].size() / 4)) probs[i][j] = 0;
// if (j == 0 && (i < 2000 && i < probs.size() / 4)) probs[i][j] = 0;
prevs[i][j] = 6;
if (j > 0) {
double np = probs[i][j-1] + comp[j-1].n - pen;
if (np > probs[i][j]) {
prevs[i][j] = 4;
probs[i][j] = np;
}
}
if (i > 0) {
double np = probs[i-1][j] + temp[i-1].n - pen;
if (np > probs[i][j]) {
prevs[i][j] = 5;
probs[i][j] = np;
}
}
if (i > 0 && j > 0) {
for (int k = 0; k < 4; k++) {
double np = probs[i-1][j-1] + (temp[i-1].p[k] + comp[j-1].p[k]) - 2*pen;
if (np > probs[i][j]) {
prevs[i][j] = k;
probs[i][j] = np;
}
}
}
}
}
fprintf(stderr, "%lf\n", probs.back().back());
char alph[] = "ACGT";
string seq;
int ipos = temp.size();
int jpos = comp.size();
/* int margin = min(2000, (int)temp.size() / 4);
for (int i = temp.size(); i >= temp.size() - margin && i >= 0; i--) {
if (probs[i][comp.size()] > probs[ipos][jpos]) {
ipos = i;
jpos = comp.size();
}
}
margin = min(2000, (int)comp.size() / 4);
for (int j = comp.size(); j >= comp.size() - margin && j >= 0; j--) {
if (probs[temp.size()][j] > probs[ipos][jpos]) {
ipos = temp.size();
jpos = j;
}
}*/
for (int i = 0; i < temp.size(); i++) {
for (int j = 0; j < comp.size(); j++) {
if (probs[i][j] > probs[ipos][jpos]) {
ipos = i;
jpos = j;
}
}
}
vector<pair<int, int>> trace;
while (ipos > 0 && jpos > 0) {
if (prevs[ipos][jpos] == 6) {
break;
}
trace.push_back(make_pair(ipos, jpos));
if (prevs[ipos][jpos] == 4) {
jpos--;
} else if (prevs[ipos][jpos] == 5) {
ipos--;
} else {
seq += alph[prevs[ipos][jpos]];
ipos--;
jpos--;
}
}
reverse(trace.begin(), trace.end());
fprintf(stderr, "%d\n", seq.size());
reverse(seq.begin(), seq.end());
printf("%s\n", seq.c_str());
int last_temp = -47;
int last_comp = -47;
for (int i = 10; i + 10 < trace.size(); i++) {
auto t = trace[i];
int temp_al = -1;
int comp_al = -1;
if (t.first != last_temp && t.first % 2 == 0) {
temp_al = (t.first - 1) / 2;
}
if (t.second != last_comp && t.second % 2 == 1) {
comp_al = comp.size() / 2 - 2 - t.second / 2;
}
if (temp_al != -1 || comp_al != -1) {
printf("%d %d\n", temp_al, comp_al);
}
last_temp = t.first;
last_comp = t.second;
}
}
import argparse
from rnn_fin import RnnPredictor
import h5py
import sys
import numpy as np
import theano as th
import os
import re
import dateutil.parser
import datetime
from helpers import *
def load_read_data(read_file):
h5 = h5py.File(read_file, "r")
ret = {}
extract_timing(h5, ret)
base_loc = get_base_loc(h5)
try:
ret["called_template"] = h5[base_loc+"/BaseCalled_template/Fastq"][()].split('\n')[1]
ret["called_complement"] = h5[base_loc+"/BaseCalled_complement/Fastq"][()].split('\n')[1]
ret["called_2d"] = h5["Analyses/Basecall_2D_000/BaseCalled_2D/Fastq"][()].split('\n')[1]
except Exception as e:
pass
try:
events = h5[base_loc+"/BaseCalled_template/Events"]
tscale, tscale_sd, tshift, tdrift = extract_scaling(h5, "template", base_loc)
ret["temp_events"] = extract_1d_event_data(
h5, "template", base_loc, tscale, tscale_sd, tshift, tdrift)
except:
pass
try:
cscale, cscale_sd, cshift, cdrift = extract_scaling(h5, "complement", base_loc)
ret["comp_events"] = extract_1d_event_data(
h5, "complement", base_loc, cscale, cscale_sd, cshift, cdrift)
except Exception as e:
pass
try:
al = h5["Analyses/Basecall_2D_000/BaseCalled_2D/Alignment"]
temp_events = h5[base_loc+"/BaseCalled_template/Events"]
comp_events = h5[base_loc+"/BaseCalled_complement/Events"]
ret["2d_events"] = []
for a in al:
ev = []
if a[0] == -1:
ev += [0, 0, 0, 0, 0]
else:
e = temp_events[a[0]]
mean = (e["mean"] - tshift) / cscale
stdv = e["stdv"] / tscale_sd
length = e["length"]
ev += [1] + preproc_event(mean, stdv, length)
if a[1] == -1:
ev += [0, 0, 0, 0, 0]
else:
e = comp_events[a[1]]
mean = (e["mean"] - cshift) / cscale
stdv = e["stdv"] / cscale_sd
length = e["length"]
ev += [1] + preproc_event(mean, stdv, length)
ret["2d_events"].append(ev)
ret["2d_events"] = np.array(ret["2d_events"], dtype=np.float32)
except Exception as e:
print e
pass
h5.close()
return ret
parser = argparse.ArgumentParser()
parser.add_argument('--template_net', type=str, default="nets_data/map6temp.npz")
parser.add_argument('--complement_net', type=str, default="nets_data/map6comp.npz")
parser.add_argument('--big_net', type=str, default="nets_data/map6-2d-big.npz")
parser.add_argument('reads', type=str, nargs='*')
parser.add_argument('--timing', action='store_true', default=False)
parser.add_argument('--type', type=str, default="all", help="One of: template, complement, 2d, all, use comma to separate multiple options, eg.: template,complement")
parser.add_argument('--output', type=str, default="output.fasta")
parser.add_argument('--output_orig', action='store_true', default=False)
parser.add_argument('--directory', type=str, default='', help="Directory where read files are stored")
args = parser.parse_args()
types = args.type.split(',')
do_template = False
do_complement = False
do_2d = False
if "all" in types or "template" in types:
do_template = True
if "all" in types or "complement" in types:
do_complement = True
if "all" in types or "2d" in types:
do_2d = True
assert do_template or do_complement or do_2d, "Nothing to do"
assert len(args.reads) != 0 or len(args.directory) != 0, "Nothing to basecall"
if do_template:
print "loading template net"
temp_net = RnnPredictor(args.template_net)
print "done"
if do_complement:
print "loading complement net"
comp_net = RnnPredictor(args.complement_net)
print "done"
if do_2d:
print "loading 2D net"
big_net = RnnPredictor(args.big_net)
print "done"
chars = "ACGT"
mapping = {"A": 0, "C": 1, "G": 2, "T": 3, "N": 4}
fo = open(args.output, "w")
total_bases = [0, 0, 0]
files = args.reads
if len(args.directory):
files += [os.path.join(args.directory, x) for x in os.listdir(args.directory)]
for i, read in enumerate(files):
basename = os.path.basename(read)
try:
data = load_read_data(read)
except Exception as e:
print "error at file", read
print e
continue
if not data:
continue
print "\rcalling read %d/%d %s" % (i, len(files), read),
sys.stdout.flush()
if args.output_orig:
try:
if "called_template" in data:
print >>fo, ">%s_template" % basename
print >>fo, data["called_template"]
if "called_complement" in data:
print >>fo, ">%s_complement" % basename
print >>fo, data["called_complement"]
if "called_2d" in data:
print >>fo, ">%s_2d" % basename
print >>fo, data["called_2d"]
except:
pass
temp_start = datetime.datetime.now()
if do_template and "temp_events" in data:
predict_and_write(data["temp_events"], temp_net, fo, "%s_template_rnn" % basename)
temp_time = datetime.datetime.now() - temp_start
comp_start = datetime.datetime.now()
if do_complement and "comp_events" in data:
predict_and_write(data["comp_events"], comp_net, fo, "%s_complement_rnn" % basename)
comp_time = datetime.datetime.now() - comp_start
start_2d = datetime.datetime.now()
if do_2d and "2d_events" in data:
predict_and_write(data["2d_events"], big_net, fo, "%s_2d_rnn" % basename)
time_2d = datetime.datetime.now() - start_2d
if args.timing:
try:
print "Events: %d/%d" % (len(data["temp_events"]), len(data["comp_events"]))
print "Our times: %f/%f/%f" % (temp_time.total_seconds(), comp_time.total_seconds(),
time_2d.total_seconds())
print "Our times per base: %f/%f/%f" % (
temp_time.total_seconds() / len(data["temp_events"]),
comp_time.total_seconds() / len(data["comp_events"]),
time_2d.total_seconds() / (len(data["comp_events"]) + len(data["temp_events"])))
print "Their times: %f/%f/%f" % (data["temp_time"].total_seconds(), data["comp_time"].total_seconds(), data["2d_time"].total_seconds())
print "Their times per base: %f/%f/%f" % (
data["temp_time"].total_seconds() / len(data["temp_events"]),
data["comp_time"].total_seconds() / len(data["comp_events"]),
data["2d_time"].total_seconds() / (len(data["comp_events"]) + len(data["temp_events"])))
except:
# Don't let timing throw us out
pass
fo.flush()
fo.close()
import argparse
from rnn_fin import RnnPredictor
import h5py
import sys
import numpy as np
import theano as th
import os
import re
import dateutil.parser
import datetime
from helpers import *
import subprocess
import time
def get_scaling_template(events, has_std):
down = 48.4631279889
up = 65.7312554591
our_down = np.percentile(events["mean"], 10)
our_up = np.percentile(events["mean"], 90)
scale = (our_up - our_down) / (up - down)
shift = (our_up / scale - up) * scale
sd = 0.807981325017
if has_std:
return scale, np.percentile(events["stdv"], 50) / sd, shift
else:
return scale, np.sqrt(np.percentile(events["variance"], 50)) / sd, shift
def get_scaling_complement(events, has_std):
down = 49.2638926877
up = 69.0192568072
our_down = np.percentile(events["mean"], 10)
our_up = np.percentile(events["mean"], 90)
scale = (our_up - our_down) / (up - down)
shift = (our_up / scale - up) * scale
sd = 1.04324844612
if has_std:
return scale, np.percentile(events["stdv"], 50) / sd, shift
else:
return scale, np.sqrt(np.percentile(events["variance"], 50)) / sd, shift
def template_complement_loc(events):
abasic_level = np.percentile(events["mean"], 99) + 5
abasic_locs = (events["mean"] > abasic_level).nonzero()[0]
last = -47
run_len = 1
runs = []
for x in abasic_locs:
if x - last == 1:
run_len += 1
else:
if run_len >= 5:
if len(runs) and last - runs[-1][0] < 50:
run_len = last - runs[-1][0]
run_len += runs[-1][1]
runs[-1] = (last, run_len)
else:
runs.append((last, run_len))
run_len = 1
last = x
to_sort = []
mid = len(events) / 2
low_third = len(events) / 3
high_third = len(events) / 3 * 2
for r in runs:
if r[0] < low_third:
continue
if r[0] > high_third:
continue
to_sort.append((abs(r[0] - mid), r[0] - r[1], r[0]))
to_sort.sort()
if len(to_sort) == 0:
return None
trim_size = 10
return {"temp": (trim_size, to_sort[0][1] - trim_size),
"comp": (to_sort[0][2] + trim_size, len(events) - trim_size)}
def load_read_data(read_file):
h5 = h5py.File(read_file, "r")
ret = {}
read_key = h5["Analyses/EventDetection_000/Reads"].keys()[0]
base_events = h5["Analyses/EventDetection_000/Reads"][read_key]["Events"]
temp_comp_loc = template_complement_loc(base_events)
sampling_rate = h5["UniqueGlobalKey/channel_id"].attrs["sampling_rate"]
if temp_comp_loc:
events = base_events[temp_comp_loc["temp"][0]:temp_comp_loc["temp"][1]]
else:
events = base_events
has_std = True
try:
std = events[0]["stdv"]
except:
has_std = False
tscale2, tscale_sd2, tshift2 = get_scaling_template(events, has_std)
index = 0.0
ret["temp_events2"] = []
for e in events:
mean = (e["mean"] - tshift2) / tscale2
if has_std:
stdv = e["stdv"] / tscale_sd2
else:
stdv = np.sqrt(e["variance"]) / tscale_sd2
length = e["length"] / sampling_rate
ret["temp_events2"].append(preproc_event(mean, stdv, length))
ret["temp_events2"] = np.array(ret["temp_events2"], dtype=np.float32)
if not temp_comp_loc:
return ret
events = base_events[temp_comp_loc["comp"][0]:temp_comp_loc["comp"][1]]
cscale2, cscale_sd2, cshift2 = get_scaling_complement(events, has_std)
index = 0.0
ret["comp_events2"] = []
for e in events:
mean = (e["mean"] - cshift2) / cscale2
if has_std:
stdv = e["stdv"] / cscale_sd2
else:
stdv = np.sqrt(e["variance"]) / cscale_sd2
length = e["length"] / sampling_rate
ret["comp_events2"].append(preproc_event(mean, stdv, length))
ret["comp_events2"] = np.array(ret["comp_events2"], dtype=np.float32)
return ret
def basecall(read_file_name, fo):
basename = os.path.basename(read_file_name)
try:
data = load_read_data(read_file_name)
except Exception as e:
print e
print "error at file", read_file_name
return
if do_template or do_2d:
o1, o2 = predict_and_write(
data["temp_events2"], temp_net,
fo if do_template else None,
"%s_template_rnn" % basename)
if (do_complement or do_2d) and "comp_events2" in data:
o1c, o2c = predict_and_write(
data["comp_events2"], comp_net,
fo if do_complement else None,
"%s_complement_rnn" % basename)
if do_2d and "comp_events2" in data and\
len(data["comp_events2"]) <= args.max_2d_length and\
len(data["temp_events2"]) <= args.max_2d_length:
p = subprocess.Popen("./align_2d", stdin=subprocess.PIPE, stdout=subprocess.PIPE)
f2d = p.stdin
print >>f2d, len(o1)+len(o2)
for a, b in zip(o1, o2):
print >>f2d, " ".join(map(str, a))
print >>f2d, " ".join(map(str, b))
print >>f2d, len(o1c)+len(o2c)
for a, b in zip(o1c, o2c):
print >>f2d, " ".join(map(str, a))
print >>f2d, " ".join(map(str, b))
f2do, f2de = p.communicate()
if p.returncode != 0:
return
lines = f2do.strip().split('\n')
print >>fo, ">%s_2d_rnn_simple" % basename
print >>fo, lines[0].strip()
events_2d = []
for l in lines[1:]:
temp_ind, comp_ind = map(int, l.strip().split())
e = []
if temp_ind == -1:
e += [0, 0, 0, 0, 0]
else:
e += [1] + list(data["temp_events2"][temp_ind])
if comp_ind == -1:
e += [0, 0, 0, 0, 0]
else:
e += [1] + list(data["comp_events2"][comp_ind])
events_2d.append(e)
events_2d = np.array(events_2d, dtype=np.float32)
predict_and_write(events_2d, big_net, fo, "%s_2d_rnn" % basename)
parser = argparse.ArgumentParser()
parser.add_argument('--template_net', type=str, default="nets_data/map6temp.npz")
parser.add_argument('--complement_net', type=str, default="nets_data/map6comp.npz")
parser.add_argument('--big_net', type=str, default="nets_data/map6-2d-no-metr23.npz")
parser.add_argument('--max_2d_length', type=int, default=10000, help='Max length for 2d basecall')
parser.add_argument('reads', type=str, nargs='*')
parser.add_argument('--type', type=str, default="all", help="One of: template, complement, 2d, all, use comma to separate multiple options, eg.: template,complement")
parser.add_argument('--output', type=str, default="output.fasta")
parser.add_argument('--directory', type=str, default='', help="Directory where read files are stored")
parser.add_argument('--watch', type=str, default='', help='Watched directory')
args = parser.parse_args()
types = args.type.split(',')
do_template = False
do_complement = False
do_2d = False
if "all" in types or "template" in types:
do_template = True
if "all" in types or "complement" in types:
do_complement = True
if "all" in types or "2d" in types:
do_2d = True
assert do_template or do_complement or do_2d, "Nothing to do"
assert len(args.reads) != 0 or len(args.directory) != 0 or len(args.watch) != 0, "Nothing to basecall"
if do_template or do_2d:
print "loading template net"
temp_net = RnnPredictor(args.template_net)
print "done"
if do_complement or do_2d:
print "loading complement net"
comp_net = RnnPredictor(args.complement_net)
print "done"
if do_2d:
print "loading 2D net"
big_net = RnnPredictor(args.big_net)
print "done"
chars = "ACGT"
mapping = {"A": 0, "C": 1, "G": 2, "T": 3, "N": 4}
if len(args.reads) or len(args.directory) != 0:
fo = open(args.output, "w")
files = args.reads
if len(args.directory):
files += [os.path.join(args.directory, x) for x in os.listdir(args.directory)]
for i, read in enumerate(files):