[med-svn] [nanopolish] 01/03: New upstream version 0.8.1
Steffen Möller
moeller at moszumanska.debian.org
Sat Sep 16 00:14:55 UTC 2017
This is an automated email from the git hooks/post-receive script.
moeller pushed a commit to branch master
in repository nanopolish.
commit 204a942588af8a379ea261b3cf22f54e9dc903a4
Author: Steffen Moeller <moeller at debian.org>
Date: Thu Sep 14 16:58:39 2017 +0200
New upstream version 0.8.1
---
Makefile | 6 +-
README.md | 28 +-
scripts/nanopolish_makerange.py | 15 +-
src/alignment/nanopolish_alignment_db.cpp | 15 +-
src/alignment/nanopolish_alignment_db.h | 2 +-
src/alignment/nanopolish_anchor.cpp | 223 ------
src/alignment/nanopolish_anchor.h | 15 -
src/alignment/nanopolish_eventalign.cpp | 12 +-
src/common/{fs_support.hpp => fs_support.cpp} | 28 +-
src/common/fs_support.hpp | 34 +-
src/common/nanopolish_common.h | 2 +-
src/common/nanopolish_fast5_map.cpp | 138 ----
src/common/nanopolish_fast5_map.h | 39 -
src/common/nanopolish_variant.cpp | 13 +-
src/main/nanopolish.cpp | 7 +-
src/nanopolish_call_methylation.cpp | 12 +-
src/nanopolish_call_variants.cpp | 15 +-
src/nanopolish_consensus.cpp | 897 ----------------------
src/nanopolish_consensus.h | 14 -
src/nanopolish_extract.cpp | 39 +-
src/nanopolish_getmodel.cpp | 4 +-
src/nanopolish_haplotype.cpp | 20 +
src/nanopolish_haplotype.h | 5 +
src/nanopolish_index.cpp | 197 +++++
src/nanopolish_index.h | 14 +
src/nanopolish_methyltrain.cpp | 16 +-
src/nanopolish_phase_reads.cpp | 12 +-
src/nanopolish_raw_loader.cpp | 282 +++++++
src/nanopolish_raw_loader.h | 25 +
src/nanopolish_read_db.cpp | 231 ++++++
src/nanopolish_read_db.h | 81 ++
src/nanopolish_scorereads.cpp | 8 +-
src/nanopolish_squiggle_read.cpp | 218 +++++-
src/nanopolish_squiggle_read.h | 34 +-
src/nanopolish_train_poremodel_from_basecalls.cpp | 4 +-
src/test/nanopolish_test.cpp | 19 +-
src/thirdparty/scrappie/event_detection.c | 319 ++++++++
src/thirdparty/scrappie/event_detection.h | 27 +
src/thirdparty/scrappie/scrappie_common.c | 73 ++
src/thirdparty/scrappie/scrappie_common.h | 10 +
src/thirdparty/scrappie/scrappie_stdlib.h | 41 +
src/thirdparty/scrappie/scrappie_structures.h | 31 +
src/thirdparty/scrappie/sse_mathfun.h | 722 +++++++++++++++++
src/thirdparty/scrappie/util.c | 288 +++++++
src/thirdparty/scrappie/util.h | 197 +++++
45 files changed, 2953 insertions(+), 1479 deletions(-)
diff --git a/Makefile b/Makefile
index 560f279..b735f2e 100644
--- a/Makefile
+++ b/Makefile
@@ -1,7 +1,7 @@
#
# Sub directories containing source code, except for the main programs
-SUBDIRS := src src/hmm src/thirdparty src/common src/alignment
+SUBDIRS := src src/hmm src/thirdparty src/thirdparty/scrappie src/common src/alignment
#
# Set libraries, paths, flags and options
@@ -11,7 +11,7 @@ SUBDIRS := src src/hmm src/thirdparty src/common src/alignment
LIBS=-lz
CXXFLAGS ?= -g -O3
CXXFLAGS += -std=c++11 -fopenmp
-CFLAGS ?= -O3
+CFLAGS ?= -O3 -std=c99
CXX ?= g++
CC ?= gcc
@@ -116,7 +116,7 @@ include .depend
$(CXX) -o $@ -c $(CXXFLAGS) $(CPPFLAGS) -fPIC $<
.c.o:
- $(CC) -o $@ -c $(CFLAGS) -fPIC $<
+ $(CC) -o $@ -c $(CFLAGS) $(H5_INCLUDE) -fPIC $<
# Link main executable
$(PROGRAM): src/main/nanopolish.o $(CPP_OBJ) $(C_OBJ) $(HTS_LIB) $(H5_LIB) $(EIGEN_CHECK)
diff --git a/README.md b/README.md
index 7466e42..ac308a2 100644
--- a/README.md
+++ b/README.md
@@ -53,16 +53,29 @@ nanopolish eventalign: align signal-level events to k-mers of a reference genome
## Analysis workflow examples
-### Computing a new consensus sequence for a draft assembly
+### Data preprocessing
-The original purpose of nanopolish was to compute an improved consensus sequence for a draft genome assembly produced by a long-read assembly like [canu](https://github.com/marbl/canu). This section describes how to do this, starting with your draft assembly which should have megabase-sized contigs.
+Nanopolish needs access to the signal-level data measured by the nanopore sequencer. The first step of any nanopolish workflow is to prepare the input data by telling nanopolish where to find the signal files. If you ran Albacore 2.0 on your data you should run `nanopolish index` on your input reads:
-First we prepare the data by extracting the reads from the FAST5 files, and aligning the basecalls to our draft assembly (`draft.fa`).
+```
+# Only run this if you used Albacore 2.0 or later
+nanopolish index -d /path/to/raw_fast5s/ albacore_output.fastq
+```
+If you basecalled your reads with Albacore 1.2 or earlier, you should run `nanopolish extract` on your input reads instead:
+
+```
+# Only run this if you used Albacore 1.2 or later
+nanopolish extract --type template directory/pass/ > reads.fa
```
-# Extract the QC-passed reads from a directory of FAST5 files
-nanopolish extract --type [2d|template] directory/pass/ > reads.fa
+Note these two commands are mutually exclusive - you only need to run one of them. You need to decide what command to run depending on the version of Albacore that you used. In the following sections we assume you have preprocessed the data by following the instructions above and that your reads are in a file named `reads.fa`.
+
+### Computing a new consensus sequence for a draft assembly
+
+The original purpose of nanopolish was to compute an improved consensus sequence for a draft genome assembly produced by a long-read assembly like [canu](https://github.com/marbl/canu). This section describes how to do this, starting with your draft assembly which should have megabase-sized contigs.
+
+```
# Index the draft genome
bwa index draft.fa
@@ -91,9 +104,6 @@ python nanopolish_merge.py polished.*.fa > polished_genome.fa
nanopolish can use the signal-level information measured by the sequencer to detect 5-mC as described [here](http://www.nature.com/nmeth/journal/vaop/ncurrent/full/nmeth.4184.html). Here's how you run it:
```
-# Extract all reads from a directory of FAST5 files
-nanopolish extract -r --type template directory/ > reads.fa
-
# Align the basecalled reads to a reference genome
bwa mem -x ont2d -t 8 reference.fa reads.fa | samtools sort -o reads.sorted.bam -T reads.tmp -
samtools index reads.sorted.bam
@@ -124,4 +134,4 @@ docker run -v /path/to/local/data/data/:/data/ -it :image_id ./nanopolish event
## Credits and Thanks
-The fast table-driven logsum implementation was provided by Sean Eddy as public domain code. This code was originally part of [hmmer3](http://hmmer.janelia.org/).
+The fast table-driven logsum implementation was provided by Sean Eddy as public domain code. This code was originally part of [hmmer3](http://hmmer.janelia.org/). Nanopolish also includes code from Oxford Nanopore's [scrappie](https://github.com/nanoporetech/scrappie) basecaller. This code is licensed under the MPL.
diff --git a/scripts/nanopolish_makerange.py b/scripts/nanopolish_makerange.py
index ee4b477..06418dd 100644
--- a/scripts/nanopolish_makerange.py
+++ b/scripts/nanopolish_makerange.py
@@ -14,12 +14,19 @@ recs = [ (rec.name, len(rec.seq)) for rec in SeqIO.parse(open(filename), "fasta"
SEGMENT_LENGTH = args.segment_length
OVERLAP_LENGTH = args.overlap_length
+MIN_SEGMENT_LENGTH = 5 * OVERLAP_LENGTH
for name, length in recs:
n_segments = (length / SEGMENT_LENGTH) + 1
- for n in range(0, length, SEGMENT_LENGTH):
- if ( n + SEGMENT_LENGTH) > length:
- print("%s:%d-%d" % (name, n, length - 1))
+ start = 0
+ while start < length:
+ end = start + SEGMENT_LENGTH
+
+ # If this segment will end near the end of the contig, extend it to end
+ if length - end < MIN_SEGMENT_LENGTH:
+ print("%s:%d-%d" % (name, start, length - 1))
+ start = length
else:
- print("%s:%d-%d" % (name, n, n + SEGMENT_LENGTH + OVERLAP_LENGTH))
+ print("%s:%d-%d" % (name, start, end + OVERLAP_LENGTH))
+ start = end
diff --git a/src/alignment/nanopolish_alignment_db.cpp b/src/alignment/nanopolish_alignment_db.cpp
index 85cc0d3..94e2766 100644
--- a/src/alignment/nanopolish_alignment_db.cpp
+++ b/src/alignment/nanopolish_alignment_db.cpp
@@ -85,9 +85,9 @@ AlignmentDB::AlignmentDB(const std::string& reads_file,
m_reference_file(reference_file),
m_sequence_bam(sequence_bam),
m_event_bam(event_bam),
- m_fast5_name_map(reads_file),
m_calibrate_on_load(calibrate_reads)
{
+ m_read_db.load(reads_file);
_clear_region();
}
@@ -250,6 +250,8 @@ std::vector<Variant> AlignmentDB::get_variants_in_region(const std::string& cont
{
std::vector<Variant> variants;
std::map<std::string, std::pair<Variant, int>> map;
+ assert(stop_position >= start_position);
+
std::vector<int> depth(stop_position - start_position + 1, 0);
for(size_t i = 0; i < m_sequence_records.size(); ++i) {
@@ -357,7 +359,13 @@ void AlignmentDB::load_region(const std::string& contig,
// Adjust end position to make sure we don't go out-of-range
m_region_contig = contig;
m_region_start = start_position;
- m_region_end = std::min(stop_position, faidx_seq_len(fai, contig.c_str()));
+ int contig_length = faidx_seq_len(fai, contig.c_str());
+ if(contig_length == -1) {
+ fprintf(stderr, "Error: could not retrieve length of contig %s from the faidx.\n", contig.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ m_region_end = std::min(stop_position, contig_length);
assert(!m_region_contig.empty());
assert(m_region_start >= 0);
@@ -601,8 +609,7 @@ void AlignmentDB::_load_squiggle_read(const std::string& read_name)
{
// Do we need to load this fast5 file?
if(m_squiggle_read_map.find(read_name) == m_squiggle_read_map.end()) {
- std::string fast5_path = m_fast5_name_map.get_path(read_name);
- SquiggleRead* sr = new SquiggleRead(read_name, fast5_path);
+ SquiggleRead* sr = new SquiggleRead(read_name, m_read_db);
m_squiggle_read_map[read_name] = sr;
}
}
diff --git a/src/alignment/nanopolish_alignment_db.h b/src/alignment/nanopolish_alignment_db.h
index 42b5f9a..0354b3a 100644
--- a/src/alignment/nanopolish_alignment_db.h
+++ b/src/alignment/nanopolish_alignment_db.h
@@ -147,7 +147,7 @@ class AlignmentDB
int m_region_end;
// cached alignments for a region
- Fast5Map m_fast5_name_map;
+ ReadDB m_read_db;
std::vector<SequenceAlignmentRecord> m_sequence_records;
std::vector<EventAlignmentRecord> m_event_records;
SquiggleReadMap m_squiggle_read_map;
diff --git a/src/alignment/nanopolish_anchor.cpp b/src/alignment/nanopolish_anchor.cpp
index 04dab90..7aa8f04 100644
--- a/src/alignment/nanopolish_anchor.cpp
+++ b/src/alignment/nanopolish_anchor.cpp
@@ -17,229 +17,6 @@
#include "nanopolish_methyltrain.h"
#include "nanopolish_squiggle_read.h"
-HMMRealignmentInput build_input_for_region(const std::string& bam_filename,
- const std::string& ref_filename,
- const Fast5Map& read_name_map,
- const std::string& contig_name,
- int start,
- int end,
- int stride)
-{
- // Initialize return data
- HMMRealignmentInput ret;
-
- // load bam file
- htsFile* bam_fh = sam_open(bam_filename.c_str(), "r");
- assert(bam_fh != NULL);
-
- // load bam index file
- std::string index_filename = bam_filename + ".bai";
- hts_idx_t* bam_idx = bam_index_load(index_filename.c_str());
- assert(bam_idx != NULL);
-
- // read the bam header
- bam_hdr_t* hdr = sam_hdr_read(bam_fh);
- int contig_id = bam_name2id(hdr, contig_name.c_str());
-
- // load reference fai file
- faidx_t *fai = fai_load(ref_filename.c_str());
-
- // Adjust end position to make sure we don't go out-of-range
- end = std::min(end, faidx_seq_len(fai, contig_name.c_str()));
-
- // load the reference sequence for this region
- int fetched_len = 0;
- char* ref_segment = faidx_fetch_seq(fai, contig_name.c_str(), start, end, &fetched_len);
- if(ref_segment == NULL) {
- exit(EXIT_FAILURE);
- }
- ret.original_sequence = ref_segment;
-
- // Initialize iteration
- bam1_t* record = bam_init1();
- hts_itr_t* itr = sam_itr_queryi(bam_idx, contig_id, start, end);
-
- // Iterate over reads aligned here
- std::vector<HMMReadAnchorSet> read_anchors;
- std::vector<std::vector<std::string>> read_substrings;
-
- // kmer size of pore model
- uint32_t k = 0;
-
- // Load the SquiggleReads aligned to this region and the bases
- // that are mapped to our reference anchoring positions
- int result;
- while((result = sam_itr_next(bam_fh, itr, record)) >= 0) {
-
- // Load a squiggle read for the mapped read
- std::string read_name = bam_get_qname(record);
- std::string fast5_path = read_name_map.get_path(read_name);
-
- // load read
- ret.reads.push_back(std::unique_ptr<SquiggleRead>(new SquiggleRead(read_name, fast5_path)));
- SquiggleRead& sr = *ret.reads.back();
-
- // Recalibrate each strand
- for(size_t strand_idx = 0; strand_idx < NUM_STRANDS; strand_idx++) {
- if(!sr.has_events_for_strand(strand_idx)) {
- continue;
- }
-
- std::vector<EventAlignment> ao = alignment_from_read(sr, strand_idx, -1,
- "", fai, hdr,
- record, -1, -1);
- recalibrate_model(sr, strand_idx, ao, &gDNAAlphabet, true, true);
- }
-
- k = sr.pore_model[T_IDX].k;
-
- // parse alignments to reference
- std::vector<AlignedPair> aligned_pairs = get_aligned_pairs(record);
- std::vector<int> read_bases_for_anchors =
- uniformally_sample_read_positions(aligned_pairs, start, end, stride);
-
- /*
- for(size_t i = 0; i < read_bases_for_anchors.size(); ++i) {
- printf("%d ", read_bases_for_anchors[i]);
- }
- printf("\n");
- */
-
- // Convert the read base positions into event indices for both strands
- HMMReadAnchorSet event_anchors;
- event_anchors.strand_anchors[T_IDX].resize(read_bases_for_anchors.size());
- event_anchors.strand_anchors[C_IDX].resize(read_bases_for_anchors.size());
-
- bool do_base_rc = bam_is_rev(record);
- bool template_rc = do_base_rc;
- bool complement_rc = !do_base_rc;
-
- for(size_t ai = 0; ai < read_bases_for_anchors.size(); ++ai) {
-
- int read_kidx = read_bases_for_anchors[ai];
-
- // read not aligned to this reference position
- if(read_kidx == -1) {
- continue;
- }
-
- if(do_base_rc)
- read_kidx = sr.flip_k_strand(read_kidx);
-
- // If the aligned base is beyong the start of the last k-mer of the read, skip
- if(read_kidx < 0 || read_kidx >= (int)(sr.read_sequence.size() - k + 1)) {
- continue;
- }
-
- int template_idx = sr.get_closest_event_to(read_kidx, T_IDX);
- int complement_idx = sr.has_events_for_strand(C_IDX) ? sr.get_closest_event_to(read_kidx, C_IDX) : -1;
-
- assert(template_idx != -1 && (!sr.has_events_for_strand(C_IDX) || complement_idx != -1));
- assert(template_idx < (int)sr.events[T_IDX].size());
- assert(sr.events[C_IDX].empty() || complement_idx < (int)sr.events[C_IDX].size());
-
- event_anchors.strand_anchors[T_IDX][ai] = { template_idx, template_rc };
- event_anchors.strand_anchors[C_IDX][ai] = { complement_idx, complement_rc };
-
- // If this is not the last anchor, extract the sequence of the read
- // from this anchor to the next anchor as an alternative assembly
- if(ai < read_bases_for_anchors.size() - 1) {
- int start_kidx = read_bases_for_anchors[ai];
- int end_kidx = read_bases_for_anchors[ai + 1];
- int max_kidx = sr.read_sequence.size() - k;
-
- // flip
- if(do_base_rc) {
- start_kidx = sr.flip_k_strand(start_kidx);
- end_kidx = sr.flip_k_strand(end_kidx);
-
- // swap
- int tmp = end_kidx;
- end_kidx = start_kidx;
- start_kidx = tmp;
- }
-
- // clamp values within range
- start_kidx = start_kidx >= 0 ? start_kidx : 0;
- end_kidx = end_kidx <= max_kidx ? end_kidx : max_kidx;
-
- std::string s = sr.read_sequence.substr(start_kidx, end_kidx - start_kidx + k);
-
- if(do_base_rc) {
- s = gDNAAlphabet.reverse_complement(s);
- }
-
- if(ai >= read_substrings.size())
- read_substrings.resize(ai + 1);
-
- read_substrings[ai].push_back(s);
- }
- }
-
- read_anchors.push_back(event_anchors);
- }
-
- // If there are no reads aligned to this segment return empty data
- if(read_anchors.empty()) {
- return ret;
- }
-
- // The HMMReadAnchorSet contains anchors for each strand of a read
- // laid out in a vector. Transpose this data so we have one anchor
- // for every read column-wise.
- assert(read_anchors.front().strand_anchors[T_IDX].size() ==
- read_anchors.back().strand_anchors[T_IDX].size());
-
- // same for all anchors, so we use the first
- size_t num_anchors = read_anchors.front().strand_anchors[T_IDX].size();
- size_t num_strands = read_anchors.size() * 2;
-
- ret.anchored_columns.resize(num_anchors);
- for(size_t ai = 0; ai < num_anchors; ++ai) {
-
- HMMAnchoredColumn& column = ret.anchored_columns[ai];
-
- for(size_t rai = 0; rai < read_anchors.size(); ++rai) {
- HMMReadAnchorSet& ras = read_anchors[rai];
- assert(ai < ras.strand_anchors[T_IDX].size());
- assert(ai < ras.strand_anchors[C_IDX].size());
-
- column.anchors.push_back(ras.strand_anchors[T_IDX][ai]);
- column.anchors.push_back(ras.strand_anchors[C_IDX][ai]);
- }
- assert(column.anchors.size() == num_strands);
-
- // Add sequences except for last anchor
- if(ai != num_anchors - 1) {
-
- // base, these sequences need to overlap by k - 1 bases
- int base_length = stride + k;
- if((int)ai * stride + base_length > fetched_len)
- base_length = fetched_len - ai * stride;
-
- column.base_sequence = std::string(ref_segment + ai * stride, base_length);
- column.base_contig = contig_name;
- column.base_start_position = start + ai * stride;
- assert(column.base_sequence.back() != '\0');
-
- // alts
- if(ai < read_substrings.size())
- column.alt_sequences = read_substrings[ai];
- }
- }
-
- // cleanup
- sam_itr_destroy(itr);
- bam_hdr_destroy(hdr);
- bam_destroy1(record);
- fai_destroy(fai);
- sam_close(bam_fh);
- hts_idx_destroy(bam_idx);
- free(ref_segment);
-
- return ret;
-}
-
std::vector<AlignedPair> get_aligned_pairs(const bam1_t* record, int read_stride)
{
std::vector<AlignedPair> out;
diff --git a/src/alignment/nanopolish_anchor.h b/src/alignment/nanopolish_anchor.h
index 2a3ff22..60666bc 100644
--- a/src/alignment/nanopolish_anchor.h
+++ b/src/alignment/nanopolish_anchor.h
@@ -14,7 +14,6 @@
#include "htslib/sam.h"
#include "nanopolish_common.h"
#include "nanopolish_squiggle_read.h"
-#include "nanopolish_fast5_map.h"
struct AlignedPair
{
@@ -80,24 +79,10 @@ struct HMMRealignmentInput
std::string original_sequence;
};
-// functions
-HMMRealignmentInput build_input_for_region(const std::string& bam_filename,
- const std::string& ref_filename,
- const Fast5Map& read_name_map,
- const std::string& contig_name,
- int start,
- int end,
- int stride);
-
// Return a vector specifying pairs of bases that have been aligned to each other
// This function can handle an "event cigar" bam record, which requires the ability
// for event indices to be in ascending or descending order. In the latter case
// read_stride should be -1
std::vector<AlignedPair> get_aligned_pairs(const bam1_t* record, int read_stride = 1);
-std::vector<int> uniformally_sample_read_positions(const std::vector<AlignedPair>& aligned_pairs,
- int ref_start,
- int ref_end,
- int ref_stride);
-
#endif
diff --git a/src/alignment/nanopolish_eventalign.cpp b/src/alignment/nanopolish_eventalign.cpp
index bc135ab..3ade148 100644
--- a/src/alignment/nanopolish_eventalign.cpp
+++ b/src/alignment/nanopolish_eventalign.cpp
@@ -28,7 +28,7 @@
#include "nanopolish_matrix.h"
#include "nanopolish_profile_hmm.h"
#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
+#include "nanopolish_read_db.h"
#include "nanopolish_hmm_input_sequence.h"
#include "nanopolish_pore_model_set.h"
#include "H5pubconf.h"
@@ -500,7 +500,7 @@ EventalignSummary summarize_alignment(const SquiggleRead& sr,
// Realign the read in event space
void realign_read(EventalignWriter writer,
- const Fast5Map& name_map,
+ const ReadDB& read_db,
const faidx_t* fai,
const bam_hdr_t* hdr,
const bam1_t* record,
@@ -510,10 +510,9 @@ void realign_read(EventalignWriter writer,
{
// Load a squiggle read for the mapped read
std::string read_name = bam_get_qname(record);
- std::string fast5_path = name_map.get_path(read_name);
// load read
- SquiggleRead sr(read_name, fast5_path, opt::write_samples ? SRF_LOAD_RAW_SAMPLES : 0);
+ SquiggleRead sr(read_name, read_db, opt::write_samples ? SRF_LOAD_RAW_SAMPLES : 0);
if(opt::verbose > 1) {
fprintf(stderr, "Realigning %s [%zu %zu]\n",
@@ -845,7 +844,8 @@ int eventalign_main(int argc, char** argv)
parse_eventalign_options(argc, argv);
omp_set_num_threads(opt::num_threads);
- Fast5Map name_map(opt::reads_file);
+ ReadDB read_db;
+ read_db.load(opt::reads_file);
// Open the BAM and iterate over reads
@@ -933,7 +933,7 @@ int eventalign_main(int argc, char** argv)
bam1_t* record = records[i];
size_t read_idx = num_reads_realigned + i;
if( (record->core.flag & BAM_FUNMAP) == 0) {
- realign_read(writer, name_map, fai, hdr, record, read_idx, clip_start, clip_end);
+ realign_read(writer, read_db, fai, hdr, record, read_idx, clip_start, clip_end);
}
}
diff --git a/src/common/fs_support.hpp b/src/common/fs_support.cpp
similarity index 67%
copy from src/common/fs_support.hpp
copy to src/common/fs_support.cpp
index 8ea2b0f..9b903ac 100644
--- a/src/common/fs_support.hpp
+++ b/src/common/fs_support.cpp
@@ -3,27 +3,22 @@
// Written by Matei David (matei.david at oicr.on.ca)
//---------------------------------------------------------
//
-#ifndef __FS_SUPPORT_HPP
-#define __FS_SUPPORT_HPP
-
-#include <string>
-#include <vector>
-
-#include <sys/types.h>
+#include "fs_support.hpp"
+#include <sys/stat.h>
#include <dirent.h>
-// This should work in windows.
-// Ref:
-// http://stackoverflow.com/a/612176/717706
-
+//
bool is_directory(const std::string& file_name)
{
auto dir = opendir(file_name.c_str());
- if (not dir) return false;
+ if(not dir) {
+ return false;
+ }
closedir(dir);
return true;
}
+//
std::vector< std::string > list_directory(const std::string& file_name)
{
std::vector< std::string > res;
@@ -31,13 +26,12 @@ std::vector< std::string > list_directory(const std::string& file_name)
struct dirent *ent;
dir = opendir(file_name.c_str());
- if (not dir) return res;
- while ((ent = readdir(dir)) != nullptr)
- {
+ if(not dir) {
+ return res;
+ }
+ while((ent = readdir(dir)) != nullptr) {
res.push_back(ent->d_name);
}
closedir(dir);
return res;
}
-
-#endif
diff --git a/src/common/fs_support.hpp b/src/common/fs_support.hpp
index 8ea2b0f..5c5a7ad 100644
--- a/src/common/fs_support.hpp
+++ b/src/common/fs_support.hpp
@@ -9,35 +9,11 @@
#include <string>
#include <vector>
-#include <sys/types.h>
-#include <dirent.h>
+// ref: http://stackoverflow.com/a/612176/717706
+// return true if the given name is a directory
+bool is_directory(const std::string& file_name);
-// This should work in windows.
-// Ref:
-// http://stackoverflow.com/a/612176/717706
-
-bool is_directory(const std::string& file_name)
-{
- auto dir = opendir(file_name.c_str());
- if (not dir) return false;
- closedir(dir);
- return true;
-}
-
-std::vector< std::string > list_directory(const std::string& file_name)
-{
- std::vector< std::string > res;
- DIR* dir;
- struct dirent *ent;
-
- dir = opendir(file_name.c_str());
- if (not dir) return res;
- while ((ent = readdir(dir)) != nullptr)
- {
- res.push_back(ent->d_name);
- }
- closedir(dir);
- return res;
-}
+// return a vector of files within the given directory
+std::vector< std::string > list_directory(const std::string& file_name);
#endif
diff --git a/src/common/nanopolish_common.h b/src/common/nanopolish_common.h
index 6518a88..7739ab5 100644
--- a/src/common/nanopolish_common.h
+++ b/src/common/nanopolish_common.h
@@ -18,7 +18,7 @@
#include "logsum.h"
#define PACKAGE_NAME "nanopolish"
-#define PACKAGE_VERSION "0.7.2"
+#define PACKAGE_VERSION "0.8.1"
#define PACKAGE_BUGREPORT "https://github.com/jts/nanopolish/issues"
//
diff --git a/src/common/nanopolish_fast5_map.cpp b/src/common/nanopolish_fast5_map.cpp
deleted file mode 100644
index a53a8a8..0000000
--- a/src/common/nanopolish_fast5_map.cpp
+++ /dev/null
@@ -1,138 +0,0 @@
-//---------------------------------------------------------
-// Copyright 2015 Ontario Institute for Cancer Research
-// Written by Jared Simpson (jared.simpson at oicr.on.ca)
-//---------------------------------------------------------
-//
-// nanopolish_fast5_map - a simple map from a read
-// name to a fast5 file path
-#include <stdio.h>
-#include <inttypes.h>
-#include <zlib.h>
-#include <fstream>
-#include <ostream>
-#include <iostream>
-#include <sys/stat.h>
-#include "nanopolish_fast5_map.h"
-#include "nanopolish_common.h"
-#include "htslib/kseq.h"
-
-//
-#define FOFN_SUFFIX ".fast5.fofn"
-
-KSEQ_INIT(gzFile, gzread)
-
-Fast5Map::Fast5Map(const std::string& fasta_filename)
-{
- // If the fofn file exists, load from it
- // otherwise parse the entire fasta file
- std::string fofn_filename = fasta_filename + FOFN_SUFFIX;
- struct stat fofn_file_s;
- struct stat fasta_file_s;
- int fofn_ret = stat(fofn_filename.c_str(), &fofn_file_s);
- stat(fasta_filename.c_str(), &fasta_file_s);
-
- // Use the stored fofn if its available and newer than the fasta
- if(fofn_ret == 0 && fofn_file_s.st_mtime > fasta_file_s.st_mtime) {
- load_from_fofn(fofn_filename);
- } else {
- load_from_fasta(fasta_filename);
- }
-}
-
-std::string Fast5Map::get_path(const std::string& read_name) const
-{
- std::map<std::string, std::string>::const_iterator
- iter = read_to_path_map.find(read_name);
-
- if(iter == read_to_path_map.end()) {
- fprintf(stderr, "error: could not find fast5 path for %s\n", read_name.c_str());
- exit(EXIT_FAILURE);
- }
-
- return iter->second;
-}
-
-//
-void Fast5Map::load_from_fasta(std::string fasta_filename)
-{
- gzFile gz_fp;
-
- FILE* fp = fopen(fasta_filename.c_str(), "r");
- if(fp == NULL) {
- fprintf(stderr, "error: could not open %s for read\n", fasta_filename.c_str());
- exit(EXIT_FAILURE);
- }
-
- gz_fp = gzdopen(fileno(fp), "r");
- if(gz_fp == NULL) {
- fprintf(stderr, "error: could not open %s using gzdopen\n", fasta_filename.c_str());
- exit(EXIT_FAILURE);
- }
-
- kseq_t* seq = kseq_init(gz_fp);
-
- while(kseq_read(seq) >= 0) {
- if(seq->comment.l == 0) {
- fprintf(stderr, "error: no path associated with read %s\n", seq->name.s);
- exit(EXIT_FAILURE);
- }
-
- // This splitting code implicitly handles both the 2 and 3 field
- // fasta format that poretools will output. The FAST5 path
- // is always the last field.
- std::vector<std::string> fields = split(seq->comment.s, ' ');
- if(read_to_path_map.find(seq->name.s) != read_to_path_map.end()) {
- fprintf(stderr, "Error: duplicate read name %s found in fasta file\n", seq->name.s);
- exit(EXIT_FAILURE);
- }
- read_to_path_map[seq->name.s] = fields.back();
- }
-
- kseq_destroy(seq);
- gzclose(gz_fp);
- fclose(fp);
-
- // Sanity check that the first path actually points to a file
- if(read_to_path_map.size() > 0) {
- std::string first_read = read_to_path_map.begin()->first;
- std::string first_path = read_to_path_map.begin()->second;
- struct stat file_s;
- int ret = stat(first_path.c_str(), &file_s);
- if(ret != 0) {
- fprintf(stderr, "Error: could not find path to FAST5 for read %s\n", first_read.c_str());
- fprintf(stderr, "Please make sure that this path is accessible: %s\n", first_path.c_str());
- exit(EXIT_FAILURE);
- }
- }
-
- // Write the map as a fofn file so next time we don't have to parse
- // the entire fasta
- write_to_fofn(fasta_filename + FOFN_SUFFIX);
-}
-
-void Fast5Map::write_to_fofn(std::string fofn_filename)
-{
- std::ofstream outfile(fofn_filename.c_str());
-
- for(std::map<std::string, std::string>::iterator iter = read_to_path_map.begin();
- iter != read_to_path_map.end(); ++iter) {
- outfile << iter->first << "\t" << iter->second << "\n";
- }
-}
-
-//
-void Fast5Map::load_from_fofn(std::string fofn_filename)
-{
- std::ifstream infile(fofn_filename.c_str());
-
- if(infile.bad()) {
- fprintf(stderr, "error: could not read fofn %s\n", fofn_filename.c_str());
- exit(EXIT_FAILURE);
- }
-
- std::string name;
- std::string path;
- while(infile >> name >> path) {
- read_to_path_map[name] = path;
- }
-}
diff --git a/src/common/nanopolish_fast5_map.h b/src/common/nanopolish_fast5_map.h
deleted file mode 100644
index 614829a..0000000
--- a/src/common/nanopolish_fast5_map.h
+++ /dev/null
@@ -1,39 +0,0 @@
-//---------------------------------------------------------
-// Copyright 2015 Ontario Institute for Cancer Research
-// Written by Jared Simpson (jared.simpson at oicr.on.ca)
-//---------------------------------------------------------
-//
-// nanopolish_fast5_map - a simple map from a read
-// name to a fast5 file path
-#ifndef NANOPOLISH_FAST5_MAP
-#define NANOPOLISH_FAST5_MAP
-
-#include <string>
-#include <map>
-
-class Fast5Map
-{
- public:
- Fast5Map(const std::string& fasta_filename);
-
- // return the path for the given read name
- // if the read does not exist in the map, emits an error
- // and exits
- std::string get_path(const std::string& read_name) const;
-
- private:
-
- // Read the read -> path map from the header of a fasta file
- void load_from_fasta(std::string fasta_filename);
-
- // Read the map from a pre-computed .fofn file
- void load_from_fofn(std::string fofn_filename);
-
- // Write the map to the .fofn file
- void write_to_fofn(std::string fofn_filename);
-
- std::map<std::string, std::string> read_to_path_map;
-
-};
-
-#endif
diff --git a/src/common/nanopolish_variant.cpp b/src/common/nanopolish_variant.cpp
index e412e24..0a441de 100644
--- a/src/common/nanopolish_variant.cpp
+++ b/src/common/nanopolish_variant.cpp
@@ -714,19 +714,22 @@ Variant score_variant_thresholded(const Variant& input_variant,
void annotate_variants_with_all_support(std::vector<Variant>& input, const AlignmentDB& alignments, int min_flanking_sequence, const uint32_t alignment_flags)
{
+ Haplotype ref_haplotype(alignments.get_region_contig(), alignments.get_region_start(), alignments.get_reference());
+
for(size_t vi = 0; vi < input.size(); vi++) {
Variant& v = input[vi];
std::string ref_name = v.ref_name;
int calling_start = v.ref_position - min_flanking_sequence;
int calling_end = v.ref_position + min_flanking_sequence;
-
+
// Construct haplotype set for A,C,G,T at this position
std::vector<Haplotype> haplotypes;
- Haplotype ref_haplotype(alignments.get_region_contig(), alignments.get_region_start(), alignments.get_reference());
+ Haplotype calling_haplotype = ref_haplotype.substr_by_reference(calling_start, calling_end);
+
Variant tmp_variant = v;
for(size_t bi = 0; bi < 4; bi++) {
- Haplotype variant_haplotype = ref_haplotype;
+ Haplotype variant_haplotype = calling_haplotype;
tmp_variant.alt_seq = "ACGT"[bi];
variant_haplotype.apply_variant(tmp_variant);
haplotypes.push_back(variant_haplotype);
@@ -746,7 +749,7 @@ void annotate_variants_with_all_support(std::vector<Variant>& input, const Align
scores.push_back(s);
if(hi > 0) {
sum_score = add_logs(s, sum_score);
- } else {
+ } else {
sum_score = s;
}
}
@@ -755,7 +758,7 @@ void annotate_variants_with_all_support(std::vector<Variant>& input, const Align
support_fraction[hi] += exp(scores[hi] - sum_score);
}
}
-
+
std::stringstream ss;
ss << std::fixed << std::setprecision(3);
for(size_t hi = 0; hi < haplotypes.size(); ++hi) {
diff --git a/src/main/nanopolish.cpp b/src/main/nanopolish.cpp
index 5746390..097edc6 100644
--- a/src/main/nanopolish.cpp
+++ b/src/main/nanopolish.cpp
@@ -9,9 +9,9 @@
#include <map>
#include <functional>
#include "logsum.h"
+#include "nanopolish_index.h"
#include "nanopolish_extract.h"
#include "nanopolish_call_variants.h"
-#include "nanopolish_consensus.h"
#include "nanopolish_eventalign.h"
#include "nanopolish_getmodel.h"
#include "nanopolish_methyltrain.h"
@@ -27,16 +27,15 @@ static std::map< std::string, std::function<int(int, char**)> > programs = {
{"help", print_usage},
{"--help", print_usage},
{"--version", print_version},
+ {"index", index_main},
{"extract", extract_main},
- {"consensus", consensus_main},
{"eventalign", eventalign_main},
{"getmodel", getmodel_main},
{"variants", call_variants_main},
{"methyltrain", methyltrain_main},
{"scorereads", scorereads_main} ,
{"phase-reads", phase_reads_main} ,
- {"call-methylation", call_methylation_main},
- {"train-poremodel-from-basecalls", train_poremodel_from_basecalls_main}
+ {"call-methylation", call_methylation_main}
};
int print_usage(int, char **)
diff --git a/src/nanopolish_call_methylation.cpp b/src/nanopolish_call_methylation.cpp
index 62eb5ef..5a20aa3 100644
--- a/src/nanopolish_call_methylation.cpp
+++ b/src/nanopolish_call_methylation.cpp
@@ -28,11 +28,11 @@
#include "nanopolish_matrix.h"
#include "nanopolish_profile_hmm.h"
#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
#include "nanopolish_methyltrain.h"
#include "nanopolish_pore_model_set.h"
#include "nanopolish_bam_processor.h"
#include "nanopolish_alignment_db.h"
+#include "nanopolish_read_db.h"
#include "H5pubconf.h"
#include "profiler.h"
#include "progress.h"
@@ -136,7 +136,7 @@ static const struct option longopts[] = {
// Test CpG sites in this read for methylation
void calculate_methylation_for_read(const OutputHandles& handles,
- const Fast5Map& name_map,
+ const ReadDB& read_db,
const faidx_t* fai,
const bam_hdr_t* hdr,
const bam1_t* record,
@@ -146,8 +146,7 @@ void calculate_methylation_for_read(const OutputHandles& handles,
{
// Load a squiggle read for the mapped read
std::string read_name = bam_get_qname(record);
- std::string fast5_path = name_map.get_path(read_name);
- SquiggleRead sr(read_name, fast5_path);
+ SquiggleRead sr(read_name, read_db);
// An output map from reference positions to scored CpG sites
std::map<int, ScoredSite> site_score_map;
@@ -401,7 +400,8 @@ void parse_call_methylation_options(int argc, char** argv)
int call_methylation_main(int argc, char** argv)
{
parse_call_methylation_options(argc, argv);
- Fast5Map name_map(opt::reads_file);
+ ReadDB read_db;
+ read_db.load(opt::reads_file);
// load reference fai file
faidx_t *fai = fai_load(opt::genome_file.c_str());
@@ -426,7 +426,7 @@ int call_methylation_main(int argc, char** argv)
// the BamProcessor framework calls the input function with the
// bam record, read index, etc passed as parameters
// bind the other parameters the worker function needs here
- auto f = std::bind(calculate_methylation_for_read, std::ref(handles), name_map, fai, _1, _2, _3, _4, _5);
+ auto f = std::bind(calculate_methylation_for_read, std::ref(handles), std::ref(read_db), fai, _1, _2, _3, _4, _5);
BamProcessor processor(opt::bam_file, opt::region, opt::num_threads);
processor.parallel_run(f);
diff --git a/src/nanopolish_call_variants.cpp b/src/nanopolish_call_variants.cpp
index ce005a9..9b59280 100644
--- a/src/nanopolish_call_variants.cpp
+++ b/src/nanopolish_call_variants.cpp
@@ -30,7 +30,6 @@
#include "nanopolish_profile_hmm.h"
#include "nanopolish_alignment_db.h"
#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
#include "nanopolish_variant.h"
#include "nanopolish_haplotype.h"
#include "nanopolish_pore_model_set.h"
@@ -1060,17 +1059,27 @@ int call_variants_main(int argc, char** argv)
std::string contig;
int start_base;
int end_base;
+ int contig_length = -1;
// If a window has been specified, only call variants/polish in that range
if(!opt::window.empty()) {
// Parse the window string
parse_region_string(opt::window, contig, start_base, end_base);
- end_base = std::min(end_base, get_contig_length(contig) - 1);
+ contig_length = get_contig_length(contig);
+ end_base = std::min(end_base, contig_length - 1);
} else {
// otherwise, run on the whole genome
contig = get_single_contig_or_fail();
+ contig_length = get_contig_length(contig);
start_base = 0;
- end_base = get_contig_length(contig) - 1;
+ end_base = contig_length - 1;
+ }
+
+ int MIN_DISTANCE_TO_END = 40;
+ if(contig_length - start_base < MIN_DISTANCE_TO_END) {
+ fprintf(stderr, "Invalid polishing window: [%d %d] - please adjust -w parameter.\n", start_base, end_base);
+ fprintf(stderr, "The starting coordinate of the polishing window must be at least %dbp from the contig end\n", MIN_DISTANCE_TO_END);
+ exit(EXIT_FAILURE);
}
FILE* out_fp;
diff --git a/src/nanopolish_consensus.cpp b/src/nanopolish_consensus.cpp
deleted file mode 100644
index f565ebe..0000000
--- a/src/nanopolish_consensus.cpp
+++ /dev/null
@@ -1,897 +0,0 @@
-//---------------------------------------------------------
-// Copyright 2015 Ontario Institute for Cancer Research
-// Written by Jared Simpson (jared.simpson at oicr.on.ca)
-//---------------------------------------------------------
-//
-// nanopolish_consensus.cpp -- entry point to consensus functions
-//
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <string>
-#include <vector>
-#include <inttypes.h>
-#include <assert.h>
-#include <math.h>
-#include <sys/time.h>
-#include <algorithm>
-#include <sstream>
-#include <set>
-#include <omp.h>
-#include <getopt.h>
-#include "nanopolish_poremodel.h"
-#include "nanopolish_transition_parameters.h"
-#include "nanopolish_matrix.h"
-#include "nanopolish_klcs.h"
-#include "nanopolish_profile_hmm.h"
-#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
-#include "nanopolish_hmm_input_sequence.h"
-#include "nanopolish_pore_model_set.h"
-#include "profiler.h"
-#include "progress.h"
-#include "stdaln.h"
-
-// Macros
-#define max3(x,y,z) std::max(std::max(x,y), z)
-
-// Flags to turn on/off debugging information
-
-//#define DEBUG_HMM_UPDATE 1
-//#define DEBUG_HMM_EMISSION 1
-//#define DEBUG_TRANSITION 1
-//#define DEBUG_PATH_SELECTION 1
-//#define DEBUG_SINGLE_SEGMENT 1
-//#define DEBUG_SHOW_TOP_TWO 1
-//#define DEBUG_SEGMENT_ID 193
-//#define DEBUG_BENCHMARK 1
-
-//
-// Getopt
-//
-#define SUBPROGRAM "consensus"
-
-static const char *CONSENSUS_VERSION_MESSAGE =
-SUBPROGRAM " Version " PACKAGE_VERSION "\n"
-"Written by Jared Simpson.\n"
-"\n"
-"Copyright 2015 Ontario Institute for Cancer Research\n";
-
-static const char *CONSENSUS_USAGE_MESSAGE =
-"Usage: " PACKAGE_NAME " " SUBPROGRAM " [OPTIONS] -w contig:start-end --reads reads.fa --bam alignments.bam --genome genome.fa\n"
-"Compute a new consensus sequence for an assembly using a signal-level HMM\n"
-"\n"
-" -v, --verbose display verbose output\n"
-" --version display version\n"
-" --help display this help and exit\n"
-" -w, --window=STR compute the consensus for window STR (format: ctg:start_id-end_id)\n"
-" -r, --reads=FILE the 2D ONT reads are in fasta FILE\n"
-" -b, --bam=FILE the reads aligned to the genome assembly are in bam FILE\n"
-" -g, --genome=FILE the genome we are computing a consensus for is in FILE\n"
-" -o, --outfile=FILE write result to FILE [default: stdout]\n"
-" -t, --threads=NUM use NUM threads (default: 1)\n"
-" --models-fofn=FILE read alternative k-mer models from FILE\n"
-"\nReport bugs to " PACKAGE_BUGREPORT "\n\n";
-
-namespace opt
-{
- static unsigned int verbose;
- static std::string reads_file;
- static std::string bam_file;
- static std::string genome_file;
- static std::string output_file;
- static std::string window;
- static std::string models_fofn;
- static int show_progress = 0;
- static int num_threads = 1;
-}
-
-static const char* shortopts = "r:b:g:t:w:o:v";
-
-enum { OPT_HELP = 1, OPT_VERSION, OPT_VCF, OPT_PROGRESS, OPT_MODELS_FOFN };
-
-static const struct option longopts[] = {
- { "verbose", no_argument, NULL, 'v' },
- { "reads", required_argument, NULL, 'r' },
- { "bam", required_argument, NULL, 'b' },
- { "genome", required_argument, NULL, 'g' },
- { "window", required_argument, NULL, 'w' },
- { "outfile", required_argument, NULL, 'o' },
- { "threads", required_argument, NULL, 't' },
- { "models-fofn", required_argument, NULL, OPT_MODELS_FOFN },
- { "progress", no_argument, NULL, OPT_PROGRESS },
- { "help", no_argument, NULL, OPT_HELP },
- { "version", no_argument, NULL, OPT_VERSION },
- { NULL, 0, NULL, 0 }
-};
-
-std::vector<HMMInputData> get_input_for_columns(HMMRealignmentInput& window,
- const HMMAnchoredColumn& start_column,
- const HMMAnchoredColumn& end_column)
-{
- assert(start_column.anchors.size() == end_column.anchors.size());
-
- std::vector<HMMInputData> input;
- for(uint32_t rsi = 0; rsi < start_column.anchors.size(); ++rsi) {
-
- HMMStrandAnchor start_sa = start_column.anchors[rsi];
- HMMStrandAnchor end_sa = end_column.anchors[rsi];
-
- // sanity checks
- // This read strand does not have events at both anchors
- if(start_sa.event_idx == -1 || end_sa.event_idx == -1)
- continue;
-
- // require a minimum number of events
- uint32_t n_events = abs(start_sa.event_idx - end_sa.event_idx);
- if(n_events < 20 || n_events > 500)
- continue;
-
- if(start_sa.rc != end_sa.rc)
- continue;
-
-
- HMMInputData data;
-
- uint32_t read_idx = rsi / 2;
- assert(read_idx < window.reads.size());
- data.anchor_index = rsi;
- data.read = window.reads[read_idx].get();
- data.strand = rsi % 2;
- data.event_start_idx = start_sa.event_idx;
- data.event_stop_idx = end_sa.event_idx;
- if(data.event_start_idx < data.event_stop_idx)
- data.event_stride = 1;
- else
- data.event_stride = -1;
- data.rc = start_sa.rc;
-
- input.push_back(data);
- }
- return input;
-}
-
-// Handy wrappers for scoring/debugging functions
-// The consensus algorithms call into these so we can switch
-// scoring functions without writing a bunch of code
-double score_sequence(const std::string& sequence, const HMMInputData& data)
-{
- return profile_hmm_score(sequence, data);
-}
-
-void update_training_with_segment(const HMMInputSequence& sequence, const HMMInputData& data)
-{
- std::vector<HMMAlignmentState> alignment = profile_hmm_align(sequence, data);
- data.read->parameters[data.strand].add_training_from_alignment(sequence, data, alignment);
-}
-
-struct PathCons
-{
- // default constructor
- PathCons(const std::string& s) : path(s), score(0.0f), sum_rank(0) {}
-
- std::string path;
-
- double score;
- size_t sum_rank;
- size_t num_improved;
- size_t num_scored;
-
- std::string mutdesc;
-
-};
-typedef std::vector<PathCons> PathConsVector;
-
-bool sortPathConsScoreDesc(const PathCons& a, const PathCons& b)
-{
- return a.score > b.score;
-}
-
-bool sortPathConsRankAsc(const PathCons& a, const PathCons& b)
-{
- return a.sum_rank < b.sum_rank;
-}
-
-bool sortPathConsRankDesc(const PathCons& a, const PathCons& b)
-{
- return a.sum_rank > b.sum_rank;
-}
-
-struct IndexedPathScore
-{
- double score;
- uint32_t path_index;
-};
-
-bool sortIndexedPathScoreDesc(const IndexedPathScore& a, const IndexedPathScore& b)
-{
- return a.score > b.score;
-}
-
-// This scores each path using the HMM and
-// sorts the paths into ascending order by score
-void score_paths(PathConsVector& paths, const std::vector<HMMInputData>& input)
-{
- PROFILE_FUNC("score_paths")
- size_t CULL_RATE = 5;
- double CULL_MIN_SCORE = -30.0f;
- double CULL_MIN_IMPROVED_FRACTION = 0.2f;
-
- // cache the initial sequence
- std::string first = paths[0].path;
-
- PathConsVector dedup_paths;
-
- // initialize and deduplicate paths to avoid redundant computation
- std::set<std::string> path_string_set;
- for(size_t pi = 0; pi < paths.size(); ++pi) {
-
- if(path_string_set.find(paths[pi].path) == path_string_set.end()) {
- paths[pi].score = 0;
- paths[pi].sum_rank = 0;
- paths[pi].num_improved = 0;
- paths[pi].num_scored = 0;
- dedup_paths.push_back(paths[pi]);
- path_string_set.insert(paths[pi].path);
- }
- }
- paths.clear();
- paths.swap(dedup_paths);
-
-
- // Score all reads
- for(uint32_t ri = 0; ri < input.size(); ++ri) {
-
- if(opt::verbose > 2) {
- fprintf(stderr, "Scoring %d\n", ri);
- }
-
- //const HMMInputData& data = input[ri];
- std::vector<IndexedPathScore> result(paths.size());
-
- // Score all paths
- #pragma omp parallel for
- for(size_t pi = 0; pi < paths.size(); ++pi) {
- double curr = score_sequence(paths[pi].path, input[ri]);
- result[pi].score = curr;
- result[pi].path_index = pi;
- }
-
- // Save score of first path
- double first_path_score = result[0].score;
-
- // Sort result by score
- std::stable_sort(result.begin(), result.end(), sortIndexedPathScoreDesc);
-
- for(size_t pri = 0; pri < result.size(); ++pri) {
- size_t pi = result[pri].path_index;
-
- paths[pi].score += (result[pri].score - first_path_score);
- uint32_t rank_score = pri;
- paths[pi].sum_rank += rank_score;
- paths[pi].num_improved += (result[pri].score > first_path_score);
- paths[pi].num_scored += 1;
- }
-
- // Cull paths
- if(ri > 0 && ri % CULL_RATE == 0) {
- PathConsVector retained_paths;
- for(size_t pi = 0; pi < paths.size(); ++pi) {
-
- // We keep a path if any of these conditions are met:
- // 1) it is the original unmodified sequence
- // 2) its score is greater than CULL_MIN_SCORE
- // 3) the fraction of reads that score better on this
- // path compared to the original sequence is greater
- // than CULL_MIN_IMPROVED_FRACTION
- double f = (double)paths[pi].num_improved / (double)paths[pi].num_scored;
- if(pi == 0 || paths[pi].score > CULL_MIN_SCORE || f >= CULL_MIN_IMPROVED_FRACTION) {
- retained_paths.push_back(paths[pi]);
- }
- }
- paths.swap(retained_paths);
- }
- }
-
- // select new sequence
- //std::stable_sort(paths.begin(), paths.end(), sortPathConsRankAsc);
- std::stable_sort(paths.begin(), paths.end(), sortPathConsScoreDesc);
-
-#if DEBUG_PATH_SELECTION
- double MIN_FIT = INFINITY;
- for(size_t pi = 0; pi < paths.size(); ++pi) {
-
- // Calculate the length of the matching prefix with the initial sequence
- const std::string& s = paths[pi].path;
-
- char initial = s == first ? 'I' : ' ';
-
- fprintf(stderr, "%zu\t%s\t%.1lf\t%zu %c %s", pi, paths[pi].path.c_str(), paths[pi].score, paths[pi].sum_rank, initial, paths[pi].mutdesc.c_str());
- // If this is the truth path or the best path, show the scores for all reads
- if(pi <= 1 || initial == 'I') {
- for(uint32_t ri = 0; ri < input.size(); ++ri) {
- const HMMInputData& data = input[ri];
- double curr = score_sequence(paths[pi].path, input[ri]);
- fprintf(stderr, "%.2lf ", curr);
- }
- }
- fprintf(stderr, "\n");
- }
-#endif
-
-}
-
-void extend_paths(PathConsVector& paths, int maxk = 2)
-{
- // Insert all possible extensions into the path sequence
- // for k in 1 to maxk
- PathConsVector new_paths;
-
- for(int k = 1; k <= maxk; ++k) {
-
- for(unsigned pi = 0; pi < paths.size(); ++pi) {
-
- std::string first(k, 'A');
- std::string extension = first;
-
- do {
- std::string current = paths[pi].path;
- std::string ns = current.insert(current.size() - 5, extension);
- PathCons ps(ns);
- new_paths.push_back(ps);
- gDNAAlphabet.lexicographic_next(extension);
- } while(extension != first);
- }
- }
-
- paths.swap(new_paths);
-}
-
-PathConsVector generate_mutations(const std::string& sequence, const uint32_t k)
-{
- PathConsVector mutations;
-
- // Add the unmutated sequence
- {
- PathCons pc(sequence);
- mutations.push_back(pc);
- }
-
- // Mutate every base except for in the first/last k-mer
- for(size_t si = k; si < sequence.size() - k; ++si) {
-
- // All subs
- for(size_t bi = 0; bi < 4; bi++) {
- char b = "ACGT"[bi];
- if(sequence[si] == b)
- continue;
- PathCons pc(sequence);
- pc.path[si] = b;
- std::stringstream ss;
- ss << "sub-" << si << "-" << b;
- pc.mutdesc = ss.str();
- mutations.push_back(pc);
- }
-
- // 1bp del at this position
- {
- PathCons pc(sequence);
- pc.path.erase(si, 1);
-
- std::stringstream ss;
- ss << "del-" << si;
- pc.mutdesc = ss.str();
-
- mutations.push_back(pc);
- }
-
- // All 1bp ins before this position
- for(size_t bi = 0; bi < 4; bi++) {
- char b = "ACGT"[bi];
- PathCons pc(sequence);
- pc.path.insert(si, 1, b);
-
- std::stringstream ss;
- ss << "ins-" << si << "-" << b;
- pc.mutdesc = ss.str();
-
- mutations.push_back(pc);
- }
- }
-
- return mutations;
-}
-
-// Run the mutation algorithm to generate an improved consensus sequence
-std::string run_mutation(const std::string& base, const std::vector<HMMInputData>& input)
-{
- PROFILE_FUNC("run_mutation")
- std::string result = base;
-
- // assume models for all the reads have the same k
- assert(!input.empty());
- const uint32_t k = input[0].read->pore_model[input[0].strand].k;
-
- int iteration = 0;
- while(iteration++ < 10) {
-
- // Generate possible sequences
- PathConsVector paths = generate_mutations(result, k);
-
- // score them in the HMM
- score_paths(paths, input);
-
- // check if no improvement was made
- if(paths[0].path == result)
- break;
- result = paths[0].path;
- }
-
- return result;
-}
-
-void generate_alt_paths(PathConsVector& paths, const std::string& base, const std::vector<std::string>& alts,
- const uint32_t k)
-{
- // Generate alternatives
- for(uint32_t ai = 0; ai < alts.size(); ++ai) {
- const std::string& alt = alts[ai];
-
- if(alt.size() < k)
- continue;
-
- kLCSResult result = kLCS(base, alt, k);
-
-#ifdef DEBUG_ALT_GENERATION
- printf("Match to alt %s\n", alt.c_str());
- for(size_t mi = 0; mi < result.size(); ++mi) {
- std::string extend = "";
- if(mi < result.size() - 1 && result[mi].j + 1 != result[mi + 1].j) {
- extend = alt.substr(result[mi].j, result[mi + 1].j - result[mi].j + k);
- }
- printf("\t%zu %zu %s %s\n", result[mi].i, result[mi].j, base.substr(result[mi].i, k).c_str(), extend.c_str());
- }
-#endif
-
- uint32_t match_idx = 0;
- while(match_idx < result.size()) {
- uint32_t last_idx = result.size() - 1;
-
- // advance the match to the next point of divergence
- while(match_idx != last_idx &&
- result[match_idx].i == result[match_idx + 1].i - 1 &&
- result[match_idx].j == result[match_idx + 1].j - 1) {
- match_idx++;
- }
- // no more divergences to process
- if(match_idx == last_idx)
- break;
-
- uint32_t bl = result[match_idx + 1].i - result[match_idx].i;
- uint32_t rl = result[match_idx + 1].j - result[match_idx].j;
-
- std::string base_subseq = base.substr(result[match_idx].i, bl);
- std::string alt_subseq = alt.substr(result[match_idx].j, rl);
-
- // Perform the splice
- PathCons new_path(base);
- new_path.path.replace(result[match_idx].i, bl, alt_subseq);
- paths.push_back(new_path);
-
- match_idx += 1;
- }
- }
-}
-
-// Run the block substitution algorithm to generate an improved consensus sequence
-std::string run_block_substitution(const std::string& base,
- const std::vector<HMMInputData>& input,
- const std::vector<std::string>& alts)
-{
- std::string result = base;
-
- // assume models for all the reads have the same k
- assert(!input.empty());
- const uint32_t k = input[0].read->pore_model[input[0].strand].k;
-
- uint32_t max_rounds = 6;
- uint32_t round = 0;
- while(round++ < max_rounds) {
-
- PathConsVector paths;
- PathCons initial_path(result);
- paths.push_back(initial_path);
-
- generate_alt_paths(paths, result, alts, k);
- score_paths(paths, input);
-
- if(paths[0].path == result)
- break;
- result = paths[0].path;
- }
- return result;
-}
-
-//
-// Outlier filtering
-//
-void filter_outlier_data(std::vector<HMMInputData>& input, const std::string& sequence)
-{
- std::vector<HMMInputData> out_rs;
- for(uint32_t ri = 0; ri < input.size(); ++ri) {
- const HMMInputData& rs = input[ri];
-
- double curr = score_sequence(sequence, rs);
- double n_events = abs((int)rs.event_start_idx - (int)rs.event_stop_idx) + 1.0f;
- double lp_per_event = curr / n_events;
-
- if(opt::verbose >= 1) {
- fprintf(stderr, "OUTLIER_FILTER %d %.2lf %.2lf %.2lf\n", ri, curr, n_events, lp_per_event);
- }
-
- // R9 thresholds
- double threshold = model_stdv() ? 8.0f : 4.0f; // TODO: check
- if(fabs(lp_per_event) < threshold) {
- out_rs.push_back(rs);
- }
- }
- input.swap(out_rs);
-}
-
-std::string join_sequences_at_kmer(const std::string& a, const std::string& b, const uint32_t k)
-{
- // this is a special case to make the calling code cleaner
- if(a.empty())
- return b;
-
- // These sequences must have a k-mer match at the start/end
- std::string a_last_kmer = a.substr(a.size() - k);
- std::string b_last_kmer = b.substr(0, k);
- assert(a_last_kmer == b_last_kmer);
- return a + b.substr(k);
-}
-
-void run_splice_segment(HMMRealignmentInput& window, uint32_t segment_id, const uint32_t k)
-{
- // The structure of the data looks like this:
-
- // --------------------------------------------------------
- // S M E
- // where is the start column, M is the middle column and E
- // is the end column. We want to call a new consensus from S
- // to E. We do this by generating the base sequence from S to E
- // and then applying all of the alternatives indicated by the
- // start and middle column. We score these alternatives using
- // the read strands spanning from S to E. After a new consensus
- // has been selected, we re-calculate the alignments of events to
- // the middle anchor.
-
- // Get the segments
- assert(segment_id + 2 < window.anchored_columns.size());
- HMMAnchoredColumn& start_column = window.anchored_columns[segment_id];
- HMMAnchoredColumn& middle_column = window.anchored_columns[segment_id + 1];
- HMMAnchoredColumn& end_column = window.anchored_columns[segment_id + 2];
-
- std::string s_m_base = start_column.base_sequence;
- std::string m_e_base = middle_column.base_sequence;
-
- // The collection of alternative sequences
- std::vector<std::string> alts;
-
- for(uint32_t ai = 0; ai < start_column.alt_sequences.size(); ++ai) {
- alts.push_back(start_column.alt_sequences[ai]);
- }
-
- // set up the input data for the HMM
- std::vector<HMMInputData> data = get_input_for_columns(window, start_column, end_column);
-
- if(opt::verbose > 0) {
- fprintf(stderr, "correcting segment %u with %zu reads\n", segment_id, data.size());
- }
-
- // The current consensus sequence
- std::string original = join_sequences_at_kmer(s_m_base, m_e_base, k);
- std::string base = original;
-
- // filter out poor quality reads
- filter_outlier_data(data, base);
-
- // Only attempt correction if there are any reads here
- if(!data.empty()) {
-
- std::string bs_result = run_block_substitution(base, data, alts);
- std::string mut_result = run_mutation(bs_result, data);
- base = mut_result;
- }
-
- if(opt::verbose > 0) {
- fprintf(stderr, "ORIGINAL[%d] %s\n", segment_id, original.c_str());
- fprintf(stderr, "RESULT[%d] %s\n", segment_id, base.c_str());
- }
-
- // Update the sequences for the start and middle segments
- // by cutting the new consensus in the middle
- // We maintain the k-mer match invariant by requiring the
- // sequences to overlap by k-bp
- assert(base.length() >= k);
- uint32_t midpoint_kmer = (base.length() - k + 1) / 2;
-
- std::string s_m_fixed = base.substr(0, midpoint_kmer + k);
- std::string m_e_fixed = base.substr(midpoint_kmer);
-
- assert(s_m_fixed.substr(s_m_fixed.size() - k) == m_e_fixed.substr(0, k));
-
- start_column.base_sequence = s_m_fixed;
- middle_column.base_sequence = m_e_fixed;
-
- // Update the event indices in the first column to match
- for(uint32_t ri = 0; ri < data.size(); ++ri) {
-
- // Realign to the consensus sequence
- std::vector<HMMAlignmentState> decodes = profile_hmm_align(base, data[ri]);
-
- // Get the closest event aligned to the target kmer
- int32_t min_k_dist = base.length();
- uint32_t event_idx = 0;
- for(uint32_t di = 0; di < decodes.size(); ++di) {
- int32_t dist = abs((int)decodes[di].kmer_idx - (int)midpoint_kmer);
- if(dist <= min_k_dist) {
- min_k_dist = dist;
- event_idx = decodes[di].event_idx;
- }
- }
-
- middle_column.anchors[data[ri].anchor_index].event_idx = event_idx;
- }
-}
-
-// update the training data on the current segment
-void train_segment(HMMRealignmentInput& window, uint32_t segment_id)
-{
- // Get the segments
- assert(segment_id + 2 < window.anchored_columns.size());
- HMMAnchoredColumn& start_column = window.anchored_columns[segment_id];
- HMMAnchoredColumn& middle_column = window.anchored_columns[segment_id + 1];
- HMMAnchoredColumn& end_column = window.anchored_columns[segment_id + 2];
-
- std::string s_m_base = start_column.base_sequence;
- std::string m_e_base = middle_column.base_sequence;
-
- // Set up the the input data for the HMM
- std::vector<HMMInputData> input = get_input_for_columns(window, start_column, end_column);
-
- // no training can be performed if there are no reads for this segment
- if(input.empty()) {
- return;
- }
-
- // assume models for all the reads have the same k
- const uint32_t k = input[0].read->pore_model[input[0].strand].k;
-
- std::string segment_sequence = join_sequences_at_kmer(s_m_base, m_e_base, k);
-
- for(uint32_t ri = 0; ri < input.size(); ++ri) {
- std::vector<HMMAlignmentState> decodes = profile_hmm_align(segment_sequence, input[ri]);
- update_training_with_segment(segment_sequence, input[ri]);
- }
-}
-
-void train(HMMRealignmentInput& window)
-{
- // train on current consensus
- uint32_t num_segments = window.anchored_columns.size();
- for(uint32_t segment_id = 0; segment_id < num_segments - 2; ++segment_id) {
- train_segment(window, segment_id);
- }
-
- // Update model parameters
- for(uint32_t ri = 0; ri < window.reads.size(); ++ri) {
- window.reads[ri]->parameters[0].train();
- window.reads[ri]->parameters[1].train();
- }
-}
-
-std::string call_consensus_for_window(const Fast5Map& name_map, const std::string& contig, int start_base, int end_base)
-{
- const int minor_segment_stride = 50;
- HMMRealignmentInput window = build_input_for_region(opt::bam_file,
- opt::genome_file,
- name_map,
- contig,
- start_base,
- end_base,
- minor_segment_stride);
- uint32_t num_segments = window.anchored_columns.size();
-
- // If there are not reads or not enough segments do not try to call a consensus sequence
- if(window.reads.empty() || num_segments < 3) {
- // No data for this window, just return the original sequence as the consensus
- assert(!window.original_sequence.empty());
- return window.original_sequence;
- }
-
- if(opt::verbose > 0) {
- fprintf(stderr, "correcting window %s:%d-%d with %zu reads\n", contig.c_str(), start_base, end_base, window.reads.size());
- }
-
- //
- // Train the HMM
- //
- WARN_ONCE("Debug: using default transition parameters");
- //train(window);
-
- // assume models for all the reads have the same k
- const uint32_t k = window.reads[0]->pore_model[T_IDX].k;
-
- //
- // Compute the new consensus sequence
- //
- std::string reference = "";
- std::string consensus = "";
-
- uint32_t start_segment_id = 0;
-
- // Copy the base segments before they are updated
- // by the consensus algorithm
- std::vector<std::string> ref_segments;
- for(uint32_t segment_id = 0; segment_id < num_segments; ++segment_id) {
- ref_segments.push_back(window.anchored_columns[segment_id].base_sequence);
- }
-
- // Initialize progress status
- std::stringstream message;
- message << "[consensus] " << contig << ":" << start_base << "-" << end_base;
- Progress progress(message.str());
-
- for(uint32_t segment_id = start_segment_id; segment_id < num_segments - 2; ++segment_id) {
-
- // update progress
- if(opt::show_progress) {
- progress.print((float)segment_id / (num_segments - 2));
- }
-
- // run the consensus algorithm for this segment
- run_splice_segment(window, segment_id, k);
-
- // run_splice_segment updates the base_sequence of the current anchor, grab it and append
- std::string base = window.anchored_columns[segment_id].base_sequence;
-
- // append the new sequences in, respecting the K overlap
- reference = join_sequences_at_kmer(reference, ref_segments[segment_id], k);
- consensus = join_sequences_at_kmer(consensus, base, k);
-
- if(opt::verbose > 0) {
- fprintf(stderr, "UNCORRECT[%d]: %s\n", segment_id, reference.c_str());
- fprintf(stderr, "CONSENSUS[%d]: %s\n", segment_id, consensus.c_str());
- }
- }
-
- // Append segment that ends at the last anchor
- reference = join_sequences_at_kmer(reference, ref_segments[num_segments - 2], k);
- const std::string& last_segment =
- window.anchored_columns[num_segments - 2].base_sequence;
- consensus = join_sequences_at_kmer(consensus, last_segment, k);
-
- if(opt::show_progress) {
- progress.end();
- }
-
- return consensus;
-}
-
-void parse_consensus_options(int argc, char** argv)
-{
- bool die = false;
- for (char c; (c = getopt_long(argc, argv, shortopts, longopts, NULL)) != -1;) {
- std::istringstream arg(optarg != NULL ? optarg : "");
- switch (c) {
- case 'r': arg >> opt::reads_file; break;
- case 'g': arg >> opt::genome_file; break;
- case 'b': arg >> opt::bam_file; break;
- case 'w': arg >> opt::window; break;
- case 'o': arg >> opt::output_file; break;
- case '?': die = true; break;
- case 't': arg >> opt::num_threads; break;
- case 'v': opt::verbose++; break;
- case OPT_MODELS_FOFN: arg >> opt::models_fofn; break;
- case OPT_PROGRESS: opt::show_progress = 1; break;
- case OPT_HELP:
- std::cout << CONSENSUS_USAGE_MESSAGE;
- exit(EXIT_SUCCESS);
- case OPT_VERSION:
- std::cout << CONSENSUS_VERSION_MESSAGE;
- exit(EXIT_SUCCESS);
- }
- }
-
- if (argc - optind < 0) {
- std::cerr << SUBPROGRAM ": missing arguments\n";
- die = true;
- } else if (argc - optind > 0) {
- std::cerr << SUBPROGRAM ": too many arguments\n";
- die = true;
- }
-
- if(opt::num_threads <= 0) {
- std::cerr << SUBPROGRAM ": invalid number of threads: " << opt::num_threads << "\n";
- die = true;
- }
-
- if(opt::reads_file.empty()) {
- std::cerr << SUBPROGRAM ": a --reads file must be provided\n";
- die = true;
- }
-
- if(opt::genome_file.empty()) {
- std::cerr << SUBPROGRAM ": a --genome file must be provided\n";
- die = true;
- }
-
- if(opt::bam_file.empty()) {
- std::cerr << SUBPROGRAM ": a --bam file must be provided\n";
- die = true;
- }
-
- if(opt::models_fofn.empty()) {
- std::cerr << SUBPROGRAM ": a --models file must be provided\n";
- die = true;
- } else {
- // initialize the model set from the fofn
- PoreModelSet::initialize(opt::models_fofn);
- }
-
- if(opt::window.empty()) {
- std::cerr << SUBPROGRAM ": the -w (window) parameter must be provided\n";
- die = true;
- }
-
- if (die)
- {
- std::cout << "\n" << CONSENSUS_USAGE_MESSAGE;
- exit(EXIT_FAILURE);
- }
-}
-
-int consensus_main(int argc, char** argv)
-{
- parse_consensus_options(argc, argv);
- omp_set_num_threads(opt::num_threads);
-
- Fast5Map name_map(opt::reads_file);
-
- // Parse the window string
- // Replace ":" and "-" with spaces to make it parseable with stringstream
- std::replace(opt::window.begin(), opt::window.end(), ':', ' ');
- std::replace(opt::window.begin(), opt::window.end(), '-', ' ');
-
- const int WINDOW_LENGTH = 10000;
- const int WINDOW_OVERLAP = 200;
-
- std::stringstream parser(opt::window);
- std::string contig;
- int start_window_id;
- int end_window_id;
-
- parser >> contig >> start_window_id >> end_window_id;
-
- FILE* out_fp = NULL;
-
- if(!opt::output_file.empty()) {
- out_fp = fopen(opt::output_file.c_str(), "w");
- } else {
- out_fp = stdout;
- }
-
- for(int window_id = start_window_id; window_id < end_window_id; ++window_id) {
- int start_base = window_id * WINDOW_LENGTH;
- int end_base = start_base + WINDOW_LENGTH + WINDOW_OVERLAP;
-
- std::string window_consensus = call_consensus_for_window(name_map, contig, start_base, end_base);
- fprintf(out_fp, ">%s:%d\n%s\n", contig.c_str(), window_id, window_consensus.c_str());
- }
-
- if(out_fp != stdout) {
- fclose(out_fp);
- }
- return 0;
-}
diff --git a/src/nanopolish_consensus.h b/src/nanopolish_consensus.h
deleted file mode 100644
index 2cefc03..0000000
--- a/src/nanopolish_consensus.h
+++ /dev/null
@@ -1,14 +0,0 @@
-//---------------------------------------------------------
-// Copyright 2015 Ontario Institute for Cancer Research
-// Written by Jared Simpson (jared.simpson at oicr.on.ca)
-//---------------------------------------------------------
-//
-// nanopolish_consensus.cpp -- compute a new consensus sequence
-// for an assembly
-//
-#ifndef NANOPOLISH_CONSENSUS_H
-#define NANOPOLISH_CONSENSUS_H
-
-int consensus_main(int argc, char** argv);
-
-#endif
diff --git a/src/nanopolish_extract.cpp b/src/nanopolish_extract.cpp
index 9abac1b..e565865 100644
--- a/src/nanopolish_extract.cpp
+++ b/src/nanopolish_extract.cpp
@@ -14,6 +14,7 @@
#include <sstream>
#include <getopt.h>
+#include "nanopolish_read_db.h"
#include "nanopolish_extract.h"
#include "nanopolish_common.h"
#include "fs_support.hpp"
@@ -39,7 +40,7 @@ static const char *EXTRACT_USAGE_MESSAGE =
" (default: 2d-or-template)\n"
" -b, --basecaller=NAME[:VERSION] consider only data produced by basecaller NAME,\n"
" optionally with given exact VERSION\n"
-" -o, --output=FILE write output to FILE (default: stdout)\n"
+" -o, --output=FILE write output to FILE\n"
"\nReport bugs to " PACKAGE_BUGREPORT "\n\n";
namespace opt
@@ -320,6 +321,12 @@ void parse_extract_options(int argc, char** argv)
auto default_level = (int)logger::warning + opt::verbose;
logger::Logger::set_default_level(default_level);
logger::Logger::set_levels_from_options(log_level, &std::clog);
+
+ if(opt::output_file.empty()) {
+ std::cerr << SUBPROGRAM ": an output file (-o) is required\n";
+ die = true;
+ }
+
// parse paths to process
while (optind < argc)
{
@@ -330,11 +337,13 @@ void parse_extract_options(int argc, char** argv)
}
opt::paths.push_back(path);
}
+
if (opt::paths.empty())
{
std::cerr << SUBPROGRAM ": no paths to process\n";
die = true;
}
+
// check read type
if (not (opt::read_type == "template"
or opt::read_type == "complement"
@@ -345,12 +354,14 @@ void parse_extract_options(int argc, char** argv)
std::cerr << SUBPROGRAM ": invalid read type: " << opt::read_type << "\n";
die = true;
}
+
// die if errors
if (die)
{
std::cerr << "\n" << EXTRACT_USAGE_MESSAGE;
exit(EXIT_FAILURE);
}
+
LOG(info) << "paths: " << alg::os_join(opt::paths, " ") << "\n";
LOG(info) << "recurse: " << (opt::recurse? "yes" : "no") << "\n";
LOG(info) << "read_type: " << opt::read_type << "\n";
@@ -361,20 +372,26 @@ void parse_extract_options(int argc, char** argv)
int extract_main(int argc, char** argv)
{
parse_extract_options(argc, argv);
- std::ofstream ofs;
- if (not opt::output_file.empty())
+
+ // Iterate over fast5 collection extracting the sequence reads
+ // We do this in a block so the file is automatically closed
+ // when ofs goes out of scope.
{
+ std::ofstream ofs;
ofs.open(opt::output_file);
os_p = &ofs;
+
+ for (unsigned i = 0; i < opt::paths.size(); ++i)
+ {
+ process_path(opt::paths[i]);
+ }
}
- else
- {
- os_p = &std::cout;
- }
- for (unsigned i = 0; i < opt::paths.size(); ++i)
- {
- process_path(opt::paths[i]);
- }
+
+ // Build the ReadDB from the output file
+ ReadDB read_db;
+ read_db.build(opt::output_file);
+ read_db.save();
+
std::clog << "[extract] found " << opt::total_files_count
<< " files, extracted " << opt::total_files_used_count
<< " reads\n";
diff --git a/src/nanopolish_getmodel.cpp b/src/nanopolish_getmodel.cpp
index 014427e..06f8b06 100644
--- a/src/nanopolish_getmodel.cpp
+++ b/src/nanopolish_getmodel.cpp
@@ -89,13 +89,13 @@ void parse_getmodel_options(int argc, char** argv)
die = true;
}
- opt::input_file = argv[optind++];
-
if (die)
{
std::cout << "\n" << GETMODEL_USAGE_MESSAGE;
exit(EXIT_FAILURE);
}
+
+ opt::input_file = argv[optind++];
}
int getmodel_main(int argc, char** argv)
diff --git a/src/nanopolish_haplotype.cpp b/src/nanopolish_haplotype.cpp
index 9d2d91d..58019cc 100644
--- a/src/nanopolish_haplotype.cpp
+++ b/src/nanopolish_haplotype.cpp
@@ -24,6 +24,10 @@ Haplotype::Haplotype(const std::string& ref_name,
m_coordinate_map[i] = m_ref_position + i;
}
}
+
+Haplotype::~Haplotype()
+{
+}
//
bool Haplotype::apply_variant(const Variant& v)
@@ -100,6 +104,12 @@ Haplotype Haplotype::substr_by_reference(size_t start, size_t end) const
}
assert(derived_base_start != m_coordinate_map.size());
+
+ // JTS: temporary debug dump for #199
+ if(derived_base_end == m_coordinate_map.size()) {
+ print_debug_info();
+ }
+
assert(derived_base_end != m_coordinate_map.size());
assert(m_coordinate_map[derived_base_start] <= start);
assert(m_coordinate_map[derived_base_end] >= end);
@@ -157,3 +167,13 @@ size_t Haplotype::_find_derived_index_by_ref_lower_bound(size_t ref_index) const
return m_coordinate_map.size();
}
+void Haplotype::print_debug_info() const
+{
+ fprintf(stderr, "[haplotype-debug] ctg: %s, position: %d\n", m_ref_name.c_str(), m_ref_position);
+ fprintf(stderr, "[haplotype-debug] r-sequence: %s\n", m_reference.c_str());
+ fprintf(stderr, "[haplotype-debug] h-sequence: %s\n", m_sequence.c_str());
+ for(size_t i = 0; i < m_variants.size(); ++i) {
+ fprintf(stderr, "[haplotype-debug] variant[%zu]: ", i);
+ m_variants[i].write_vcf(stderr);
+ }
+}
diff --git a/src/nanopolish_haplotype.h b/src/nanopolish_haplotype.h
index fa29ad7..9bd0861 100644
--- a/src/nanopolish_haplotype.h
+++ b/src/nanopolish_haplotype.h
@@ -19,6 +19,9 @@ class Haplotype
Haplotype(const std::string& ref_name,
const size_t ref_position,
const std::string& ref_sequence);
+
+ //
+ ~Haplotype();
// get the sequence of the haplotype
const std::string& get_sequence() const { return m_sequence; }
@@ -70,6 +73,8 @@ class Haplotype
// This mimics std::lower_bound
size_t _find_derived_index_by_ref_lower_bound(size_t ref_index) const;
+ void print_debug_info() const;
+
//
// data
//
diff --git a/src/nanopolish_index.cpp b/src/nanopolish_index.cpp
new file mode 100644
index 0000000..b997c3c
--- /dev/null
+++ b/src/nanopolish_index.cpp
@@ -0,0 +1,197 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//
+// nanopolish index - build an index from FASTA/FASTQ file
+// and the associated signal-level data
+//
+//---------------------------------------------------------
+//
+
+//
+// Getopt
+//
+#define SUBPROGRAM "index"
+
+#include <iostream>
+#include <sstream>
+#include <getopt.h>
+
+#include "nanopolish_index.h"
+#include "nanopolish_common.h"
+#include "nanopolish_read_db.h"
+#include "fs_support.hpp"
+#include "logger.hpp"
+#include "fast5.hpp"
+#include "profiler.h"
+
+static const char *INDEX_VERSION_MESSAGE =
+SUBPROGRAM " Version " PACKAGE_VERSION "\n"
+"Written by Jared Simpson.\n"
+"\n"
+"Copyright 2017 Ontario Institute for Cancer Research\n";
+
+static const char *INDEX_USAGE_MESSAGE =
+"Usage: " PACKAGE_NAME " " SUBPROGRAM " [OPTIONS] -d nanopore_raw_file_directory reads.fastq\n"
+"Build an index mapping from basecalled reads to the signals measured by the sequencer\n"
+"\n"
+" --help display this help and exit\n"
+" --version display version\n"
+" -v, --verbose display verbose output\n"
+" -d, --directory path to the directory containing the raw ONT signal files\n"
+" -f, --fast5-fofn file containing the paths to each fast5 file for the run\n"
+"\nReport bugs to " PACKAGE_BUGREPORT "\n\n";
+
+namespace opt
+{
+ static unsigned int verbose = 0;
+ static std::string raw_file_directory;
+ static std::string fast5_fofn;
+ static std::string reads_file;
+}
+static std::ostream* os_p;
+
+void index_file(ReadDB& read_db, const std::string& fn)
+{
+ PROFILE_FUNC("index_file")
+ fast5::File* fp = new fast5::File(fn);
+ if(fp->is_open()) {
+ fast5::Raw_Samples_Params params = fp->get_raw_samples_params();
+ std::string read_id = params.read_id;
+ read_db.add_signal_path(read_id, fn);
+ }
+ delete fp;
+} // process_file
+
+void index_path(ReadDB& read_db, const std::string& path)
+{
+ fprintf(stderr, "Indexing %s\n", path.c_str());
+ if (is_directory(path)) {
+ auto dir_list = list_directory(path);
+ for (const auto& fn : dir_list) {
+ if(fn == "." or fn == "..") {
+ continue;
+ }
+
+ std::string full_fn = path + "/" + fn;
+ if(is_directory(full_fn)) {
+ // recurse
+ index_path(read_db, full_fn);
+ read_db.print_stats();
+ } else if (full_fn.find(".fast5") != -1 && fast5::File::is_valid_file(full_fn)) {
+ index_file(read_db, full_fn);
+ }
+ }
+ }
+} // process_path
+
+void process_fast5_fofn(ReadDB& read_db, const std::string& fast5_fofn)
+{
+ //
+ std::ifstream in_file(fast5_fofn.c_str());
+ if(in_file.bad()) {
+ fprintf(stderr, "error: could not file %s\n", fast5_fofn.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ // read
+ size_t count = 0;
+ std::string filename;
+ while(getline(in_file, filename)) {
+ if(count++ % 1000 == 0) {
+ fprintf(stderr, "Loaded %zu from fofn\n", count);
+ }
+ index_file(read_db, filename);
+ }
+}
+
+static const char* shortopts = "vd:f:";
+
+enum {
+ OPT_HELP = 1,
+ OPT_VERSION,
+ OPT_LOG_LEVEL,
+};
+
+static const struct option longopts[] = {
+ { "help", no_argument, NULL, OPT_HELP },
+ { "version", no_argument, NULL, OPT_VERSION },
+ { "log-level", required_argument, NULL, OPT_LOG_LEVEL },
+ { "verbose", no_argument, NULL, 'v' },
+ { "directory", required_argument, NULL, 'd' },
+ { "fast5-fofn", required_argument, NULL, 'f' },
+ { NULL, 0, NULL, 0 }
+};
+
+void parse_index_options(int argc, char** argv)
+{
+ bool die = false;
+ std::vector< std::string> log_level;
+ for (char c; (c = getopt_long(argc, argv, shortopts, longopts, NULL)) != -1;) {
+ std::istringstream arg(optarg != NULL ? optarg : "");
+ switch (c) {
+ case OPT_HELP:
+ std::cout << INDEX_USAGE_MESSAGE;
+ exit(EXIT_SUCCESS);
+ case OPT_VERSION:
+ std::cout << INDEX_VERSION_MESSAGE;
+ exit(EXIT_SUCCESS);
+ case OPT_LOG_LEVEL:
+ log_level.push_back(arg.str());
+ break;
+ case 'v': opt::verbose++; break;
+ case 'd': arg >> opt::raw_file_directory; break;
+ case 'f': arg >> opt::fast5_fofn; break;
+ }
+ }
+
+ // set log levels
+ auto default_level = (int)logger::warning + opt::verbose;
+ logger::Logger::set_default_level(default_level);
+ logger::Logger::set_levels_from_options(log_level, &std::clog);
+
+ if (argc - optind < 1) {
+ std::cerr << SUBPROGRAM ": not enough arguments\n";
+ die = true;
+ }
+
+ if (argc - optind > 1) {
+ std::cerr << SUBPROGRAM ": too many arguments\n";
+ die = true;
+ }
+
+ if (die)
+ {
+ std::cout << "\n" << INDEX_USAGE_MESSAGE;
+ exit(EXIT_FAILURE);
+ }
+
+ opt::reads_file = argv[optind++];
+}
+
+int index_main(int argc, char** argv)
+{
+ parse_index_options(argc, argv);
+
+ // import read names, and possibly fast5 paths, from the fasta/fastq file
+ ReadDB read_db;
+ read_db.build(opt::reads_file);
+
+ bool all_reads_have_paths = read_db.check_signal_paths();
+
+ // if the input fastq did not contain a complete set of paths
+ // use the fofn/directory provided to augment the index
+ if(!all_reads_have_paths) {
+ if(!opt::raw_file_directory.empty()) {
+ index_path(read_db, opt::raw_file_directory);
+ }
+
+ if(!opt::fast5_fofn.empty()) {
+ process_fast5_fofn(read_db, opt::fast5_fofn);
+ }
+ }
+
+ read_db.print_stats();
+ read_db.save();
+ return 0;
+}
diff --git a/src/nanopolish_index.h b/src/nanopolish_index.h
new file mode 100644
index 0000000..fae2dfb
--- /dev/null
+++ b/src/nanopolish_index.h
@@ -0,0 +1,14 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//---------------------------------------------------------
+//
+#ifndef NANOPOLISH_INDEX_H
+#define NANOPOLISH_INDEX_H
+
+#include <string>
+#include <vector>
+
+int index_main(int argc, char** argv);
+
+#endif
diff --git a/src/nanopolish_methyltrain.cpp b/src/nanopolish_methyltrain.cpp
index b797364..64829b5 100644
--- a/src/nanopolish_methyltrain.cpp
+++ b/src/nanopolish_methyltrain.cpp
@@ -32,9 +32,9 @@
#include "nanopolish_matrix.h"
#include "nanopolish_profile_hmm.h"
#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
#include "nanopolish_model_names.h"
#include "nanopolish_pore_model_set.h"
+#include "nanopolish_read_db.h"
#include "training_core.hpp"
#include "H5pubconf.h"
#include "profiler.h"
@@ -310,7 +310,7 @@ bool recalibrate_model(SquiggleRead &sr,
}
// Update the training data with aligned events from a read
-void add_aligned_events(const Fast5Map& name_map,
+void add_aligned_events(const ReadDB& read_db,
const faidx_t* fai,
const bam_hdr_t* hdr,
const bam1_t* record,
@@ -325,10 +325,9 @@ void add_aligned_events(const Fast5Map& name_map,
{
// Load a squiggle read for the mapped read
std::string read_name = bam_get_qname(record);
- std::string fast5_path = name_map.get_path(read_name);
// load read
- SquiggleRead sr(read_name, fast5_path);
+ SquiggleRead sr(read_name, read_db);
// replace the models that are built into the read with the model we are training
sr.replace_models(training_kit, training_alphabet, training_k);
@@ -698,7 +697,7 @@ TrainingResult retrain_model_from_events(const PoreModel& current_model,
return result;
}
-void train_one_round(const Fast5Map& name_map,
+void train_one_round(const ReadDB& read_db,
const std::string& kit_name,
const std::string& alphabet,
size_t k,
@@ -781,7 +780,7 @@ void train_one_round(const Fast5Map& name_map,
bam1_t* record = records[i];
size_t read_idx = num_reads_realigned + i;
if( (record->core.flag & BAM_FUNMAP) == 0) {
- add_aligned_events(name_map, fai, hdr, record, read_idx,
+ add_aligned_events(read_db, fai, hdr, record, read_idx,
clip_start, clip_end,
kit_name, alphabet, k,
round, model_training_data);
@@ -873,7 +872,8 @@ int methyltrain_main(int argc, char** argv)
parse_methyltrain_options(argc, argv);
omp_set_num_threads(opt::num_threads);
- Fast5Map name_map(opt::reads_file);
+ ReadDB read_db;
+ read_db.load(opt::reads_file);
// Import the models to train into the pore model set
assert(!opt::models_fofn.empty());
@@ -890,7 +890,7 @@ int methyltrain_main(int argc, char** argv)
for(size_t round = 0; round < opt::num_training_rounds; round++) {
fprintf(stderr, "Starting round %zu\n", round);
- train_one_round(name_map, training_kit, mtrain_alphabet->get_name(), training_k, round);
+ train_one_round(read_db, training_kit, mtrain_alphabet->get_name(), training_k, round);
/*
if(opt::write_models) {
write_models(training_kit, mtrain_alphabet->get_name(), training_k, round);
diff --git a/src/nanopolish_phase_reads.cpp b/src/nanopolish_phase_reads.cpp
index 43c62e3..a3e9a88 100644
--- a/src/nanopolish_phase_reads.cpp
+++ b/src/nanopolish_phase_reads.cpp
@@ -28,13 +28,13 @@
#include "nanopolish_poremodel.h"
#include "nanopolish_transition_parameters.h"
#include "nanopolish_profile_hmm.h"
-#include "nanopolish_fast5_map.h"
#include "nanopolish_pore_model_set.h"
#include "nanopolish_variant.h"
#include "nanopolish_haplotype.h"
#include "nanopolish_alignment_db.h"
#include "nanopolish_bam_processor.h"
#include "nanopolish_bam_utils.h"
+#include "nanopolish_index.h"
#include "H5pubconf.h"
#include "profiler.h"
#include "progress.h"
@@ -173,7 +173,7 @@ void parse_phase_reads_options(int argc, char** argv)
}
}
-void phase_single_read(const Fast5Map& name_map,
+void phase_single_read(const ReadDB& read_db,
const faidx_t* fai,
const std::vector<Variant>& variants,
samFile* sam_fp,
@@ -190,10 +190,9 @@ void phase_single_read(const Fast5Map& name_map,
// Load a squiggle read for the mapped read
std::string read_name = bam_get_qname(record);
- std::string fast5_path = name_map.get_path(read_name);
// load read
- SquiggleRead sr(read_name, fast5_path);
+ SquiggleRead sr(read_name, read_db);
std::string ref_name = hdr->target_name[record->core.tid];
int alignment_start_pos = record->core.pos;
@@ -341,7 +340,8 @@ int phase_reads_main(int argc, char** argv)
parse_phase_reads_options(argc, argv);
omp_set_num_threads(opt::num_threads);
- Fast5Map name_map(opt::reads_file);
+ ReadDB read_db;
+ read_db.load(opt::reads_file);
// load reference fai file
faidx_t *fai = fai_load(opt::genome_file.c_str());
@@ -371,7 +371,7 @@ int phase_reads_main(int argc, char** argv)
// the BamProcessor framework calls the input function with the
// bam record, read index, etc passed as parameters
// bind the other parameters the worker function needs here
- auto f = std::bind(phase_single_read, name_map, fai, std::ref(variants), sam_out, _1, _2, _3, _4, _5);
+ auto f = std::bind(phase_single_read, std::ref(read_db), std::ref(fai), std::ref(variants), sam_out, _1, _2, _3, _4, _5);
BamProcessor processor(opt::bam_file, opt::region, opt::num_threads);
// Copy the bam header to std
diff --git a/src/nanopolish_raw_loader.cpp b/src/nanopolish_raw_loader.cpp
new file mode 100644
index 0000000..33f772a
--- /dev/null
+++ b/src/nanopolish_raw_loader.cpp
@@ -0,0 +1,282 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//---------------------------------------------------------
+//
+// nanopolish_raw_loader - utilities and helpers for loading
+// data directly from raw nanopore files without events
+//
+#include "nanopolish_raw_loader.h"
+#include "nanopolish_profile_hmm.h"
+
+//
+void estimate_scalings_using_mom(const std::string& sequence,
+ const PoreModel& pore_model,
+ const event_table& et,
+ double& out_shift,
+ double& out_scale)
+{
+ size_t k = pore_model.k;
+ size_t n_kmers = sequence.size() - k + 1;
+ const Alphabet* alphabet = pore_model.pmalphabet;
+
+ // Calculate summary statistics over the events and
+ // the model implied by the read
+ double event_level_sum = 0.0f;
+ for(size_t i = 0; i < et.n; ++i) {
+ event_level_sum += et.event[i].mean;
+ }
+
+ double kmer_level_sum = 0.0f;
+ double kmer_level_sq_sum = 0.0f;
+ for(size_t i = 0; i < n_kmers; ++i) {
+ size_t kmer_rank = alphabet->kmer_rank(sequence.substr(i, k).c_str(), k);
+ double l = pore_model.get_parameters(kmer_rank).level_mean;
+ kmer_level_sum += l;
+ kmer_level_sq_sum += pow(l, 2.0f);
+ }
+ out_shift = event_level_sum / et.n - kmer_level_sum / n_kmers;
+
+ // estimate scale
+ double event_level_sq_sum = 0.0f;
+ for(size_t i = 0; i < et.n; ++i) {
+ event_level_sq_sum += pow(et.event[i].mean - out_shift, 2.0);
+ }
+
+ out_scale = (event_level_sq_sum / et.n) / (kmer_level_sq_sum / n_kmers);
+
+#if DEBUG_PRINT_STATS
+ fprintf(stderr, "event mean: %.2lf kmer mean: %.2lf shift: %.2lf\n", event_level_sum / et.n, kmer_level_sum / n_kmers, out_shift);
+ fprintf(stderr, "event sq-mean: %.2lf kmer sq-mean: %.2lf scale: %.2lf\n", event_level_sq_sum / et.n, kmer_level_sq_sum / n_kmers, out_scale);
+ fprintf(stderr, "truth shift: %.2lf scale: %.2lf\n", pore_model.shift, pore_model.scale);
+#endif
+}
+
+std::vector<AlignedPair> banded_simple_event_align(SquiggleRead& read, const std::string& sequence)
+{
+ size_t strand_idx = 0;
+ size_t k = read.pore_model[strand_idx].k;
+ const Alphabet* alphabet = read.pore_model[strand_idx].pmalphabet;
+
+#if DEBUG_PRINT_STATS
+ // Build a debug event->kmer map
+ std::vector<size_t> kmer_for_event(read.events[strand_idx].size());
+ for(size_t ki = 0; ki < read.base_to_event_map.size(); ++ki) {
+ IndexPair& elem = read.base_to_event_map[ki].indices[0];
+ if(elem.start == -1) {
+ continue;
+ }
+
+ for(int j = elem.start; j <= elem.stop; ++j) {
+ if(j >= 0 && j < kmer_for_event.size()) {
+ kmer_for_event[j] = ki;
+ }
+ }
+ }
+#endif
+
+ const uint8_t FROM_D = 0;
+ const uint8_t FROM_U = 1;
+ const uint8_t FROM_L = 2;
+
+ // qc
+ double min_average_log_emission = -5.0;
+
+ // banding
+ int bandwidth = 500;
+ int half_band = bandwidth / 2;
+
+ // transitions
+ double lp_skip = log(0.001);
+ double lp_stay = log(0.5);
+ double lp_step = log(1.0 - exp(lp_skip) - exp(lp_stay));
+ double lp_trim = log(0.1);
+
+ size_t n_events = read.events[strand_idx].size();
+ size_t n_kmers = sequence.size() - k + 1;
+
+ // Calculate the minimum event index that is within the band for each read kmer
+ // We determine this using the expected number of events observed per kmer
+ double events_per_kmer = (double)n_events / n_kmers;
+ std::vector<int> min_event_idx_by_kmer(n_kmers);
+ for(size_t ki = 0; ki < n_kmers; ++ki) {
+ int expected_event_idx = (double)(ki * events_per_kmer);
+ min_event_idx_by_kmer[ki] = std::max(expected_event_idx - half_band, 0);
+ }
+ // Initialize DP matrices
+ DoubleMatrix viterbi_matrix;
+ UInt8Matrix backtrack_matrix;
+
+ size_t n_rows = bandwidth;
+ size_t n_cols = n_kmers + 1;
+ allocate_matrix(viterbi_matrix, n_rows, n_cols);
+ allocate_matrix(backtrack_matrix, n_rows, n_cols);
+
+ for(size_t i = 0; i < bandwidth; ++i) {
+ set(viterbi_matrix, i, 0, i * lp_trim);
+ set(backtrack_matrix, i, 0, 0);
+ }
+
+ // Fill in the matrix
+ for(int col = 1; col < n_cols; ++col) {
+ int kmer_idx = col - 1;
+ int min_event_idx = min_event_idx_by_kmer[kmer_idx];
+ int min_event_idx_prev_col = kmer_idx > 0 ? min_event_idx_by_kmer[kmer_idx - 1] : 0;
+ size_t kmer_rank = alphabet->kmer_rank(sequence.substr(kmer_idx, k).c_str(), k);
+
+ for(int row = 0; row < n_rows; ++row) {
+
+ int event_idx = min_event_idx + row;
+ if(event_idx >= n_events) {
+ set(viterbi_matrix, row, col, -INFINITY);
+ continue;
+ }
+
+ // dp update
+ // here we are calculating whether the event for each neighboring cell is within the band
+ // and calculating its position within the column
+ int row_up = event_idx - min_event_idx - 1;
+ int row_diag = event_idx - min_event_idx_prev_col - 1;
+ int row_left = event_idx - min_event_idx_prev_col;
+
+ double up = row_up >= 0 && row_up < n_rows ? get(viterbi_matrix, row_up, col) : -INFINITY;
+ double diag = row_diag >= 0 && row_diag < n_rows ? get(viterbi_matrix, row_diag, col - 1) : -INFINITY;
+ double left = row_left >= 0 && row_left < n_rows ? get(viterbi_matrix, row_left, col - 1) : -INFINITY;
+
+ float lp_emission = log_probability_match_r9(read, kmer_rank, event_idx, strand_idx);
+
+ double score_d = diag + lp_step + lp_emission;
+ double score_u = up + lp_stay + lp_emission;
+ double score_l = left + lp_skip;
+
+ double max_score = score_d;
+ uint8_t from = FROM_D;
+
+ max_score = score_u > max_score ? score_u : max_score;
+ from = max_score == score_u ? FROM_U : from;
+
+ max_score = score_l > max_score ? score_l : max_score;
+ from = max_score == score_l ? FROM_L : from;
+
+ set(viterbi_matrix, row, col, max_score);
+ set(backtrack_matrix, row, col, from);
+ }
+ }
+
+ // Backtrack
+
+ // Initialize by finding best alignment between an event and the last kmer
+ int curr_k_idx = n_kmers - 1;
+ int curr_event_idx = 0;
+ double max_score = -INFINITY;
+ for(size_t row = 0; row < n_rows; ++row) {
+ size_t col = curr_k_idx + 1;
+ int ei = row + min_event_idx_by_kmer[curr_k_idx];
+ double s = get(viterbi_matrix, row, col) + (n_events - ei - 1) * lp_trim;
+ if(s > max_score && ei < n_events) {
+ max_score = s;
+ curr_event_idx = ei;
+ }
+ }
+
+ // debug stats
+ double sum_emission = 0;
+ double n_aligned_events = 0;
+ std::vector<AlignedPair> out;
+
+#if DEBUG_PRINT_STATS
+ int sum_distance_from_debug = 0;
+ int max_distance_from_debug = 0;
+ int num_exact_matches = 0;
+ int max_distance_from_expected = 0;
+ int sum_distance_from_expected = 0;
+ std::vector<int> debug_event_for_kmer(n_kmers, -1);
+#endif
+
+ while(curr_k_idx >= 0) {
+ // emit alignment
+ out.push_back({curr_k_idx, curr_event_idx});
+
+ size_t kmer_rank = alphabet->kmer_rank(sequence.substr(curr_k_idx, k).c_str(), k);
+ sum_emission += log_probability_match_r9(read, kmer_rank, curr_event_idx, strand_idx);
+ n_aligned_events += 1;
+
+#if DEBUG_PRINT_STATS
+ // update debug stats
+ debug_event_for_kmer[curr_k_idx] = curr_event_idx;
+
+ int kd = abs(curr_k_idx - kmer_for_event[curr_event_idx]);
+ sum_distance_from_debug += kd;
+ max_distance_from_debug = std::max(kd, max_distance_from_debug);
+ num_exact_matches += curr_k_idx == kmer_for_event[curr_event_idx];
+
+ int expected_event = (double)(curr_k_idx * events_per_kmer);
+ int ed = curr_event_idx - expected_event;
+ max_distance_from_expected = std::max(abs(ed), max_distance_from_expected);
+ sum_distance_from_expected += ed;
+#endif
+
+ // update indices using backtrack pointers
+ int row = curr_event_idx - min_event_idx_by_kmer[curr_k_idx];
+ int col = curr_k_idx + 1;
+
+ uint8_t from = get(backtrack_matrix, row, col);
+ if(from == FROM_D) {
+ curr_k_idx -= 1;
+ curr_event_idx -= 1;
+ } else if(from == FROM_U) {
+ curr_event_idx -= 1;
+ } else {
+ curr_k_idx -= 1;
+ }
+ }
+ std::reverse(out.begin(), out.end());
+
+#if DEBUG_PRINT_MATRIX
+ // Columm-wise debugging
+ for(size_t ci = 1; ci < n_cols; ++ci) {
+ size_t curr_k_idx = ci - 1;
+ size_t expected_event = (double)(curr_k_idx * events_per_kmer);
+
+ // Find the maximum value in this column
+ double max_row = -INFINITY;
+ size_t max_event_idx = 0;
+
+ std::stringstream debug_out;
+ for(size_t ri = 0; ri < n_rows; ++ri) {
+ size_t curr_event_idx = ri + min_event_idx_by_kmer[curr_k_idx];
+ double s = get(viterbi_matrix, ri, ci);
+ if(s > max_row) {
+ max_row = s;
+ max_event_idx = curr_event_idx;
+ }
+ debug_out << s << ",";
+ }
+
+ std::string debug_str = debug_out.str();
+ fprintf(stderr, "DEBUG_MATRIX\t%s\n", debug_str.substr(0, debug_str.size() - 1).c_str());
+
+ size_t aligned_event = debug_event_for_kmer[curr_k_idx];
+ fprintf(stderr, "DEBUG BAND: k: %zu ee: %zu me: %zu ae: %zu d: %d\n",
+ curr_k_idx, expected_event, max_event_idx, aligned_event, (int)max_event_idx - (int)aligned_event);
+ }
+#endif
+
+ // QC results
+ double avg_log_emission = sum_emission / n_aligned_events;
+ bool spanned = out.front().ref_pos == 0 && out.back().ref_pos == n_kmers - 1;
+
+ if(avg_log_emission < min_average_log_emission || !spanned) {
+ out.clear();
+ }
+
+#if DEBUG_PRINT_STATS
+ fprintf(stderr, "events per base: %.2lf\n", events_per_kmer);
+ fprintf(stderr, "truth stats -- avg: %.2lf max: %d exact: %d\n", (double)sum_distance_from_debug / n_kmers, max_distance_from_debug, num_exact_matches);
+ fprintf(stderr, "event stats -- avg: %.2lf max: %d\n", (double)sum_distance_from_expected / n_events, max_distance_from_expected);
+ fprintf(stderr, "emission stats -- avg: %.2lf\n", sum_emission / n_aligned_events);
+#endif
+ free_matrix(viterbi_matrix);
+ free_matrix(backtrack_matrix);
+ return out;
+}
diff --git a/src/nanopolish_raw_loader.h b/src/nanopolish_raw_loader.h
new file mode 100644
index 0000000..6c8bfe3
--- /dev/null
+++ b/src/nanopolish_raw_loader.h
@@ -0,0 +1,25 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//---------------------------------------------------------
+//
+// nanopolish_raw_loader - utilities and helpers for loading
+// data directly from raw nanopore files without events
+//
+#ifndef NANOPOLISH_RAW_LOADER_H
+#define NANOPOLISH_RAW_LOADER_H
+
+#include "nanopolish_squiggle_read.h"
+#include "nanopolish_anchor.h"
+#include "scrappie_structures.h"
+
+void estimate_scalings_using_mom(const std::string& sequence,
+ const PoreModel& pore_model,
+ const event_table& et,
+ double& out_shift,
+ double& out_scale);
+
+std::vector<AlignedPair> banded_simple_event_align(SquiggleRead& read,
+ const std::string& sequence);
+
+#endif
diff --git a/src/nanopolish_read_db.cpp b/src/nanopolish_read_db.cpp
new file mode 100644
index 0000000..c8039a6
--- /dev/null
+++ b/src/nanopolish_read_db.cpp
@@ -0,0 +1,231 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//---------------------------------------------------------
+//
+// nanopolish_read_db -- database of reads and their
+// associated signal data
+//
+#include <zlib.h>
+#include <fstream>
+#include <ostream>
+#include <iostream>
+#include <sys/stat.h>
+#include "nanopolish_common.h"
+#include "htslib/kseq.h"
+#include "htslib/bgzf.h"
+#include "nanopolish_read_db.h"
+
+#define READ_DB_SUFFIX ".readdb"
+#define GZIPPED_READS_SUFFIX ".fa.gz"
+
+// Tell KSEQ what functions to use to open/read files
+KSEQ_INIT(gzFile, gzread)
+
+//
+ReadDB::ReadDB() : m_fai(NULL)
+{
+
+}
+
+//
+void ReadDB::build(const std::string& input_reads_filename)
+{
+ // generate output filename
+ m_indexed_reads_filename = input_reads_filename + GZIPPED_READS_SUFFIX;
+
+ // Populate database with read names and convert the fastq
+ // input into fasta for faidx
+ import_reads(input_reads_filename, m_indexed_reads_filename);
+
+ // build faidx
+ int ret = fai_build(m_indexed_reads_filename.c_str());
+ if(ret != 0) {
+ fprintf(stderr, "Error running faidx_build on %s\n", m_indexed_reads_filename.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ m_fai = NULL;
+}
+
+//
+void ReadDB::load(const std::string& input_reads_filename)
+{
+ // generate input filenames
+ m_indexed_reads_filename = input_reads_filename + GZIPPED_READS_SUFFIX;
+ std::string in_filename = m_indexed_reads_filename + READ_DB_SUFFIX;
+
+ //
+ std::ifstream in_file(in_filename.c_str());
+ if(in_file.bad()) {
+ fprintf(stderr, "error: could not read db %s\n", in_filename.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ // read the database
+ std::string line;
+ while(getline(in_file, line)) {
+ std::vector<std::string> fields = split(line, '\t');
+
+ std::string name = "";
+ std::string path = "";
+ if(fields.size() == 2) {
+ name = fields[0];
+ path = fields[1];
+ m_data[name].signal_data_path = path;
+ }
+ }
+
+ // load faidx
+ m_fai = fai_load(m_indexed_reads_filename.c_str());
+}
+
+ReadDB::~ReadDB()
+{
+ if(m_fai != NULL) {
+ fai_destroy(m_fai);
+ }
+}
+
+//
+void ReadDB::import_reads(const std::string& input_filename, const std::string& out_fasta_filename)
+{
+ // Open readers
+ FILE* read_fp = fopen(input_filename.c_str(), "r");
+ if(read_fp == NULL) {
+ fprintf(stderr, "error: could not open %s for read\n", input_filename.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ gzFile gz_read_fp = gzdopen(fileno(read_fp), "r");
+ if(gz_read_fp == NULL) {
+ fprintf(stderr, "error: could not open %s using gzdopen\n", input_filename.c_str());
+ exit(EXIT_FAILURE);
+ }
+
+ // Open writers
+ FILE* write_fp = fopen(out_fasta_filename.c_str(), "w");
+ BGZF* bgzf_write_fp = bgzf_dopen(fileno(write_fp), "w");
+
+ // read input sequences, add to DB and convert to fasta
+ kseq_t* seq = kseq_init(gz_read_fp);
+ while(kseq_read(seq) >= 0) {
+ // Check for a path to the fast5 file in the comment of the read
+ std::string path = "";
+ if(seq->comment.l > 0) {
+
+ // This splitting code implicitly handles both the 2 and 3 field
+ // fasta format that poretools will output. The FAST5 path
+ // is always the last field.
+ std::vector<std::string> fields = split(seq->comment.s, ' ');
+ path = fields.back();
+
+ // as a sanity check we require the path name to end in ".fast5"
+ if(path.substr(path.length() - 6) != ".fast5") {
+ path = "";
+ }
+ }
+
+ // sanity check that the read does not exist in the database
+ auto iter = m_data.find(seq->name.s);
+ if(iter != m_data.end()) {
+ fprintf(stderr, "Error: duplicate read name %s found in fasta file\n", seq->name.s);
+ exit(EXIT_FAILURE);
+ }
+
+ // add path
+ add_signal_path(seq->name.s, path);
+
+ // write sequence in gzipped fasta for fai indexing later
+ std::string out_record;
+ out_record += ">";
+ out_record += seq->name.s;
+ out_record += "\n";
+ out_record += seq->seq.s;
+ out_record += "\n";
+ bgzf_write(bgzf_write_fp, out_record.c_str(), out_record.length());
+ }
+
+ // cleanup
+ kseq_destroy(seq);
+
+ gzclose(gz_read_fp);
+ fclose(read_fp);
+
+ bgzf_close(bgzf_write_fp);
+ fclose(write_fp);
+}
+
+//
+void ReadDB::add_signal_path(const std::string& read_id, const std::string& path)
+{
+ m_data[read_id].signal_data_path = path;
+}
+
+//
+std::string ReadDB::get_signal_path(const std::string& read_id) const
+{
+ const auto& iter = m_data.find(read_id);
+ if(iter == m_data.end()) {
+ return "";
+ } else {
+ return iter->second.signal_data_path;
+ }
+}
+
+//
+std::string ReadDB::get_read_sequence(const std::string& read_id) const
+{
+ assert(m_fai != NULL);
+
+ int length;
+ char* seq;
+
+ // this call is not threadsafe
+ #pragma omp critical
+ seq = fai_fetch(m_fai, read_id.c_str(), &length);
+
+ if(seq == NULL) {
+ return "";
+ }
+
+ std::string out(seq);
+ free(seq);
+ return out;
+}
+
+//
+void ReadDB::save() const
+{
+ std::string out_filename = m_indexed_reads_filename + READ_DB_SUFFIX;
+
+ std::ofstream out_file(out_filename.c_str());
+
+ for(const auto& iter : m_data) {
+ const ReadDBData& entry = iter.second;
+ out_file << iter.first << "\t" << entry.signal_data_path << "\n";
+ }
+}
+
+
+
+//
+bool ReadDB::check_signal_paths() const
+{
+ for(const auto& iter : m_data) {
+ if(iter.second.signal_data_path == "") {
+ return false;
+ }
+ }
+ return true;
+}
+
+//
+void ReadDB::print_stats() const
+{
+ size_t num_reads_with_path = 0;
+ for(const auto& iter : m_data) {
+ num_reads_with_path += iter.second.signal_data_path != "";
+ }
+ fprintf(stderr, "[readdb] num reads: %zu, num reads with path: %zu\n", m_data.size(), num_reads_with_path);
+}
diff --git a/src/nanopolish_read_db.h b/src/nanopolish_read_db.h
new file mode 100644
index 0000000..6ce44c6
--- /dev/null
+++ b/src/nanopolish_read_db.h
@@ -0,0 +1,81 @@
+//---------------------------------------------------------
+// Copyright 2017 Ontario Institute for Cancer Research
+// Written by Jared Simpson (jared.simpson at oicr.on.ca)
+//---------------------------------------------------------
+//
+// nanopolish_read_db -- database of reads and their
+// associated signal data
+//
+#ifndef NANOPOLISH_READ_DB
+#define NANOPOLISH_READ_DB
+
+#include <map>
+#include "htslib/faidx.h"
+
+struct ReadDBData
+{
+ // path to the signal-level data for this read
+ std::string signal_data_path;
+};
+
+class ReadDB
+{
+ public:
+ ReadDB();
+ ~ReadDB();
+
+ //
+ // I/O
+ //
+
+ // construct the database from an input reads file
+ void build(const std::string& reads_filename);
+
+ // save the database to disk
+ void save() const;
+
+ // restore the database from disk
+ void load(const std::string& reads_filename);
+
+ //
+ // Data Access
+ //
+
+ // set the signal path for the given read
+ void add_signal_path(const std::string& read_id, const std::string& path);
+
+ // returns the path to the signal data for the given read
+ std::string get_signal_path(const std::string& read_id) const;
+
+ // returns the basecalled sequence for the given read
+ std::string get_read_sequence(const std::string& read_id) const;
+
+ // returns the number of reads in the database
+ size_t get_num_reads() const { return m_data.size(); }
+
+ //
+ // Summaries and sanity checks
+ //
+
+ // returns true if all reads in the database have paths to their signal-level data
+ bool check_signal_paths() const;
+
+ // print some summary stats about the database
+ void print_stats() const;
+
+ private:
+
+ //
+ void import_reads(const std::string& input_filename, const std::string& output_fasta_filename);
+
+ // the filename of the indexed data, after converting to fasta
+ std::string m_indexed_reads_filename;
+
+ //
+ std::map<std::string, ReadDBData> m_data;
+
+ //
+ faidx_t* m_fai;
+};
+
+#endif
diff --git a/src/nanopolish_scorereads.cpp b/src/nanopolish_scorereads.cpp
index 8a37326..64c419c 100644
--- a/src/nanopolish_scorereads.cpp
+++ b/src/nanopolish_scorereads.cpp
@@ -32,7 +32,7 @@
#include "nanopolish_matrix.h"
#include "nanopolish_profile_hmm.h"
#include "nanopolish_anchor.h"
-#include "nanopolish_fast5_map.h"
+#include "nanopolish_read_db.h"
#include "nanopolish_pore_model_set.h"
#include "H5pubconf.h"
@@ -427,7 +427,8 @@ int scorereads_main(int argc, char** argv)
parse_scorereads_options(argc, argv);
omp_set_num_threads(opt::num_threads);
- Fast5Map name_map(opt::reads_file);
+ ReadDB read_db;
+ read_db.load(opt::reads_file);
// Open the BAM and iterate over reads
@@ -515,8 +516,7 @@ int scorereads_main(int argc, char** argv)
//load read
std::string read_name = bam_get_qname(record);
- std::string fast5_path = name_map.get_path(read_name);
- SquiggleRead sr(read_name, fast5_path);
+ SquiggleRead sr(read_name, read_db);
// TODO: early exit when have processed all of the reads in readnames
if (!opt::readnames.empty() &&
diff --git a/src/nanopolish_squiggle_read.cpp b/src/nanopolish_squiggle_read.cpp
index 8171c12..7f9583f 100644
--- a/src/nanopolish_squiggle_read.cpp
+++ b/src/nanopolish_squiggle_read.cpp
@@ -12,6 +12,13 @@
#include "nanopolish_pore_model_set.h"
#include "nanopolish_methyltrain.h"
#include "nanopolish_extract.h"
+#include "nanopolish_raw_loader.h"
+
+extern "C" {
+#include "event_detection.h"
+#include "scrappie_common.h"
+}
+
#include <fast5.hpp>
//#define DEBUG_MODEL_SELECTION 1
@@ -24,27 +31,36 @@ int g_unparseable_reads = 0;
int g_qc_fail_reads = 0;
int g_failed_calibration_reads = 0;
+const double MIN_CALIBRATION_VAR = 2.5;
+
//
-SquiggleRead::SquiggleRead(const std::string& name, const std::string& path, const uint32_t flags) :
+SquiggleRead::SquiggleRead(const std::string& name, const ReadDB& read_db, const uint32_t flags) :
read_name(name),
pore_type(PT_UNKNOWN),
- fast5_path(path),
drift_correction_performed(false),
f_p(nullptr)
{
- events_per_base[0] = events_per_base[1] = 0.0f;
+ this->events_per_base[0] = events_per_base[1] = 0.0f;
+ this->fast5_path = read_db.get_signal_path(this->read_name);
#pragma omp critical(sr_load_fast5)
{
- load_from_fast5(flags);
+ bool is_event_read = is_extract_read_name(this->read_name);
+ if(is_event_read) {
+ load_from_events(flags);
+ } else {
+ this->read_sequence = read_db.get_read_sequence(read_name);
+ load_from_raw(flags);
+ }
+
+ // perform drift correction and other scalings
+ transform();
}
-
- // perform drift correction and other scalings
- transform();
}
SquiggleRead::~SquiggleRead()
{
+
}
// helper for get_closest_event_to
@@ -95,7 +111,7 @@ void SquiggleRead::transform()
}
//
-void SquiggleRead::load_from_fast5(const uint32_t flags)
+void SquiggleRead::load_from_events(const uint32_t flags)
{
f_p = new fast5::File(fast5_path);
assert(f_p->is_open());
@@ -204,6 +220,181 @@ void SquiggleRead::load_from_fast5(const uint32_t flags)
f_p = nullptr;
}
+//
+void SquiggleRead::load_from_raw(const uint32_t flags)
+{
+ // File not in db, can't load
+ if(this->fast5_path == "" || this->read_sequence == "") {
+ return;
+ }
+
+ // Hardcoded parameters, for now we can only do template with the main R9.4 model
+ size_t strand_idx = 0;
+ std::string kit = "r9.4_450bps";
+ std::string alphabet = "nucleotide";
+ std::string strand_str = "template";
+ size_t k = 6;
+
+ // Open file for read
+ this->f_p = new fast5::File(fast5_path);
+ assert(f_p->is_open());
+
+ // Read the sample rate
+ auto channel_params = f_p->get_channel_id_params();
+ this->sample_rate = channel_params.sampling_rate;
+
+ // Read the actual samples
+ auto& sample_read_names = f_p->get_raw_samples_read_name_list();
+ if(sample_read_names.empty()) {
+ fprintf(stderr, "Error, no raw samples found\n");
+ exit(EXIT_FAILURE);
+ }
+
+ // we assume the first raw sample read is the one we're after
+ std::string sample_read_name = sample_read_names.front();
+ std::vector<float> samples = f_p->get_raw_samples(sample_read_name);
+
+ // convert samples to scrappie's format (for event detection)
+ raw_table rt;
+ rt.n = samples.size();
+ rt.start = 0;
+ rt.end = samples.size() - 1;
+ rt.raw = (float*)malloc(sizeof(float) * samples.size());
+ for(size_t i = 0; i < samples.size(); ++i) {
+ rt.raw[i] = samples[i];
+ }
+
+ // trim using scrappie's internal method
+ // parameters taken directly from scrappie defaults
+ int trim_start = 200;
+ int trim_end = 10;
+ int varseg_chunk = 100;
+ float varseg_thresh = 0.0;
+ trim_and_segment_raw(rt, trim_start, trim_end, varseg_chunk, varseg_thresh);
+ event_table et = detect_events(rt, event_detection_defaults);
+ assert(rt.n > 0);
+ assert(et.n > 0);
+
+ // Load pore model and scale to events using method-of-moments
+ this->pore_model[strand_idx] = PoreModelSet::get_model(kit,
+ alphabet,
+ strand_str,
+ 6);
+ double shift, scale;
+ estimate_scalings_using_mom(this->read_sequence,
+ this->pore_model[strand_idx],
+ et,
+ shift,
+ scale);
+
+ // apply parameters to pore model
+ this->pore_model[strand_idx].shift = shift;
+ this->pore_model[strand_idx].scale = scale;
+ this->pore_model[strand_idx].drift = 0.0f;
+ this->pore_model[strand_idx].var = 1.0f;
+ transform();
+ this->pore_model[strand_idx].bake_gaussian_parameters();
+
+ // copy events into nanopolish's format
+ this->events[strand_idx].resize(et.n);
+ double start_time = 0;
+ for(size_t i = 0; i < et.n; ++i) {
+ float length_in_seconds = et.event[i].length / this->sample_rate;
+ this->events[strand_idx][i] = { et.event[i].mean, et.event[i].stdv, start_time, length_in_seconds, logf(et.event[i].stdv) };
+ start_time += length_in_seconds;
+ }
+
+ // clean up scrappie raw and event tables
+ assert(rt.raw != NULL);
+ assert(et.event != NULL);
+ free(rt.raw);
+ free(et.event);
+
+ // align events to the basecalled read
+ std::vector<AlignedPair> event_alignment = banded_simple_event_align(*this, read_sequence);
+
+ // transform alignment into the base-to-event map
+ if(event_alignment.size() > 0) {
+
+ // create base-to-event map
+ size_t n_kmers = read_sequence.size() - k + 1;
+ this->base_to_event_map.clear();
+ this->base_to_event_map.resize(n_kmers);
+
+ size_t max_event = 0;
+ size_t min_event = std::numeric_limits<size_t>::max();
+
+ size_t prev_event_idx = -1;
+ for(size_t i = 0; i < event_alignment.size(); ++i) {
+
+ size_t k_idx = event_alignment[i].ref_pos;
+ size_t event_idx = event_alignment[i].read_pos;
+ IndexPair& elem = this->base_to_event_map[k_idx].indices[strand_idx];
+ if(event_idx != prev_event_idx) {
+ if(elem.start == -1) {
+ elem.start = event_idx;
+ }
+ elem.stop = event_idx;
+ }
+
+ max_event = std::max(max_event, event_idx);
+ min_event = std::min(min_event, event_idx);
+ prev_event_idx = event_idx;
+ }
+
+ events_per_base[strand_idx] = (double)(max_event - min_event) / n_kmers;
+
+ // prepare data structures for the final calibration
+ std::vector<EventAlignment> alignment =
+ get_eventalignment_for_1d_basecalls(read_sequence, this->base_to_event_map, k, strand_idx, 0);
+
+ // reset default scaling parameters
+ this->pore_model[strand_idx].shift = 0.0;
+ this->pore_model[strand_idx].scale = 1.0;
+ this->pore_model[strand_idx].drift = 0.0;
+ this->pore_model[strand_idx].var = 1.0;
+ this->pore_model[strand_idx].scale_sd = 1.0;
+ this->pore_model[strand_idx].var_sd = 1.0;
+ this->pore_model[strand_idx].bake_gaussian_parameters();
+
+ // run recalibration to get the best set of scaling parameters and the residual
+ // between the (scaled) event levels and the model.
+ // internally this function will set shift/scale/etc of the pore model
+ bool calibrated = recalibrate_model(*this, strand_idx, alignment, this->pore_model[strand_idx].pmalphabet, true, false);
+
+#ifdef DEBUG_MODEL_SELECTION
+ fprintf(stderr, "[calibration] read: %s"
+ "scale: %.2lf shift: %.2lf drift: %.5lf var: %.2lf\n",
+ read_name.substr(0, 6).c_str(), pore_model[strand_idx].scale,
+ pore_model[strand_idx].shift, pore_model[strand_idx].drift, pore_model[strand_idx].var);
+#endif
+
+ // QC calibration
+ if(!calibrated || this->pore_model[strand_idx].var > MIN_CALIBRATION_VAR) {
+ events[strand_idx].clear();
+ g_failed_calibration_reads += 1;
+ }
+ } else {
+ // Could not align, fail this read
+ this->events[strand_idx].clear();
+ this->events_per_base[strand_idx] = 0.0f;
+ g_failed_calibration_reads += 1;
+ }
+ g_total_reads += 1;
+
+ // Filter poor quality reads that have too many "stays"
+ if(!this->events[strand_idx].empty() && this->events_per_base[strand_idx] > 5.0) {
+ g_qc_fail_reads += 1;
+ events[0].clear();
+ events[1].clear();
+ }
+
+ delete f_p;
+ f_p = nullptr;
+}
+
+
+
void SquiggleRead::_load_R7(uint32_t si)
{
assert(f_p and f_p->is_open());
@@ -404,7 +595,7 @@ void SquiggleRead::_load_R9(uint32_t si,
#endif
}
- if(best_model_var < 2.5) {
+ if(best_model_var < MIN_CALIBRATION_VAR) {
#ifdef DEBUG_MODEL_SELECTION
fprintf(stderr, "[calibration] selected model with var %.4lf\n", best_model_var);
#endif
@@ -976,6 +1167,15 @@ bool SquiggleRead::check_basecall_group() const
return true;
}
+//
+bool SquiggleRead::is_extract_read_name(std::string& name) const
+{
+ // albacore read names are uuids with hex characters separated
+ // by underscores. If the read name contains three colon-delimited fields
+ // we infer it is output by nanopolish extract.
+ return std::count(name.begin(), name.end(), ':') == 2;
+}
+
void SquiggleRead::detect_basecall_group()
{
assert(f_p and f_p->is_open());
diff --git a/src/nanopolish_squiggle_read.h b/src/nanopolish_squiggle_read.h
index 1bd8176..7a3666a 100644
--- a/src/nanopolish_squiggle_read.h
+++ b/src/nanopolish_squiggle_read.h
@@ -13,6 +13,7 @@
#include "nanopolish_poremodel.h"
#include "nanopolish_transition_parameters.h"
#include "nanopolish_eventalign.h"
+#include "nanopolish_read_db.h"
#include <string>
enum PoreType
@@ -41,11 +42,11 @@ enum SquiggleReadFlags
// The raw event data for a read
struct SquiggleEvent
{
- float mean; // current level mean in picoamps
- float stdv; // current level stdv
+ float mean; // current level mean in picoamps
+ float stdv; // current level stdv
double start_time; // start time of the event in seconds
- float duration; // duration of the event in seconds
- float log_stdv; // precompute for efficiency
+ float duration; // duration of the event in seconds
+ float log_stdv; // precompute for efficiency
};
struct IndexPair
@@ -68,7 +69,7 @@ class SquiggleRead
public:
SquiggleRead() : drift_correction_performed(false) {} // legacy TODO remove
- SquiggleRead(const std::string& name, const std::string& path, const uint32_t flags = 0);
+ SquiggleRead(const std::string& name, const ReadDB& read_db, const uint32_t flags = 0);
~SquiggleRead();
//
@@ -195,7 +196,7 @@ class SquiggleRead
// one event sequence for each strand
std::vector<SquiggleEvent> events[2];
-
+
// optional fields holding the raw data
// this is not split into strands so there is only one vector, unlike events
std::vector<float> samples;
@@ -218,8 +219,11 @@ class SquiggleRead
SquiggleRead(const SquiggleRead&) {}
- // Load all the read data from a fast5 file
- void load_from_fast5(const uint32_t flags);
+ // Load all read data from events in a fast5 file
+ void load_from_events(const uint32_t flags);
+
+ // Load all read data from raw samples
+ void load_from_raw(const uint32_t flags);
// Version-specific intialization functions
void _load_R7(uint32_t si);
@@ -231,13 +235,13 @@ class SquiggleRead
// make a map from a base of the 1D read sequence to the range of events supporting that base
std::vector<EventRangeForBase> build_event_map_1d(const std::string& read_sequence_1d,
- uint32_t strand,
+ uint32_t strand,
std::vector<fast5::Basecall_Event>& f5_events);
-
+
std::vector<EventRangeForBase> read_reconstruction(const std::string& read_sequence_1d,
- uint32_t strand,
+ uint32_t strand,
std::vector<fast5::Basecall_Event>& f5_events);
-
+
void _find_kmer_event_pair(const std::string& read_sequence_1d,
std::vector<fast5::Basecall_Event>& events,
size_t& k_idx,
@@ -249,10 +253,16 @@ class SquiggleRead
// helper for get_closest_event_to
int get_next_event(int start, int stop, int stride, uint32_t strand) const;
+
// detect pore_type
void detect_pore_type();
+
+ // check whether the input read name conforms to nanopolish extract's signature
+ bool is_extract_read_name(std::string& name) const;
+
// detect basecall_group and read_type
void detect_basecall_group();
+
// check basecall_group and read_type
bool check_basecall_group() const;
};
diff --git a/src/nanopolish_train_poremodel_from_basecalls.cpp b/src/nanopolish_train_poremodel_from_basecalls.cpp
index 1cdf997..ce54cea 100644
--- a/src/nanopolish_train_poremodel_from_basecalls.cpp
+++ b/src/nanopolish_train_poremodel_from_basecalls.cpp
@@ -216,7 +216,7 @@ void alignment_to_training_data(const SquiggleRead* read,
int train_poremodel_from_basecalls_main(int argc, char** argv)
{
parse_train_poremodel_from_basecalls_options(argc, argv);
-
+#if 0
std::ifstream fofn_reader(opt::fofn_file);
std::string fast5_name;
@@ -397,6 +397,6 @@ int train_poremodel_from_basecalls_main(int argc, char** argv)
for(auto* read : reads) {
delete read;
}
-
+#endif
return 0;
}
diff --git a/src/test/nanopolish_test.cpp b/src/test/nanopolish_test.cpp
index b34ffd2..d17e2da 100644
--- a/src/test/nanopolish_test.cpp
+++ b/src/test/nanopolish_test.cpp
@@ -337,9 +337,14 @@ std::string event_alignment_to_string(const std::vector<HMMAlignmentState>& alig
TEST_CASE( "hmm", "[hmm]") {
- // read the FAST5
- SquiggleRead sr("01234567-0123-0123-0123-0123456789ab:2D_000:2d", "test/data/LomanLabz_PC_Ecoli_K12_R7.3_2549_1_ch8_file30_strand.fast5");
- sr.transform();
+ // load the FAST5
+ std::string read_id = "01234567-0123-0123-0123-0123456789ab:2D_000:2d";
+ std::string fast5_path = "test/data/LomanLabz_PC_Ecoli_K12_R7.3_2549_1_ch8_file30_strand.fast5";
+
+ ReadDB read_db;
+ read_db.add_signal_path(read_id, fast5_path);
+
+ SquiggleRead sr(read_id, read_db);
// The reference sequence to align to:
std::string ref_subseq = "ATCAGTAAAATAACGTAGAGCGGTAACCTTGCCATAAAGGTCGAGTTTA"
@@ -372,11 +377,11 @@ TEST_CASE( "hmm", "[hmm]") {
"MMMKMMMMMKMKMKMEMKKMKMKKMMMMMMEMMMMKMKMEEMMMMKMEEEEEM";
expected_alignment[1] =
- "MMKMMMKMEEMMKMKMKMEMMMKMMMKMEMMMKMMMKMMMMMMMMMKKMEMMMM"
- "EMMMMMMMMKMKKMMMMMMMEMMMMMKMMMMMKMEMMMMMKMMMMMEEEEEEEEM";
+ "MMKMMMKMEEMMKMKMKMEMMMKMMMKMEMMMKMMMKMMMMMMMMMKKMEMMM"
+ "EMMMMMMMMMKMKKMMMMMMMEMMMMMKMMMMMKMEMMMMMKMMMMMEEEEEEEEM";
- double expected_viterbi_last_state[2] = { -237.7690734863, -266.2348022461 };
- double expected_forward[2] = { -216.053604126, -254.2341003418 };
+ double expected_viterbi_last_state[2] = { -237.7808380127, -267.9027709961 };
+ double expected_forward[2] = { -216.053604126, -254.5881347656 };
for(int si = 0; si <= 1; ++si) {
diff --git a/src/thirdparty/scrappie/event_detection.c b/src/thirdparty/scrappie/event_detection.c
new file mode 100644
index 0000000..a228498
--- /dev/null
+++ b/src/thirdparty/scrappie/event_detection.c
@@ -0,0 +1,319 @@
+#include <float.h>
+#include <math.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <stdio.h>
+
+#include "event_detection.h"
+#include "scrappie_stdlib.h"
+
+typedef struct {
+ int DEF_PEAK_POS;
+ float DEF_PEAK_VAL;
+ float *signal;
+ size_t signal_length;
+ float threshold;
+ size_t window_length;
+ size_t masked_to;
+ int peak_pos;
+ float peak_value;
+ bool valid_peak;
+} Detector;
+typedef Detector *DetectorPtr;
+
+/**
+ * Compute cumulative sum and sum of squares for a vector of data
+ *
+ * Element i sum (sumsq) is the sum (sum of squares) up to but
+ * excluding element i of the inputy data.
+ *
+ * @param data double[d_length] Data to be summed over (in)
+ * @param sum double[d_length + 1] Vector to store sum (out)
+ * @param sumsq double[d_length + 1] Vector to store sum of squares (out)
+ * @param d_length Length of data vector
+ **/
+void compute_sum_sumsq(const float *data, double *sum,
+ double *sumsq, size_t d_length) {
+ RETURN_NULL_IF(NULL == data, );
+ RETURN_NULL_IF(NULL == sum, );
+ RETURN_NULL_IF(NULL == sumsq, );
+ assert(d_length > 0);
+
+ sum[0] = 0.0f;
+ sumsq[0] = 0.0f;
+ for (size_t i = 0; i < d_length; ++i) {
+ sum[i + 1] = sum[i] + data[i];
+ sumsq[i + 1] = sumsq[i] + data[i] * data[i];
+ }
+}
+
+/**
+ * Compute windowed t-statistic from summary information
+ *
+ * @param sum double[d_length] Cumulative sums of data (in)
+ * @param sumsq double[d_length] Cumulative sum of squares of data (in)
+ * @param d_length Length of data vector
+ * @param w_length Window length to calculate t-statistic over
+ *
+ * @returns float array containing tstats. Returns NULL on error
+ **/
+float *compute_tstat(const double *sum, const double *sumsq,
+ size_t d_length, size_t w_length) {
+ assert(d_length > 0);
+ assert(w_length > 0);
+ RETURN_NULL_IF(NULL == sum, NULL);
+ RETURN_NULL_IF(NULL == sumsq, NULL);
+
+ float *tstat = calloc(d_length, sizeof(float));
+ RETURN_NULL_IF(NULL == tstat, NULL);
+
+ const float eta = FLT_MIN;
+ const float w_lengthf = (float)w_length;
+
+ // Quick return:
+ // t-test not defined for number of points less than 2
+ // need at least as many points as twice the window length
+ if (d_length < 2 * w_length || w_length < 2) {
+ for (size_t i = 0; i < d_length; ++i) {
+ tstat[i] = 0.0f;
+ }
+ return tstat;
+ }
+ // fudge boundaries
+ for (size_t i = 0; i < w_length; ++i) {
+ tstat[i] = 0;
+ tstat[d_length - i - 1] = 0;
+ }
+
+ // get to work on the rest
+ for (size_t i = w_length; i <= d_length - w_length; ++i) {
+ double sum1 = sum[i];
+ double sumsq1 = sumsq[i];
+ if (i > w_length) {
+ sum1 -= sum[i - w_length];
+ sumsq1 -= sumsq[i - w_length];
+ }
+ float sum2 = (float)(sum[i + w_length] - sum[i]);
+ float sumsq2 = (float)(sumsq[i + w_length] - sumsq[i]);
+ float mean1 = sum1 / w_lengthf;
+ float mean2 = sum2 / w_lengthf;
+ float combined_var = sumsq1 / w_lengthf - mean1 * mean1
+ + sumsq2 / w_lengthf - mean2 * mean2;
+
+ // Prevent problem due to very small variances
+ combined_var = fmaxf(combined_var, eta);
+
+ //t-stat
+ // Formula is a simplified version of Student's t-statistic for the
+ // special case where there are two samples of equal size with
+ // differing variance
+ const float delta_mean = mean2 - mean1;
+ tstat[i] = fabs(delta_mean) / sqrt(combined_var / w_lengthf);
+ }
+
+ return tstat;
+}
+
+/**
+ *
+ * @returns array of length nsample whose elements contain peak positions
+ * Remaining elements are padded by zeros.
+ **/
+size_t *short_long_peak_detector(DetectorPtr short_detector,
+ DetectorPtr long_detector,
+ const float peak_height) {
+ assert(short_detector->signal_length == long_detector->signal_length);
+ RETURN_NULL_IF(NULL == short_detector->signal, NULL);
+ RETURN_NULL_IF(NULL == long_detector->signal, NULL);
+
+ const size_t ndetector = 2;
+ DetectorPtr detectors[] = { short_detector, long_detector };
+
+ size_t *peaks = calloc(short_detector->signal_length, sizeof(size_t));
+ RETURN_NULL_IF(NULL == peaks, NULL);
+
+ size_t peak_count = 0;
+ for (size_t i = 0; i < short_detector->signal_length; i++) {
+ for (int k = 0; k < ndetector; k++) {
+ DetectorPtr detector = detectors[k];
+ //Carry on if we've been masked out
+ if (detector->masked_to >= i) {
+ continue;
+ }
+
+ float current_value = detector->signal[i];
+
+ if (detector->peak_pos == detector->DEF_PEAK_POS) {
+ //CASE 1: We've not yet recorded a maximum
+ if (current_value < detector->peak_value) {
+ //Either record a deeper minimum...
+ detector->peak_value = current_value;
+ } else if (current_value - detector->peak_value >
+ peak_height) {
+ // ...or we've seen a qualifying maximum
+ detector->peak_value = current_value;
+ detector->peak_pos = i;
+ //otherwise, wait to rise high enough to be considered a peak
+ }
+ } else {
+ //CASE 2: In an existing peak, waiting to see if it is good
+ if (current_value > detector->peak_value) {
+ //Update the peak
+ detector->peak_value = current_value;
+ detector->peak_pos = i;
+ }
+ //Dominate other tstat signals if we're going to fire at some point
+ if (detector == short_detector) {
+ if (detector->peak_value > detector->threshold) {
+ long_detector->masked_to =
+ detector->peak_pos + detector->window_length;
+ long_detector->peak_pos =
+ long_detector->DEF_PEAK_POS;
+ long_detector->peak_value =
+ long_detector->DEF_PEAK_VAL;
+ long_detector->valid_peak = false;
+ }
+ }
+ //Have we convinced ourselves we've seen a peak
+ if (detector->peak_value - current_value > peak_height
+ && detector->peak_value > detector->threshold) {
+ detector->valid_peak = true;
+ }
+ //Finally, check the distance if this is a good peak
+ if (detector->valid_peak
+ && (i - detector->peak_pos) >
+ detector->window_length / 2) {
+ //Emit the boundary and reset
+ peaks[peak_count] = detector->peak_pos;
+ peak_count++;
+ detector->peak_pos = detector->DEF_PEAK_POS;
+ detector->peak_value = current_value;
+ detector->valid_peak = false;
+ }
+ }
+ }
+ }
+
+ return peaks;
+}
+
+/** Create an event given boundaries
+ *
+ * Note: Bounds are CADLAG (i.e. lower bound is contained in the interval but
+ * the upper bound is not).
+ *
+ * @param start Index of lower bound
+ * @param end Index of upper bound
+ * @param sums
+ * @param sumsqs
+ * @param nsample Total number of samples in read
+ *
+ * @returns An initialised event. A 'null' event is returned on error.
+ **/
+event_t create_event(size_t start, size_t end, double const *sums,
+ double const *sumsqs, size_t nsample) {
+ assert(start < nsample);
+ assert(end <= nsample);
+
+ event_t event = { 0 };
+ event.pos = -1;
+ event.state = -1;
+ RETURN_NULL_IF(NULL == sums, event);
+ RETURN_NULL_IF(NULL == sumsqs, event);
+
+ event.start = (uint64_t)start;
+ event.length = (float)(end - start);
+ event.mean = (float)(sums[end] - sums[start]) / event.length;
+ const float deltasqr = (sumsqs[end] - sumsqs[start]);
+ const float var = deltasqr / event.length - event.mean * event.mean;
+ event.stdv = sqrtf(fmaxf(var, 0.0f));
+
+ return event;
+}
+
+event_table create_events(size_t const *peaks, double const *sums,
+ double const *sumsqs, size_t nsample) {
+ event_table et = { 0 };
+ RETURN_NULL_IF(NULL == sums, et);
+ RETURN_NULL_IF(NULL == sumsqs, et);
+ RETURN_NULL_IF(NULL == peaks, et);
+
+ // Count number of events found
+ size_t n = 1;
+ for (size_t i = 0; i < nsample; ++i) {
+ if (peaks[i] > 0 && peaks[i] < nsample) {
+ n++;
+ }
+ }
+
+ et.event = calloc(n, sizeof(event_t));
+ RETURN_NULL_IF(NULL == et.event, et);
+
+ et.n = n;
+ et.end = et.n;
+
+
+ // First event -- starts at zero
+ et.event[0] = create_event(0, peaks[0], sums, sumsqs, nsample);
+ // Other events -- peak[i-1] -> peak[i]
+ for(size_t ev=1 ; ev < n - 1 ; ev++){
+ et.event[ev] = create_event(peaks[ev - 1], peaks[ev], sums, sumsqs, nsample);
+ }
+ // Last event -- ends at nsample
+ et.event[n - 1] = create_event(peaks[n - 2], nsample, sums, sumsqs, nsample);
+
+ return et;
+}
+
+event_table detect_events(raw_table const rt, detector_param const edparam) {
+
+ event_table et = { 0 };
+ RETURN_NULL_IF(NULL == rt.raw, et);
+
+ double *sums = calloc(rt.n + 1, sizeof(double));
+ double *sumsqs = calloc(rt.n + 1, sizeof(double));
+
+ compute_sum_sumsq(rt.raw, sums, sumsqs, rt.n);
+ float *tstat1 = compute_tstat(sums, sumsqs, rt.n, edparam.window_length1);
+ float *tstat2 = compute_tstat(sums, sumsqs, rt.n, edparam.window_length2);
+
+ Detector short_detector = {
+ .DEF_PEAK_POS = -1,
+ .DEF_PEAK_VAL = FLT_MAX,
+ .signal = tstat1,
+ .signal_length = rt.n,
+ .threshold = edparam.threshold1,
+ .window_length = edparam.window_length1,
+ .masked_to = 0,
+ .peak_pos = -1,
+ .peak_value = FLT_MAX,
+ .valid_peak = false
+ };
+
+ Detector long_detector = {
+ .DEF_PEAK_POS = -1,
+ .DEF_PEAK_VAL = FLT_MAX,
+ .signal = tstat2,
+ .signal_length = rt.n,
+ .threshold = edparam.threshold2,
+ .window_length = edparam.window_length2,
+ .masked_to = 0,
+ .peak_pos = -1,
+ .peak_value = FLT_MAX,
+ .valid_peak = false
+ };
+
+ size_t *peaks =
+ short_long_peak_detector(&short_detector, &long_detector,
+ edparam.peak_height);
+
+ et = create_events(peaks, sums, sumsqs, rt.n);
+
+ free(peaks);
+ free(tstat2);
+ free(tstat1);
+ free(sumsqs);
+ free(sums);
+
+ return et;
+}
diff --git a/src/thirdparty/scrappie/event_detection.h b/src/thirdparty/scrappie/event_detection.h
new file mode 100644
index 0000000..87f9ac5
--- /dev/null
+++ b/src/thirdparty/scrappie/event_detection.h
@@ -0,0 +1,27 @@
+#ifndef EVENT_DETECTION_H
+# define EVENT_DETECTION_H
+
+# include "scrappie_structures.h"
+
+typedef struct {
+ size_t window_length1;
+ size_t window_length2;
+ float threshold1;
+ float threshold2;
+ float peak_height;
+} detector_param;
+
+
+static detector_param const event_detection_defaults = {
+ .window_length1 = 3,
+ .window_length2 = 6,
+ .threshold1 = 1.4f,
+ .threshold2 = 9.0f,
+ .peak_height = 0.2f
+};
+
+
+
+event_table detect_events(raw_table const rt, detector_param const edparam);
+
+#endif /* EVENT_DETECTION_H */
diff --git a/src/thirdparty/scrappie/scrappie_common.c b/src/thirdparty/scrappie/scrappie_common.c
new file mode 100644
index 0000000..4b148e3
--- /dev/null
+++ b/src/thirdparty/scrappie/scrappie_common.c
@@ -0,0 +1,73 @@
+#include "scrappie_common.h"
+#include "scrappie_stdlib.h"
+#include "util.h"
+
+raw_table trim_and_segment_raw(raw_table rt, int trim_start, int trim_end, int varseg_chunk, float varseg_thresh) {
+ RETURN_NULL_IF(NULL == rt.raw, (raw_table){0});
+
+ rt = trim_raw_by_mad(rt, varseg_chunk, varseg_thresh);
+ RETURN_NULL_IF(NULL == rt.raw, (raw_table){0});
+
+ rt.start += trim_start;
+ rt.end -= trim_end;
+
+ if (rt.start >= rt.end) {
+ free(rt.raw);
+ return (raw_table){0};
+ }
+
+ return rt;
+}
+
+/** Simple segmentation of a raw read by thresholding the MAD
+ *
+ * The MAD of the raw signal is calculated for non-overlapping chunks and then
+ * thresholded to find regions at the beginning and end of the signal that have
+ * unusually low variation (generally a stall or open pore). The threshhold is
+ * derived from the distribution of the calaculated MADs.
+ *
+ * The threshold is chosen to be high since a single chunk above it will trigger
+ * the end of the trimming: the threshhold is chosen so it is unlikely to be
+ * exceeded in the leader but commonly exceeded in the main read.
+ *
+ * @param rt Structure containing raw signal
+ * @param chunk_size Size of non-overlapping chunks
+ * @param perc The quantile to be calculated to use for threshholding
+ *
+ * @return A range structure containing new start and end for read
+ **/
+raw_table trim_raw_by_mad(raw_table rt, int chunk_size, float perc) {
+ assert(chunk_size > 1);
+ assert(perc >= 0.0 && perc <= 1.0);
+
+ const size_t nsample = rt.end - rt.start;
+ const size_t nchunk = nsample / chunk_size;
+ // Truncation of end to be consistent with Sloika
+ rt.end = nchunk * chunk_size;
+
+ float *madarr = malloc(nchunk * sizeof(float));
+ RETURN_NULL_IF(NULL == madarr, (raw_table){0});
+ for (size_t i = 0; i < nchunk; i++) {
+ madarr[i] = madf(rt.raw + rt.start + i * chunk_size, chunk_size, NULL);
+ }
+ quantilef(madarr, nchunk, &perc, 1);
+
+ const float thresh = perc;
+ for (size_t i = 0; i < nchunk; i++) {
+ if (madarr[i] > thresh) {
+ break;
+ }
+ rt.start += chunk_size;
+ }
+ for (size_t i = nchunk; i > 0; i--) {
+ if (madarr[i - 1] > thresh) {
+ break;
+ }
+ rt.end -= chunk_size;
+ }
+ assert(rt.end > rt.start);
+
+ free(madarr);
+
+ return rt;
+}
diff --git a/src/thirdparty/scrappie/scrappie_common.h b/src/thirdparty/scrappie/scrappie_common.h
new file mode 100644
index 0000000..c2ff301
--- /dev/null
+++ b/src/thirdparty/scrappie/scrappie_common.h
@@ -0,0 +1,10 @@
+#pragma once
+#ifndef SCRAPPIE_COMMON_H
+#define SCRAPPIE_COMMON_H
+
+#include "scrappie_structures.h"
+
+raw_table trim_and_segment_raw(raw_table rt, int trim_start, int trim_end, int varseg_chunk, float varseg_thresh);
+raw_table trim_raw_by_mad(raw_table rt, int chunk_size, float proportion);
+
+#endif /* SCRAPPIE_COMMON_H */
diff --git a/src/thirdparty/scrappie/scrappie_stdlib.h b/src/thirdparty/scrappie/scrappie_stdlib.h
new file mode 100644
index 0000000..7aff438
--- /dev/null
+++ b/src/thirdparty/scrappie/scrappie_stdlib.h
@@ -0,0 +1,41 @@
+#pragma once
+#ifndef SCRAPPIE_STDLIB
+# define SCRAPPIE_STDLIB
+# include <assert.h>
+# include <err.h>
+# include <stdlib.h>
+# include <string.h>
+
+# if defined(CHAOSMONKEY) && ! defined(BANANA)
+# define FUZZTHRESH (long int)(CHAOSMONKEY * RAND_MAX)
+# define malloc(A) (warnx("Opportunity for chaos at %s:%d", __FILE__, __LINE__), rand() > FUZZTHRESH) \
+ ? malloc(A) \
+ : (warnx("Injected NULL at %s:%d", __FILE__, __LINE__), NULL)
+
+# define calloc(A, B) (warnx("Opportunity for chaos at %s:%d", __FILE__, __LINE__), rand() > FUZZTHRESH) \
+ ? calloc(A, B) \
+ : (warnx("Injected NULL at %s:%d", __FILE__, __LINE__), NULL)
+
+# define aligned_alloc(A, B) (warnx("Opportunity for chaos at %s:%d", __FILE__, __LINE__), rand() > FUZZTHRESH) \
+ ? aligned_alloc(A, B) \
+ : (warnx("Injected NULL at %s:%d", __FILE__, __LINE__), NULL)
+# else
+# define malloc(A) malloc(A)
+# define calloc(A, B) calloc(A, B)
+# define aligned_alloc(A, B) aligned_alloc(A, B)
+# endif /* CHAOSMONKEY */
+
+# ifdef ABORT_ON_NULL
+# define RETURN_NULL_IF(A, B) \
+ if (A) { \
+ warnx("Failure at %s : %d", __FILE__, __LINE__); \
+ abort(); \
+ }
+# else
+# define RETURN_NULL_IF(A, B) if (A) { return B; }
+# endif
+
+
+// strlen on NULL is undefined. Define it.
+# define strlen(A) (NULL != A) ? strlen(A) : 0
+#endif /* SCRAPPIE_STDLIB */
diff --git a/src/thirdparty/scrappie/scrappie_structures.h b/src/thirdparty/scrappie/scrappie_structures.h
new file mode 100644
index 0000000..a836fb2
--- /dev/null
+++ b/src/thirdparty/scrappie/scrappie_structures.h
@@ -0,0 +1,31 @@
+#pragma once
+#ifndef SCRAPPIE_STRUCTURES_H
+# define SCRAPPIE_STRUCTURES_H
+
+# include <inttypes.h>
+# include <stddef.h>
+
+typedef struct {
+ uint64_t start;
+ float length;
+ float mean;
+ float stdv;
+ int pos;
+ int state;
+} event_t;
+
+typedef struct {
+ size_t n;
+ size_t start;
+ size_t end;
+ event_t *event;
+} event_table;
+
+typedef struct {
+ size_t n;
+ size_t start;
+ size_t end;
+ float *raw;
+} raw_table;
+
+#endif /* SCRAPPIE_DATA_H */
diff --git a/src/thirdparty/scrappie/sse_mathfun.h b/src/thirdparty/scrappie/sse_mathfun.h
new file mode 100644
index 0000000..c0c470a
--- /dev/null
+++ b/src/thirdparty/scrappie/sse_mathfun.h
@@ -0,0 +1,722 @@
+/* SIMD (SSE1+MMX or SSE2) implementation of sin, cos, exp and log
+
+ Inspired by Intel Approximate Math library, and based on the
+ corresponding algorithms of the cephes math library
+
+ The default is to use the SSE1 version. If you define USE_SSE2 the
+ the SSE2 intrinsics will be used in place of the MMX intrinsics. Do
+ not expect any significant performance improvement with SSE2.
+*/
+
+/* Copyright (C) 2007 Julien Pommier
+
+ This software is provided 'as-is', without any express or implied
+ warranty. In no event will the authors be held liable for any damages
+ arising from the use of this software.
+
+ Permission is granted to anyone to use this software for any purpose,
+ including commercial applications, and to alter it and redistribute it
+ freely, subject to the following restrictions:
+
+ 1. The origin of this software must not be misrepresented; you must not
+ claim that you wrote the original software. If you use this software
+ in a product, an acknowledgment in the product documentation would be
+ appreciated but is not required.
+ 2. Altered source versions must be plainly marked as such, and must not be
+ misrepresented as being the original software.
+ 3. This notice may not be removed or altered from any source distribution.
+
+ (this is the zlib license)
+*/
+
+/* Copyright (C) 2016 Oxford Nanopore technologies
+
+ This software contain modifications (C) Oxford Nanopore Technologies
+ and is not the original software, which can be obtained from
+ http://gruntthepeon.free.fr/ssemath/
+*/
+#pragma once
+#ifndef SSE_MATHFUN_H
+#define SSE_MATHFUN_H
+
+
+#include <xmmintrin.h>
+
+/* yes I know, the top of this file is quite ugly */
+
+#ifdef _MSC_VER /* visual c++ */
+# define ALIGN16_BEG __declspec(align(16))
+# define ALIGN16_END
+#else /* gcc or icc */
+# define ALIGN16_BEG
+# define ALIGN16_END __attribute__((aligned(16)))
+#endif
+
+/* __m128 is ugly to write */
+typedef __m128 v4sf; // vector of 4 float (sse1)
+
+#ifdef USE_SSE2
+# include <emmintrin.h>
+typedef __m128i v4si; // vector of 4 int (sse2)
+#else
+typedef __m64 v2si; // vector of 2 int (mmx)
+#endif
+
+/* declare some SSE constants -- why can't I figure a better way to do that? */
+#define _PS_CONST(Name, Val) \
+ static const ALIGN16_BEG float _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
+#define _PI32_CONST(Name, Val) \
+ static const ALIGN16_BEG int _pi32_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
+#define _PS_CONST_TYPE(Name, Type, Val) \
+ static const ALIGN16_BEG Type _ps_##Name[4] ALIGN16_END = { Val, Val, Val, Val }
+
+_PS_CONST(1 , 1.0f);
+_PS_CONST(0p5, 0.5f);
+/* the smallest non denormalized float number */
+_PS_CONST_TYPE(min_norm_pos, int, 0x00800000);
+_PS_CONST_TYPE(mant_mask, int, 0x7f800000);
+_PS_CONST_TYPE(inv_mant_mask, int, ~0x7f800000);
+
+_PS_CONST_TYPE(sign_mask, int, (int)0x80000000);
+_PS_CONST_TYPE(inv_sign_mask, int, ~0x80000000);
+
+_PI32_CONST(1, 1);
+_PI32_CONST(inv1, ~1);
+_PI32_CONST(2, 2);
+_PI32_CONST(4, 4);
+_PI32_CONST(0x7f, 0x7f);
+
+_PS_CONST(cephes_SQRTHF, 0.707106781186547524);
+_PS_CONST(cephes_log_p0, 7.0376836292E-2);
+_PS_CONST(cephes_log_p1, - 1.1514610310E-1);
+_PS_CONST(cephes_log_p2, 1.1676998740E-1);
+_PS_CONST(cephes_log_p3, - 1.2420140846E-1);
+_PS_CONST(cephes_log_p4, + 1.4249322787E-1);
+_PS_CONST(cephes_log_p5, - 1.6668057665E-1);
+_PS_CONST(cephes_log_p6, + 2.0000714765E-1);
+_PS_CONST(cephes_log_p7, - 2.4999993993E-1);
+_PS_CONST(cephes_log_p8, + 3.3333331174E-1);
+_PS_CONST(cephes_log_q1, -2.12194440e-4);
+_PS_CONST(cephes_log_q2, 0.693359375);
+
+#ifndef USE_SSE2
+typedef union xmm_mm_union {
+ __m128 xmm;
+ __m64 mm[2];
+} xmm_mm_union;
+
+#define COPY_XMM_TO_MM(xmm_, mm0_, mm1_) { \
+ xmm_mm_union u; u.xmm = xmm_; \
+ mm0_ = u.mm[0]; \
+ mm1_ = u.mm[1]; \
+}
+
+#define COPY_MM_TO_XMM(mm0_, mm1_, xmm_) { \
+ xmm_mm_union u; u.mm[0]=mm0_; u.mm[1]=mm1_; xmm_ = u.xmm; \
+ }
+
+#endif // USE_SSE2
+
+/* natural logarithm computed for 4 simultaneous float
+ return NaN for x <= 0
+*/
+static inline v4sf __attribute__((__always_inline__)) log_ps(v4sf x) {
+#ifdef USE_SSE2
+ v4si emm0;
+#else
+ v2si mm0, mm1;
+#endif
+ v4sf one = *(v4sf*)_ps_1;
+
+ v4sf invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
+
+ x = _mm_max_ps(x, *(v4sf*)_ps_min_norm_pos); /* cut off denormalized stuff */
+
+#ifndef USE_SSE2
+ /* part 1: x = frexpf(x, &e); */
+ COPY_XMM_TO_MM(x, mm0, mm1);
+ mm0 = _mm_srli_pi32(mm0, 23);
+ mm1 = _mm_srli_pi32(mm1, 23);
+#else
+ emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
+#endif
+ /* keep only the fractional part */
+ x = _mm_and_ps(x, *(v4sf*)_ps_inv_mant_mask);
+ x = _mm_or_ps(x, *(v4sf*)_ps_0p5);
+
+#ifndef USE_SSE2
+ /* now e=mm0:mm1 contain the really base-2 exponent */
+ mm0 = _mm_sub_pi32(mm0, *(v2si*)_pi32_0x7f);
+ mm1 = _mm_sub_pi32(mm1, *(v2si*)_pi32_0x7f);
+ v4sf e = _mm_cvtpi32x2_ps(mm0, mm1);
+ _mm_empty(); /* bye bye mmx */
+#else
+ emm0 = _mm_sub_epi32(emm0, *(v4si*)_pi32_0x7f);
+ v4sf e = _mm_cvtepi32_ps(emm0);
+#endif
+
+ e = _mm_add_ps(e, one);
+
+ /* part2:
+ if( x < SQRTHF ) {
+ e -= 1;
+ x = x + x - 1.0;
+ } else { x = x - 1.0; }
+ */
+ v4sf mask = _mm_cmplt_ps(x, *(v4sf*)_ps_cephes_SQRTHF);
+ v4sf tmp = _mm_and_ps(x, mask);
+ x = _mm_sub_ps(x, one);
+ e = _mm_sub_ps(e, _mm_and_ps(one, mask));
+ x = _mm_add_ps(x, tmp);
+
+
+ v4sf z = _mm_mul_ps(x,x);
+
+ v4sf y = *(v4sf*)_ps_cephes_log_p0;
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p1);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p2);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p3);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p4);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p5);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p6);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p7);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_log_p8);
+ y = _mm_mul_ps(y, x);
+
+ y = _mm_mul_ps(y, z);
+
+
+ tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q1);
+ y = _mm_add_ps(y, tmp);
+
+
+ tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
+ y = _mm_sub_ps(y, tmp);
+
+ tmp = _mm_mul_ps(e, *(v4sf*)_ps_cephes_log_q2);
+ x = _mm_add_ps(x, y);
+ x = _mm_add_ps(x, tmp);
+ x = _mm_or_ps(x, invalid_mask); // negative arg will be NAN
+ return x;
+}
+
+_PS_CONST(exp_hi, 88.3762626647949f);
+_PS_CONST(exp_lo, -88.3762626647949f);
+
+_PS_CONST(cephes_LOG2EF, 1.44269504088896341);
+_PS_CONST(cephes_exp_C1, 0.693359375);
+_PS_CONST(cephes_exp_C2, -2.12194440e-4);
+
+_PS_CONST(cephes_exp_p0, 1.9875691500E-4);
+_PS_CONST(cephes_exp_p1, 1.3981999507E-3);
+_PS_CONST(cephes_exp_p2, 8.3334519073E-3);
+_PS_CONST(cephes_exp_p3, 4.1665795894E-2);
+_PS_CONST(cephes_exp_p4, 1.6666665459E-1);
+_PS_CONST(cephes_exp_p5, 5.0000001201E-1);
+
+static inline __attribute__((__always_inline__)) v4sf exp_ps(v4sf x) {
+ v4sf tmp = _mm_setzero_ps(), fx;
+#ifdef USE_SSE2
+ v4si emm0;
+#else
+ v2si mm0, mm1;
+#endif
+ v4sf one = *(v4sf*)_ps_1;
+
+ x = _mm_min_ps(x, *(v4sf*)_ps_exp_hi);
+ x = _mm_max_ps(x, *(v4sf*)_ps_exp_lo);
+
+ /* express exp(x) as exp(g + n*log(2)) */
+ fx = _mm_mul_ps(x, *(v4sf*)_ps_cephes_LOG2EF);
+ fx = _mm_add_ps(fx, *(v4sf*)_ps_0p5);
+
+ /* how to perform a floorf with SSE: just below */
+#ifndef USE_SSE2
+ /* step 1 : cast to int */
+ tmp = _mm_movehl_ps(tmp, fx);
+ mm0 = _mm_cvttps_pi32(fx);
+ mm1 = _mm_cvttps_pi32(tmp);
+ /* step 2 : cast back to float */
+ tmp = _mm_cvtpi32x2_ps(mm0, mm1);
+#else
+ emm0 = _mm_cvttps_epi32(fx);
+ tmp = _mm_cvtepi32_ps(emm0);
+#endif
+ /* if greater, substract 1 */
+ v4sf mask = _mm_cmpgt_ps(tmp, fx);
+ mask = _mm_and_ps(mask, one);
+ fx = _mm_sub_ps(tmp, mask);
+
+ tmp = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C1);
+ v4sf z = _mm_mul_ps(fx, *(v4sf*)_ps_cephes_exp_C2);
+ x = _mm_sub_ps(x, tmp);
+ x = _mm_sub_ps(x, z);
+
+ z = _mm_mul_ps(x,x);
+
+ v4sf y = *(v4sf*)_ps_cephes_exp_p0;
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p1);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p2);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p3);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p4);
+ y = _mm_mul_ps(y, x);
+ y = _mm_add_ps(y, *(v4sf*)_ps_cephes_exp_p5);
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, x);
+ y = _mm_add_ps(y, one);
+
+ /* build 2^n */
+#ifndef USE_SSE2
+ z = _mm_movehl_ps(z, fx);
+ mm0 = _mm_cvttps_pi32(fx);
+ mm1 = _mm_cvttps_pi32(z);
+ mm0 = _mm_add_pi32(mm0, *(v2si*)_pi32_0x7f);
+ mm1 = _mm_add_pi32(mm1, *(v2si*)_pi32_0x7f);
+ mm0 = _mm_slli_pi32(mm0, 23);
+ mm1 = _mm_slli_pi32(mm1, 23);
+
+ v4sf pow2n;
+ COPY_MM_TO_XMM(mm0, mm1, pow2n);
+ _mm_empty();
+#else
+ emm0 = _mm_cvttps_epi32(fx);
+ emm0 = _mm_add_epi32(emm0, *(v4si*)_pi32_0x7f);
+ emm0 = _mm_slli_epi32(emm0, 23);
+ v4sf pow2n = _mm_castsi128_ps(emm0);
+#endif
+ y = _mm_mul_ps(y, pow2n);
+ return y;
+}
+
+_PS_CONST(minus_cephes_DP1, -0.78515625);
+_PS_CONST(minus_cephes_DP2, -2.4187564849853515625e-4);
+_PS_CONST(minus_cephes_DP3, -3.77489497744594108e-8);
+_PS_CONST(sincof_p0, -1.9515295891E-4);
+_PS_CONST(sincof_p1, 8.3321608736E-3);
+_PS_CONST(sincof_p2, -1.6666654611E-1);
+_PS_CONST(coscof_p0, 2.443315711809948E-005);
+_PS_CONST(coscof_p1, -1.388731625493765E-003);
+_PS_CONST(coscof_p2, 4.166664568298827E-002);
+_PS_CONST(cephes_FOPI, 1.27323954473516); // 4 / M_PI
+
+
+/* evaluation of 4 sines at onces, using only SSE1+MMX intrinsics so
+ it runs also on old athlons XPs and the pentium III of your grand
+ mother.
+
+ The code is the exact rewriting of the cephes sinf function.
+ Precision is excellent as long as x < 8192 (I did not bother to
+ take into account the special handling they have for greater values
+ -- it does not return garbage for arguments over 8192, though, but
+ the extra precision is missing).
+
+ Note that it is such that sinf((float)M_PI) = 8.74e-8, which is the
+ surprising but correct result.
+
+ Performance is also surprisingly good, 1.33 times faster than the
+ macos vsinf SSE2 function, and 1.5 times faster than the
+ __vrs4_sinf of amd's ACML (which is only available in 64 bits). Not
+ too bad for an SSE1 function (with no special tuning) !
+ However the latter libraries probably have a much better handling of NaN,
+ Inf, denormalized and other special arguments..
+
+ On my core 1 duo, the execution of this function takes approximately 95 cycles.
+
+ From what I have observed on the experiments with Intel AMath lib, switching to an
+ SSE2 version would improve the perf by only 10%.
+
+ Since it is based on SSE intrinsics, it has to be compiled at -O2 to
+ deliver full speed.
+*/
+static v4sf sin_ps(v4sf x) { // any x
+ v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, sign_bit, y;
+
+#ifdef USE_SSE2
+ v4si emm0, emm2;
+#else
+ v2si mm0, mm1, mm2, mm3;
+#endif
+ sign_bit = x;
+ /* take the absolute value */
+ x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
+ /* extract the sign bit (upper one) */
+ sign_bit = _mm_and_ps(sign_bit, *(v4sf*)_ps_sign_mask);
+
+ /* scale by 4/Pi */
+ y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
+
+#ifdef USE_SSE2
+ /* store the integer part of y in mm0 */
+ emm2 = _mm_cvttps_epi32(y);
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
+ y = _mm_cvtepi32_ps(emm2);
+
+ /* get the swap sign flag */
+ emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4);
+ emm0 = _mm_slli_epi32(emm0, 29);
+ /* get the polynom selection mask
+ there is one polynom for 0 <= x <= Pi/4
+ and another one for Pi/4<x<=Pi/2
+
+ Both branches will be computed.
+ */
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
+ emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
+
+ v4sf swap_sign_bit = _mm_castsi128_ps(emm0);
+ v4sf poly_mask = _mm_castsi128_ps(emm2);
+ sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
+
+#else
+ /* store the integer part of y in mm0:mm1 */
+ xmm2 = _mm_movehl_ps(xmm2, y);
+ mm2 = _mm_cvttps_pi32(y);
+ mm3 = _mm_cvttps_pi32(xmm2);
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
+ mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
+ y = _mm_cvtpi32x2_ps(mm2, mm3);
+ /* get the swap sign flag */
+ mm0 = _mm_and_si64(mm2, *(v2si*)_pi32_4);
+ mm1 = _mm_and_si64(mm3, *(v2si*)_pi32_4);
+ mm0 = _mm_slli_pi32(mm0, 29);
+ mm1 = _mm_slli_pi32(mm1, 29);
+ /* get the polynom selection mask */
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
+ mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
+ mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
+ v4sf swap_sign_bit, poly_mask;
+ COPY_MM_TO_XMM(mm0, mm1, swap_sign_bit);
+ COPY_MM_TO_XMM(mm2, mm3, poly_mask);
+ sign_bit = _mm_xor_ps(sign_bit, swap_sign_bit);
+ _mm_empty(); /* good-bye mmx */
+#endif
+
+ /* The magic pass: "Extended precision modular arithmetic"
+ x = ((x - y * DP1) - y * DP2) - y * DP3; */
+ xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
+ xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
+ xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
+ xmm1 = _mm_mul_ps(y, xmm1);
+ xmm2 = _mm_mul_ps(y, xmm2);
+ xmm3 = _mm_mul_ps(y, xmm3);
+ x = _mm_add_ps(x, xmm1);
+ x = _mm_add_ps(x, xmm2);
+ x = _mm_add_ps(x, xmm3);
+
+ /* Evaluate the first polynom (0 <= x <= Pi/4) */
+ y = *(v4sf*)_ps_coscof_p0;
+ v4sf z = _mm_mul_ps(x,x);
+
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
+ y = _mm_mul_ps(y, z);
+ y = _mm_mul_ps(y, z);
+ v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
+ y = _mm_sub_ps(y, tmp);
+ y = _mm_add_ps(y, *(v4sf*)_ps_1);
+
+ /* Evaluate the second polynom (Pi/4 <= x <= 0) */
+
+ v4sf y2 = *(v4sf*)_ps_sincof_p0;
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_mul_ps(y2, x);
+ y2 = _mm_add_ps(y2, x);
+
+ /* select the correct result from the two polynoms */
+ xmm3 = poly_mask;
+ y2 = _mm_and_ps(xmm3, y2); //, xmm3);
+ y = _mm_andnot_ps(xmm3, y);
+ y = _mm_add_ps(y,y2);
+ /* update the sign */
+ y = _mm_xor_ps(y, sign_bit);
+ return y;
+}
+
+/* almost the same as sin_ps */
+static v4sf cos_ps(v4sf x) { // any x
+ v4sf xmm1, xmm2 = _mm_setzero_ps(), xmm3, y;
+#ifdef USE_SSE2
+ v4si emm0, emm2;
+#else
+ v2si mm0, mm1, mm2, mm3;
+#endif
+ /* take the absolute value */
+ x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
+
+ /* scale by 4/Pi */
+ y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
+
+#ifdef USE_SSE2
+ /* store the integer part of y in mm0 */
+ emm2 = _mm_cvttps_epi32(y);
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
+ y = _mm_cvtepi32_ps(emm2);
+
+ emm2 = _mm_sub_epi32(emm2, *(v4si*)_pi32_2);
+
+ /* get the swap sign flag */
+ emm0 = _mm_andnot_si128(emm2, *(v4si*)_pi32_4);
+ emm0 = _mm_slli_epi32(emm0, 29);
+ /* get the polynom selection mask */
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
+ emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
+
+ v4sf sign_bit = _mm_castsi128_ps(emm0);
+ v4sf poly_mask = _mm_castsi128_ps(emm2);
+#else
+ /* store the integer part of y in mm0:mm1 */
+ xmm2 = _mm_movehl_ps(xmm2, y);
+ mm2 = _mm_cvttps_pi32(y);
+ mm3 = _mm_cvttps_pi32(xmm2);
+
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
+ mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
+
+ y = _mm_cvtpi32x2_ps(mm2, mm3);
+
+
+ mm2 = _mm_sub_pi32(mm2, *(v2si*)_pi32_2);
+ mm3 = _mm_sub_pi32(mm3, *(v2si*)_pi32_2);
+
+ /* get the swap sign flag in mm0:mm1 and the
+ polynom selection mask in mm2:mm3 */
+
+ mm0 = _mm_andnot_si64(mm2, *(v2si*)_pi32_4);
+ mm1 = _mm_andnot_si64(mm3, *(v2si*)_pi32_4);
+ mm0 = _mm_slli_pi32(mm0, 29);
+ mm1 = _mm_slli_pi32(mm1, 29);
+
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
+
+ mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
+ mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
+
+ v4sf sign_bit, poly_mask;
+ COPY_MM_TO_XMM(mm0, mm1, sign_bit);
+ COPY_MM_TO_XMM(mm2, mm3, poly_mask);
+ _mm_empty(); /* good-bye mmx */
+#endif
+ /* The magic pass: "Extended precision modular arithmetic"
+ x = ((x - y * DP1) - y * DP2) - y * DP3; */
+ xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
+ xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
+ xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
+ xmm1 = _mm_mul_ps(y, xmm1);
+ xmm2 = _mm_mul_ps(y, xmm2);
+ xmm3 = _mm_mul_ps(y, xmm3);
+ x = _mm_add_ps(x, xmm1);
+ x = _mm_add_ps(x, xmm2);
+ x = _mm_add_ps(x, xmm3);
+
+ /* Evaluate the first polynom (0 <= x <= Pi/4) */
+ y = *(v4sf*)_ps_coscof_p0;
+ v4sf z = _mm_mul_ps(x,x);
+
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
+ y = _mm_mul_ps(y, z);
+ y = _mm_mul_ps(y, z);
+ v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
+ y = _mm_sub_ps(y, tmp);
+ y = _mm_add_ps(y, *(v4sf*)_ps_1);
+
+ /* Evaluate the second polynom (Pi/4 <= x <= 0) */
+
+ v4sf y2 = *(v4sf*)_ps_sincof_p0;
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_mul_ps(y2, x);
+ y2 = _mm_add_ps(y2, x);
+
+ /* select the correct result from the two polynoms */
+ xmm3 = poly_mask;
+ y2 = _mm_and_ps(xmm3, y2); //, xmm3);
+ y = _mm_andnot_ps(xmm3, y);
+ y = _mm_add_ps(y,y2);
+ /* update the sign */
+ y = _mm_xor_ps(y, sign_bit);
+
+ return y;
+}
+
+/* since sin_ps and cos_ps are almost identical, sincos_ps could replace both of them..
+ it is almost as fast, and gives you a free cosine with your sine */
+static void sincos_ps(v4sf x, v4sf *s, v4sf *c) {
+ v4sf xmm1, xmm2, xmm3 = _mm_setzero_ps(), sign_bit_sin, y;
+#ifdef USE_SSE2
+ v4si emm0, emm2, emm4;
+#else
+ v2si mm0, mm1, mm2, mm3, mm4, mm5;
+#endif
+ sign_bit_sin = x;
+ /* take the absolute value */
+ x = _mm_and_ps(x, *(v4sf*)_ps_inv_sign_mask);
+ /* extract the sign bit (upper one) */
+ sign_bit_sin = _mm_and_ps(sign_bit_sin, *(v4sf*)_ps_sign_mask);
+
+ /* scale by 4/Pi */
+ y = _mm_mul_ps(x, *(v4sf*)_ps_cephes_FOPI);
+
+#ifdef USE_SSE2
+ /* store the integer part of y in emm2 */
+ emm2 = _mm_cvttps_epi32(y);
+
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ emm2 = _mm_add_epi32(emm2, *(v4si*)_pi32_1);
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_inv1);
+ y = _mm_cvtepi32_ps(emm2);
+
+ emm4 = emm2;
+
+ /* get the swap sign flag for the sine */
+ emm0 = _mm_and_si128(emm2, *(v4si*)_pi32_4);
+ emm0 = _mm_slli_epi32(emm0, 29);
+ v4sf swap_sign_bit_sin = _mm_castsi128_ps(emm0);
+
+ /* get the polynom selection mask for the sine*/
+ emm2 = _mm_and_si128(emm2, *(v4si*)_pi32_2);
+ emm2 = _mm_cmpeq_epi32(emm2, _mm_setzero_si128());
+ v4sf poly_mask = _mm_castsi128_ps(emm2);
+#else
+ /* store the integer part of y in mm2:mm3 */
+ xmm3 = _mm_movehl_ps(xmm3, y);
+ mm2 = _mm_cvttps_pi32(y);
+ mm3 = _mm_cvttps_pi32(xmm3);
+
+ /* j=(j+1) & (~1) (see the cephes sources) */
+ mm2 = _mm_add_pi32(mm2, *(v2si*)_pi32_1);
+ mm3 = _mm_add_pi32(mm3, *(v2si*)_pi32_1);
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_inv1);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_inv1);
+
+ y = _mm_cvtpi32x2_ps(mm2, mm3);
+
+ mm4 = mm2;
+ mm5 = mm3;
+
+ /* get the swap sign flag for the sine */
+ mm0 = _mm_and_si64(mm2, *(v2si*)_pi32_4);
+ mm1 = _mm_and_si64(mm3, *(v2si*)_pi32_4);
+ mm0 = _mm_slli_pi32(mm0, 29);
+ mm1 = _mm_slli_pi32(mm1, 29);
+ v4sf swap_sign_bit_sin;
+ COPY_MM_TO_XMM(mm0, mm1, swap_sign_bit_sin);
+
+ /* get the polynom selection mask for the sine */
+
+ mm2 = _mm_and_si64(mm2, *(v2si*)_pi32_2);
+ mm3 = _mm_and_si64(mm3, *(v2si*)_pi32_2);
+ mm2 = _mm_cmpeq_pi32(mm2, _mm_setzero_si64());
+ mm3 = _mm_cmpeq_pi32(mm3, _mm_setzero_si64());
+ v4sf poly_mask;
+ COPY_MM_TO_XMM(mm2, mm3, poly_mask);
+#endif
+
+ /* The magic pass: "Extended precision modular arithmetic"
+ x = ((x - y * DP1) - y * DP2) - y * DP3; */
+ xmm1 = *(v4sf*)_ps_minus_cephes_DP1;
+ xmm2 = *(v4sf*)_ps_minus_cephes_DP2;
+ xmm3 = *(v4sf*)_ps_minus_cephes_DP3;
+ xmm1 = _mm_mul_ps(y, xmm1);
+ xmm2 = _mm_mul_ps(y, xmm2);
+ xmm3 = _mm_mul_ps(y, xmm3);
+ x = _mm_add_ps(x, xmm1);
+ x = _mm_add_ps(x, xmm2);
+ x = _mm_add_ps(x, xmm3);
+
+#ifdef USE_SSE2
+ emm4 = _mm_sub_epi32(emm4, *(v4si*)_pi32_2);
+ emm4 = _mm_andnot_si128(emm4, *(v4si*)_pi32_4);
+ emm4 = _mm_slli_epi32(emm4, 29);
+ v4sf sign_bit_cos = _mm_castsi128_ps(emm4);
+#else
+ /* get the sign flag for the cosine */
+ mm4 = _mm_sub_pi32(mm4, *(v2si*)_pi32_2);
+ mm5 = _mm_sub_pi32(mm5, *(v2si*)_pi32_2);
+ mm4 = _mm_andnot_si64(mm4, *(v2si*)_pi32_4);
+ mm5 = _mm_andnot_si64(mm5, *(v2si*)_pi32_4);
+ mm4 = _mm_slli_pi32(mm4, 29);
+ mm5 = _mm_slli_pi32(mm5, 29);
+ v4sf sign_bit_cos;
+ COPY_MM_TO_XMM(mm4, mm5, sign_bit_cos);
+ _mm_empty(); /* good-bye mmx */
+#endif
+
+ sign_bit_sin = _mm_xor_ps(sign_bit_sin, swap_sign_bit_sin);
+
+
+ /* Evaluate the first polynom (0 <= x <= Pi/4) */
+ v4sf z = _mm_mul_ps(x,x);
+ y = *(v4sf*)_ps_coscof_p0;
+
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p1);
+ y = _mm_mul_ps(y, z);
+ y = _mm_add_ps(y, *(v4sf*)_ps_coscof_p2);
+ y = _mm_mul_ps(y, z);
+ y = _mm_mul_ps(y, z);
+ v4sf tmp = _mm_mul_ps(z, *(v4sf*)_ps_0p5);
+ y = _mm_sub_ps(y, tmp);
+ y = _mm_add_ps(y, *(v4sf*)_ps_1);
+
+ /* Evaluate the second polynom (Pi/4 <= x <= 0) */
+
+ v4sf y2 = *(v4sf*)_ps_sincof_p0;
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p1);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_add_ps(y2, *(v4sf*)_ps_sincof_p2);
+ y2 = _mm_mul_ps(y2, z);
+ y2 = _mm_mul_ps(y2, x);
+ y2 = _mm_add_ps(y2, x);
+
+ /* select the correct result from the two polynoms */
+ xmm3 = poly_mask;
+ v4sf ysin2 = _mm_and_ps(xmm3, y2);
+ v4sf ysin1 = _mm_andnot_ps(xmm3, y);
+ y2 = _mm_sub_ps(y2,ysin2);
+ y = _mm_sub_ps(y, ysin1);
+
+ xmm1 = _mm_add_ps(ysin1,ysin2);
+ xmm2 = _mm_add_ps(y,y2);
+
+ /* update the sign */
+ *s = _mm_xor_ps(xmm1, sign_bit_sin);
+ *c = _mm_xor_ps(xmm2, sign_bit_cos);
+}
+#endif /* SSE_MATHFUN_H */
diff --git a/src/thirdparty/scrappie/util.c b/src/thirdparty/scrappie/util.c
new file mode 100644
index 0000000..6895c4a
--- /dev/null
+++ b/src/thirdparty/scrappie/util.c
@@ -0,0 +1,288 @@
+// Too many calls to quantile lead to high failure rate for `scrappie raw`
+#define BANANA 1
+#include <assert.h>
+#include <err.h>
+#include <math.h>
+#include "scrappie_stdlib.h"
+#include "util.h"
+
+int argmaxf(const float *x, int n) {
+ assert(n > 0);
+ if (NULL == x) {
+ return -1;
+ }
+ int imax = 0;
+ float vmax = x[0];
+ for (int i = 1; i < n; i++) {
+ if (x[i] > vmax) {
+ vmax = x[i];
+ imax = i;
+ }
+ }
+ return imax;
+}
+
+int argminf(const float *x, int n) {
+ assert(n > 0);
+ if (NULL == x) {
+ return -1;
+ }
+ int imin = 0;
+ float vmin = x[0];
+ for (int i = 1; i < n; i++) {
+ if (x[i] > vmin) {
+ vmin = x[i];
+ imin = i;
+ }
+ }
+ return imin;
+}
+
+float valmaxf(const float *x, int n) {
+ assert(n > 0);
+ if (NULL == x) {
+ return NAN;
+ }
+ float vmax = x[0];
+ for (int i = 1; i < n; i++) {
+ if (x[i] > vmax) {
+ vmax = x[i];
+ }
+ }
+ return vmax;
+}
+
+float valminf(const float *x, int n) {
+ assert(n > 0);
+ if (NULL == x) {
+ return NAN;
+ }
+ float vmin = x[0];
+ for (int i = 1; i < n; i++) {
+ if (x[i] > vmin) {
+ vmin = x[i];
+ }
+ }
+ return vmin;
+}
+
+int floatcmp(const void *x, const void *y) {
+ float d = *(float *)x - *(float *)y;
+ if (d > 0) {
+ return 1;
+ }
+ return -1;
+}
+
+/** Quantiles from n array
+ *
+ * Using a relatively inefficent qsort resulting in O(n log n)
+ * performance but better performance is possible for small np.
+ * The array p is modified inplace, containing which quantiles to
+ * calculation on input and the quantiles on output; on error, p
+ * is filled with the value NAN.
+ *
+ * @param x An array to calculate quantiles from
+ * @param nx Length of array x
+ * @param p An array containing quantiles to calculate [in/out]
+ * @param np Length of array p
+ *
+ * @return void
+ **/
+void quantilef(const float *x, size_t nx, float *p, size_t np) {
+ if (NULL == p) {
+ return;
+ }
+ for (int i = 0; i < np; i++) {
+ assert(p[i] >= 0.0f && p[i] <= 1.0f);
+ }
+ if (NULL == x) {
+ for (int i = 0; i < np; i++) {
+ p[i] = NAN;
+ }
+ return;
+ }
+ // Sort array
+ float *space = malloc(nx * sizeof(float));
+ if (NULL == space) {
+ for (int i = 0; i < np; i++) {
+ p[i] = NAN;
+ }
+ return;
+ }
+ memcpy(space, x, nx * sizeof(float));
+ qsort(space, nx, sizeof(float), floatcmp);
+
+ // Extract quantiles
+ for (int i = 0; i < np; i++) {
+ const size_t idx = p[i] * (nx - 1);
+ const float remf = p[i] * (nx - 1) - idx;
+ if (idx < nx - 1) {
+ p[i] = (1.0 - remf) * space[idx] + remf * space[idx + 1];
+ } else {
+ // Should only occur when p is exactly 1.0
+ p[i] = space[idx];
+ }
+ }
+
+ free(space);
+ return;
+}
+
+/** Median of an array
+ *
+ * Using a relatively inefficent qsort resulting in O(n log n)
+ * performance but O(n) is possible.
+ *
+ * @param x An array to calculate median of
+ * @param n Length of array
+ *
+ * @return Median of array on success, NAN otherwise.
+ **/
+float medianf(const float *x, size_t n) {
+ float p = 0.5;
+ quantilef(x, n, &p, 1);
+ return p;
+}
+
+/** Median Absolute Deviation of an array
+ *
+ * @param x An array to calculate the MAD of
+ * @param n Length of array
+ * @param med Median of the array. If NAN then median is calculated.
+ *
+ * @return MAD of array on success, NAN otherwise.
+ **/
+float madf(const float *x, size_t n, const float *med) {
+ const float mad_scaling_factor = 1.4826;
+ if (NULL == x) {
+ return NAN;
+ }
+ if (1 == n) {
+ return 0.0f;
+ }
+
+ float *absdiff = malloc(n * sizeof(float));
+ if (NULL == absdiff) {
+ return NAN;
+ }
+
+ const float _med = (NULL == med) ? medianf(x, n) : *med;
+
+ for (size_t i = 0; i < n; i++) {
+ absdiff[i] = fabsf(x[i] - _med);
+ }
+
+ const float mad = medianf(absdiff, n);
+ free(absdiff);
+ return mad * mad_scaling_factor;
+}
+
+/** Med-MAD normalisation of an array
+ *
+ * Normalise an array using the median and MAD as measures of
+ * location and scale respectively. The array is updated inplace.
+ *
+ * @param x An array containing values to normalise
+ * @param n Length of array
+ * @return void
+ **/
+void medmad_normalise_array(float *x, size_t n) {
+ if (NULL == x) {
+ return;
+ }
+ if (1 == n) {
+ x[0] = 0.0;
+ return;
+ }
+
+ const float xmed = medianf(x, n);
+ const float xmad = madf(x, n, &xmed);
+ for (int i = 0; i < n; i++) {
+ x[i] = (x[i] - xmed) / xmad;
+ }
+}
+
+/** Studentise array using Kahan summation algorithm
+ *
+ * Studentise an array using the Kahan summation
+ * algorithm for numerical stability. The array is updated inplace.
+ * https://en.wikipedia.org/wiki/Kahan_summation_algorithm
+ *
+ * @param x An array to normalise
+ * @param n Length of array
+ * @return void
+ **/
+void studentise_array_kahan(float *x, size_t n) {
+ if (NULL == x) {
+ return;
+ }
+
+ double sum, sumsq, comp, compsq;
+ sumsq = sum = comp = compsq = 0.0;
+ for (int i = 0; i < n; i++) {
+ double d1 = x[i] - comp;
+ double sum_tmp = sum + d1;
+ comp = (sum_tmp - sum) - d1;
+ sum = sum_tmp;
+
+ double d2 = x[i] * x[i] - compsq;
+ double sumsq_tmp = sumsq + d2;
+ compsq = (sumsq_tmp - sumsq) - d2;
+ sumsq = sumsq_tmp;
+ }
+ sum /= n;
+ sumsq /= n;
+ sumsq -= sum * sum;
+
+ sumsq = sqrt(sumsq);
+
+ const float sumf = sum;
+ const float sumsqf = sumsq;
+ for (int i = 0; i < n; i++) {
+ x[i] = (x[i] - sumf) / sumsqf;
+ }
+}
+
+bool equality_array(double const * x, double const * y, size_t n, double const tol){
+
+ if(NULL == x || NULL == y){
+ return NULL == x && NULL == y;
+ }
+ for(size_t i=0 ; i < n ; i++){
+ if(fabs(x[i] - y[i]) > tol){
+ warnx("Failure at elt %zu: %f %f\n", i, x[i], y[i]);
+ return false;
+ }
+ }
+
+ return true;
+}
+
+bool equality_arrayf(float const * x, float const * y, size_t n, float const tol){
+ if(NULL == x || NULL == y){
+ return NULL == x && NULL == y;
+ }
+ for(size_t i=0 ; i < n ; i++){
+ if(fabsf(x[i] - y[i]) > tol){
+ warnx("Failure at elt %zu: %f %f\n", i, x[i], y[i]);
+ return false;
+ }
+ }
+
+ return true;
+}
+
+bool equality_arrayi(int const * x, int const * y, size_t n){
+ if(NULL == x || NULL == y){
+ return NULL == x && NULL == y;
+ }
+ for(size_t i=0 ; i < n ; i++){
+ if(x[i] != y[i]){
+ warnx("Failure at elt %zu: %d %d\n", i, x[i], y[i]);
+ return false;
+ }
+ }
+
+ return true;
+}
diff --git a/src/thirdparty/scrappie/util.h b/src/thirdparty/scrappie/util.h
new file mode 100644
index 0000000..7787f2a
--- /dev/null
+++ b/src/thirdparty/scrappie/util.h
@@ -0,0 +1,197 @@
+#pragma once
+#ifndef UTIL_H
+# define UTIL_H
+# include <immintrin.h>
+# include <math.h>
+# include <stdbool.h>
+# include <stdint.h>
+# include <stdio.h>
+# include "sse_mathfun.h"
+
+# ifdef FAST_LOG
+# define LOGFV fast_logfv
+# else
+# define LOGFV logfv
+# endif
+
+# ifdef FAST_EXP
+# define EXPFV fast_expfv
+# else
+# define EXPFV expfv
+# endif
+
+# ifdef FAST_TANH
+# define TANHFV fast_tanhfv
+# else
+# define TANHFV tanhfv
+# endif
+
+# ifdef FAST_LOGISTIC
+# define LOGISTICFV fast_logisticfv
+# else
+# define LOGISTICFV logisticfv
+# endif
+
+# ifdef FAST_ELU
+# define ELUFV fast_elufv
+# else
+# define ELUFV elufv
+# endif
+
+/* Create a vector of ones. */
+extern __inline __m128 __attribute__ ((__gnu_inline__, __always_inline__))
+ _mm_setone_ps(void) {
+ return __extension__(__m128) {
+ 1.0f, 1.0f, 1.0f, 1.0f};
+}
+
+extern __inline __m128d __attribute__ ((__gnu_inline__, __always_inline__))
+ _mm_setone_pd(void) {
+ return __extension__(__m128d) {
+ 1.0, 1.0};
+}
+
+
+/**
+ * Standard implementations
+ **/
+static inline float logisticf(float x) {
+ return 1.0f / (1.0f + expf(-x));
+}
+
+static inline float eluf(float x){
+ return (x >= 0.0f) ? x : expm1f(x);
+}
+
+// Constants for fast exp approximation. See Schraudolph (1999)
+# define _A 12102203.161561485f
+// Minimum maximum relative error
+//# define _B 1064986822.5027076f
+// No bias at zero
+# define _B 1065353216.0f
+// Minimum RMS relative error
+//# define _B 1064866803.6193156f
+// Minimum mean relative error
+//# define _B 1064807268.0425934f
+# define _BOUND 88.02969193111305
+static inline float fast_expf(float x) {
+ x = fmaxf(-_BOUND, fminf(_BOUND, x));
+ union {
+ uint32_t i;
+ float f;
+ } value = {
+ .i = (uint32_t) (_A * x + _B)};
+ return value.f;
+}
+
+static inline float fast_logisticf(float x) {
+ return 1.0 / (1.0 + fast_expf(-x));
+}
+
+static inline float fast_tanhf(float x) {
+ const float y = fast_logisticf(x + x);
+ return y + y - 1.0;
+}
+
+static inline float fast_eluf(float x){
+ return (x >= 0.0f) ? x : (fast_expf(x) - 1.0);
+}
+
+/**
+ * Vectorised functions
+ **/
+static inline __m128 __attribute__ ((__always_inline__)) expfv(__m128 x) {
+ __v4sf y = (__v4sf) x;
+ return (__m128) exp_ps(y);
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) logfv(__m128 x) {
+ __v4sf y = (__v4sf) x;
+ return (__m128) log_ps(y);
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) logisticfv(__m128 x) {
+ const __m128 ones = _mm_setone_ps();
+ return _mm_div_ps(ones, _mm_add_ps(ones, expfv(-x)));
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) tanhfv(__m128 x) {
+ const __m128 y = logisticfv(_mm_add_ps(x, x));
+ return _mm_sub_ps(_mm_add_ps(y, y), _mm_setone_ps());
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) elufv(__m128 x) {
+ if(0 == _mm_movemask_ps(x)){
+ // All positive, early return.
+ return x;
+ }
+ const __m128 mask = _mm_cmpge_ps(x, _mm_setzero_ps());
+ const __m128 y = expfv(x) - _mm_setone_ps();
+ return _mm_or_ps(_mm_and_ps(mask, x), _mm_andnot_ps(mask, y));
+}
+
+/**
+ * Fast vectorised approximations
+ **/
+static inline __m128 fast_expfv(__m128 x) {
+ const __m128 a = (__m128) (__v4sf) { _A, _A, _A, _A };
+ const __m128 b = (__m128) (__v4sf) { _B, _B, _B, _B };
+ const __m128 _bound = (__m128) (__v4sf) { _BOUND, _BOUND, _BOUND, _BOUND };
+ x = _mm_max_ps(-_bound, _mm_min_ps(_bound, x));
+
+ __m128 y = a * x + b;
+ return _mm_castsi128_ps(_mm_cvtps_epi32(y));
+}
+
+static inline __m128 fast_logfv(__m128 x) {
+# define _Alogfv 8.262958294867817e-08f
+# define _Blogfv 1064872507.1541044f
+ const __m128 a = (__m128) (__v4sf) { _Alogfv, _Alogfv, _Alogfv, _Alogfv };
+ const __m128 b = (__m128) (__v4sf) { _Blogfv, _Blogfv, _Blogfv, _Blogfv };
+ x = _mm_cvtepi32_ps(_mm_castps_si128(x));
+ return a * (x - b);
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) fast_logisticfv(__m128 x) {
+ return _mm_rcp_ps(_mm_add_ps(_mm_setone_ps(), fast_expfv(-x)));
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) fast_tanhfv(__m128 x) {
+ const __m128 y = fast_logisticfv(_mm_add_ps(x, x));
+ return _mm_sub_ps(_mm_add_ps(y, y), _mm_setone_ps());
+}
+
+static inline __m128 __attribute__ ((__always_inline__)) fast_elufv(__m128 x) {
+ if(0 == _mm_movemask_ps(x)){
+ // All positive, early return.
+ return x;
+ }
+ const __m128 mask = _mm_cmpge_ps(x, _mm_setzero_ps());
+ const __m128 y = fast_expfv(x) - _mm_setone_ps();
+ return _mm_or_ps(_mm_and_ps(mask, x), _mm_andnot_ps(mask, y));
+}
+
+int argmaxf(const float *x, int n);
+int argminf(const float *x, int n);
+float valmaxf(const float *x, int n);
+float valminf(const float *x, int n);
+
+static inline int iceil(int x, int y) {
+ return (x + y - 1) / y;
+}
+
+static inline int ifloor(int x, int y) {
+ return x / y;
+}
+
+void quantilef(const float *x, size_t nx, float *p, size_t np);
+float medianf(const float *x, size_t n);
+float madf(const float *x, size_t n, const float *med);
+void medmad_normalise_array(float *x, size_t n);
+void studentise_array_kahan(float *x, size_t n);
+
+bool equality_array(double const * x, double const * y, size_t n, double const tol);
+bool equality_arrayf(float const * x, float const * y, size_t n, float const tol);
+bool equality_arrayi(int const * x, int const * y, size_t n);
+
+#endif /* UTIL_H */
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/nanopolish.git
More information about the debian-med-commit
mailing list