diff --git a/.gitignore b/.gitignore
index cd2946a..1f78b82 100644
--- a/.gitignore
+++ b/.gitignore
@@ -45,3 +45,9 @@ $RECYCLE.BIN/
Network Trash Folder
Temporary Items
.apdisk
+
+# Manual backup files
+*.orig
+
+# binaries directory
+Linux-amd64
diff --git a/README.MD b/README.MD
index aab0e20..58350ac 100644
--- a/README.MD
+++ b/README.MD
@@ -1,680 +1,680 @@
-**We have released a new version [MECAT2](https://github.com/xiaochuanle/MECAT2). Please go and download that new version. This version will not be updated any more.**
-
-# Contents
-
-* [Introduction](#S-introduction)
-
-* [Installation](#S-installation)
-
-* [Quick Start](#S-quick-start)
-
-* [Input Format](#S-input-format)
-
-* [Program Descriptions](#S-program-description)
-
-* [Citation](#S-citation)
-
-* [Contact](#S-contact)
-
-* [Update Information](#S-update)
-
-# Introdction
-
-
-
-MECAT is an ultra-fast Mapping, Error Correction and de novo Assembly Tools for single molecula sequencing (SMRT) reads. MECAT employs novel alignment and error correction algorithms that are much more efficient than the state of art of aligners and error correction tools. MECAT can be used for effectively de novo assemblying large genomes. For example, on a 32-thread computer with 2.0 GHz CPU , MECAT takes 9.5 days to assemble a human genome based on 54x SMRT data, which is 40 times faster than the current [PBcR-Mhap pipeline](http://cbcb.umd.edu/software/pbcr/mhap/). We also use MECAT to assemble a diploid human genome based on 102x SMRT data only in 25 days. The latter assembly leads a great improvement of quality to the previous genome assembled from the 54x haploid SMRT data. MECAT performance were compared with [PBcR-Mhap pipeline](http://cbcb.umd.edu/software/pbcr/mhap/), [FALCON](https://github.com/PacificBiosciences/falcon) and [Canu(v1.3)](http://canu.readthedocs.io/en/latest/) in five real datasets. The quality of assembled contigs produced by MECAT is the same or better than that of the [PBcR-Mhap pipeline](http://cbcb.umd.edu/software/pbcr/mhap/) and [FALCON](https://github.com/PacificBiosciences/falcon). Here are some comparisons on the 32-thread computer with 2.0 GHz CPU and 512 GB RAM memory:
-
-
-
-
-
- | Genome |
- Pipeline |
- Thread number |
- Total running time (h) |
- NG50 of genome |
-
-
- | E.coli |
- FALCON |
- 16 |
- 1.21 |
- 4,635,129 |
-
-
- |
- PBcR-MHAP |
- 16 |
- 1.29 |
- 4,652,272 |
-
-
- |
- Canu |
- 16 |
- 0.71 |
- 4,648,002 |
-
-
- |
- MECAT |
- 16 |
- 0.24 |
- 4,649,626 |
-
-
- | Yeast |
- FALCON |
- 16 |
- 2.16 |
- 587,169 |
-
-
- |
- PBcR-MHAP |
- 16 |
- 4.2 |
- 818,229 |
-
-
- |
- Canu |
- 16 |
- 5.11 |
- 739,902 |
-
-
- |
- MECAT |
- 16 |
- 0.91 |
- 929,350 |
-
-
- | A.thaliana |
- FALCON |
- 16 |
- 223.84 |
- 7,583,032 |
-
-
- |
- PBcR-MHAP |
- 16 |
- 188.7 |
- 9,610,192 |
-
-
- |
- Canu |
- 16 |
- 118.57 |
- 8,315,338 |
-
-
- |
- MECAT |
- 16 |
- 10.73 |
- 12600961 |
-
-
- | D.melanogaster |
- FALCON |
- 16 |
- 140.72 |
- 15,664,372 |
-
-
- |
- PBcR-MHAP |
- 16 |
- 101.22 |
- 13,627,256 |
-
-
- |
- Canu |
- 16 |
- 69.34 |
- 14,179,324 |
-
-
- |
- MECAT |
- 16 |
- 9.58 |
- 18,111,159 |
-
-
- | Human(54X) |
- PBcR-MHAH(f) |
- 32 |
- 5750 |
- 1,857,788 |
-
-
- |
- PBcR-MHAH(s) |
- 32 |
- 13000 |
- 4,320,471 |
-
-
- |
- MECAT |
- 32 |
- 230.54 |
- 4,878,957 |
-
-
-
+
+
+ | Genome |
+ Pipeline |
+ Thread number |
+ Total running time (h) |
+ NG50 of genome |
+
+
+ | E.coli |
+ FALCON |
+ 16 |
+ 1.21 |
+ 4,635,129 |
+
+
+ |
+ PBcR-MHAP |
+ 16 |
+ 1.29 |
+ 4,652,272 |
+
+
+ |
+ Canu |
+ 16 |
+ 0.71 |
+ 4,648,002 |
+
+
+ |
+ MECAT |
+ 16 |
+ 0.24 |
+ 4,649,626 |
+
+
+ | Yeast |
+ FALCON |
+ 16 |
+ 2.16 |
+ 587,169 |
+
+
+ |
+ PBcR-MHAP |
+ 16 |
+ 4.2 |
+ 818,229 |
+
+
+ |
+ Canu |
+ 16 |
+ 5.11 |
+ 739,902 |
+
+
+ |
+ MECAT |
+ 16 |
+ 0.91 |
+ 929,350 |
+
+
+ | A.thaliana |
+ FALCON |
+ 16 |
+ 223.84 |
+ 7,583,032 |
+
+
+ |
+ PBcR-MHAP |
+ 16 |
+ 188.7 |
+ 9,610,192 |
+
+
+ |
+ Canu |
+ 16 |
+ 118.57 |
+ 8,315,338 |
+
+
+ |
+ MECAT |
+ 16 |
+ 10.73 |
+ 12600961 |
+
+
+ | D.melanogaster |
+ FALCON |
+ 16 |
+ 140.72 |
+ 15,664,372 |
+
+
+ |
+ PBcR-MHAP |
+ 16 |
+ 101.22 |
+ 13,627,256 |
+
+
+ |
+ Canu |
+ 16 |
+ 69.34 |
+ 14,179,324 |
+
+
+ |
+ MECAT |
+ 16 |
+ 9.58 |
+ 18,111,159 |
+
+
+ | Human(54X) |
+ PBcR-MHAH(f) |
+ 32 |
+ 5750 |
+ 1,857,788 |
+
+
+ |
+ PBcR-MHAH(s) |
+ 32 |
+ 13000 |
+ 4,320,471 |
+
+
+ |
+ MECAT |
+ 32 |
+ 230.54 |
+ 4,878,957 |
+
+
+
// list<>
+#include // ostringstream
+#include // string
+#include // chmod()
+#include // ... unlink()
+#include // pair<>
+#include // vector<>
+
+static void grid_start(const char* const prog, const ReadsCorrectionOptions &options, const int i, const std::string& exit_file) {
+ // make sure exit marker is not present
+ unlink(exit_file.c_str());
+ // create grid script, have grid run it
+ ReadsCorrectionOptions new_options(options);
+ new_options.job_index = i;
+ new_options.grid_options = NULL;
+ new_options.grid_options_split = NULL;
+ std::string name("m2cns.");
+ if (i == -1) {
+ name += "split";
+ new_options.num_threads = 0; // flag as grid spin-off
+ } else {
+ std::ostringstream x;
+ x << i;
+ name += x.str();
+ }
+ std::string script_file(name + ".sh");
+ // as we might not have write permission, delete it
+ unlink(script_file.c_str());
+ std::ofstream out;
+ open_fstream(out, script_file.c_str(), std::ios::out);
+ std::string new_options_str;
+ make_options(new_options, new_options_str);
+ out << "#!/bin/bash\nset -e\ntrap 'touch " << exit_file << "' EXIT\nulimit -c 0\ntime " << prog << new_options_str << "\n";
+ if (!out) {
+ std::cerr << "Error writing to " << script_file << "\n";
+ exit(1);
+ }
+ close_fstream(out);
+ chmod(script_file.c_str(), S_IRUSR | S_IXUSR);
+ std::string cmd(i == -1 && options.grid_options_split ? options.grid_options_split : options.grid_options);
+ cmd += " " + name + " " + script_file;
+ assert(system(cmd.c_str()) == 0);
+}
+
+// exit files get modified during loop
+
+static void wait_for_files(std::list& exit_files) {
+ const std::list::const_iterator end_a(exit_files.end());
+ while (!exit_files.empty()) {
+ sleep(60);
+ std::list::iterator a(exit_files.begin());
+ while (a != end_a) {
+ if (access(a->c_str(), F_OK) == 0) {
+ // now test to see if it worked, by removed ".exit" from end
+ a->resize(a->size() - 5);
+ if (access(a->c_str(), F_OK) != 0) { // failed!
+ ERROR("Failed: %s does not exist", a->c_str());
+ }
+ a = exit_files.erase(a);
+ } else {
+ ++a;
+ }
+ }
+ }
+}
+
+static void merge_results(const char* const output, const std::list& files) {
+ std::string out_tmp(output);
+ out_tmp += ".tmp";
+ std::ofstream out(out_tmp.c_str());
+ std::list::const_iterator a(files.begin());
+ const std::list::const_iterator end_a(files.end());
+ for (; a != end_a; ++a) {
+ std::ifstream in(a->c_str());
+ out << in.rdbuf();
+ if (!in) {
+ std::cerr << "Error reading from " << *a << "\n";
+ exit(1);
+ } else if (!out) {
+ std::cerr << "Error writing to " << out_tmp << "\n";
+ exit(1);
+ }
+ }
+ out.close();
+ if (rename(out_tmp.c_str(), output) == -1) {
+ std::cerr << "Could not rename concatenated output file: " << out_tmp << "\n";
exit(1);
}
- if (rco.print_usage_info) {
+}
+
+int main(int argc, char** argv) {
+ ReadsCorrectionOptions rco;
+ if (!parse_arguments(argc, argv, rco)) {
+ print_usage(argv[0]);
+ exit(1);
+ } else if (rco.print_usage_info) {
print_usage(argv[0]);
exit(0);
}
-
- if (rco.input_type == INPUT_TYPE_CAN)
- {
- return reads_correction_can(rco);
+ // partition once up front, to allow for grid jobs
+ if (rco.job_index == -1) {
+ if (access(rco.corrected_reads, F_OK) == 0) { // full run already done
+ return 0;
+ } else if (access("partition.done", F_OK) == 0) { // split already done
+ } else if (rco.grid_options || rco.grid_options_split) {
+ grid_start(argv[0], rco, -1, "partition.done.exit");
+ std::list partition_results;
+ partition_results.push_back("partition.done.exit");
+ wait_for_files(partition_results);
+ } else {
+ if (rco.input_type == INPUT_TYPE_CAN) {
+ // this speeds up candidate starts
+ const int64_t n_reads(PackedDB::convert_fasta_to_db(rco.reads, "fasta.db", rco.min_size));
+ if (rco.reads_to_correct <= 0 || n_reads < rco.reads_to_correct) {
+ rco.reads_to_correct = n_reads;
+ }
+ // can't reorder reads for memory efficiency - takes
+ // too much memory to hold read-read pairings
+ partition_candidates(rco.m4, "fasta.db", rco.batch_size, rco.num_partition_files, rco.reads_to_correct);
+ } else {
+ partition_m4records(rco.m4, rco.min_mapping_ratio - 0.02, rco.batch_size, rco.min_size, rco.num_partition_files);
+ }
+ if (rco.num_threads == 0) { // flag saying we're a grid run
+ return 0;
+ }
+ }
}
- else
- {
- return reads_correction_m4(rco);
+ if (rco.grid_options == NULL) { // single process run
+ if (rco.input_type == INPUT_TYPE_CAN) {
+ return reads_correction_can(rco);
+ } else {
+ return reads_correction_m4(rco);
+ }
+ } else { // grid run
+ // get number of partitions
+ std::string idx_file_name;
+ generate_partition_index_file_name(rco.m4, idx_file_name);
+ std::vector partition_file_vec;
+ load_partition_files_info(idx_file_name.c_str(), partition_file_vec);
+ std::list exit_files, results;
+ for (size_t i(0); i != partition_file_vec.size(); ++i) {
+ std::ostringstream os;
+ os << rco.corrected_reads << "." << i;
+ const std::string done_file(os.str());
+ const std::string exit_file(os.str() + ".exit");
+ results.push_back(done_file);
+ exit_files.push_back(exit_file);
+ if (access(os.str().c_str(), F_OK) != 0) {
+ grid_start(argv[0], rco, i, exit_file);
+ if (rco.grid_start_delay) {
+ sleep(rco.grid_start_delay);
+ }
+ }
+ }
+ wait_for_files(exit_files);
+ merge_results(rco.corrected_reads, results);
}
return 0;
}
diff --git a/src/mecat2cns/mecat2cns.mk b/src/mecat2cns/mecat2cns.mk
index c96a2db..abb9110 100644
--- a/src/mecat2cns/mecat2cns.mk
+++ b/src/mecat2cns/mecat2cns.mk
@@ -7,7 +7,6 @@ endif
TARGET := mecat2cns
SOURCES := main.cpp \
- argument.cpp \
dw.cpp \
MECAT_AlnGraphBoost.C \
mecat_correction.cpp \
@@ -16,11 +15,15 @@ SOURCES := main.cpp \
reads_correction_aux.cpp \
reads_correction_can.cpp \
reads_correction_m4.cpp \
+ packed_db.cpp \
SRC_INCDIRS := . libboost
-TGT_LDFLAGS := -L${TARGET_DIR}
-TGT_LDLIBS := -lmecat
-TGT_PREREQS := libmecat.a
+# make sure large files are okay (and off_t is 8 bytes);
+# requires c++11 or higher for headers
+TGT_CXXFLAGS := -static -D _FILE_OFFSET_BITS=64 -std=c++11 -O3 -mfpmath=sse -march=native -flto -fno-fat-lto-objects -fno-builtin -mmmx -msse -msse2 -mssse3 -msse4.1 -msse4.2 -mavx -maes -mpopcnt -mfxsr -mxsave -mxsaveopt
+TGT_LDFLAGS := -static -L${TARGET_DIR} ${TGT_CXXFLAGS}
+TGT_LDLIBS := -lmecat
+TGT_PREREQS := libmecat.a
SUBMAKEFILES :=
diff --git a/src/mecat2cns/mecat_correction.cpp b/src/mecat2cns/mecat_correction.cpp
index 0fc3942..3fc1389 100644
--- a/src/mecat2cns/mecat_correction.cpp
+++ b/src/mecat2cns/mecat_correction.cpp
@@ -1,75 +1,73 @@
#include "mecat_correction.h"
-#include "MECAT_AlnGraphBoost.H"
-
-using namespace ns_banded_sw;
+#include // assert()
+#include // numeric_limits::max()
+#include // uint8_t
-namespace ns_meap_cns {
+#include "MECAT_AlnGraphBoost.H"
+#include "reads_correction_aux.h" // GAP_CHAR
#define FMAT 1
#define FDEL 2
#define FINS 4
#define UNDS 8
-inline uint1
-identify_one_consensus_item(CnsTableItem& cns_item, const int min_cov)
-{
- uint1 ident = 0;
- int cov = cns_item.mat_cnt + cns_item.ins_cnt;
- if (cns_item.mat_cnt >= cov * 0.8) ident |= FMAT;
- if (cns_item.ins_cnt >= cov * 0.8) ident |= FINS;
- if (!ident) ident |= UNDS;
- if (cns_item.del_cnt >= cov * 0.4) ident |= FDEL;
+// returns type of coverage present
+
+static inline uint8_t identify_one_consensus_item(const CnsTableItem& cns_item) {
+ const int cov((cns_item.mat_cnt + cns_item.ins_cnt) * 0.8);
+ uint8_t ident;
+ if (cns_item.mat_cnt >= cov) { // coverage is 80% or more matches
+ ident = FMAT;
+ } else if (cns_item.ins_cnt >= cov) { // coverage is 80% or more inserts
+ ident = FINS;
+ } else { // neither of the above
+ ident = UNDS;
+ }
+ if (cns_item.del_cnt * 2 >= cov) { // deletes are 40% or more than the coverage
+ ident |= FDEL;
+ }
return ident;
}
-struct CompareOverlapByOverlapSize
-{
- bool operator()(const Overlap& a ,const Overlap& b)
- {
- const index_t ovlp_a = std::max(a.qend - a.qoff, a.send - a.soff);
- const index_t ovlp_b = std::max(b.qend - b.qoff, b.send - b.soff);
- return ovlp_a > ovlp_b;
+struct CompareOverlapByOverlapSize {
+ bool operator()(const Overlap& a, const Overlap& b) {
+ const int64_t ovlp_a(std::max(a.qend - a.qoff, a.send - a.soff));
+ const int64_t ovlp_b(std::max(b.qend - b.qoff, b.send - b.soff));
+ return ovlp_b < ovlp_a;
}
};
-void
-meap_add_one_aln(const std::string& qaln, const std::string& saln, index_t start_soff, CnsTableItem* cns_table, const char* org_seq)
-{
- r_assert(qaln.size() == saln.size());
- const index_t aln_size = qaln.size();
- index_t i = 0;
- const char kGap = '-';
- while (i < aln_size)
- {
- const char q = qaln[i];
- const char s = saln[i];
- if (q == kGap && s == kGap) { ++i; continue; }
-
- if (q == s) { ++cns_table[start_soff].mat_cnt; cns_table[start_soff].base = s; ++start_soff; ++i; }
- else if (q == kGap) { ++cns_table[start_soff].ins_cnt; ++start_soff; ++i; }
- else
- {
- r_assert(s == kGap);
- index_t j = i + 1;
- while (j < aln_size && saln[j] == kGap) ++j;
+static void meap_add_one_aln(const std::string& qaln, const std::string& saln, int64_t start_soff, std::vector& cns_table) {
+ assert(qaln.size() == saln.size());
+ const int64_t aln_size(qaln.size());
+ for (int64_t i(0); i < aln_size;) {
+ const char q(qaln[i]);
+ const char s(saln[i]);
+ if (q == GAP_CHAR && s == GAP_CHAR) { // skip
+ ++i;
+ } else if (q == s) { // match
+ ++cns_table[start_soff].mat_cnt;
+ cns_table[start_soff].base = s;
+ ++start_soff;
+ ++i;
+ } else if (q == GAP_CHAR) { // insert
+ ++cns_table[start_soff].ins_cnt;
+ ++start_soff;
+ ++i;
+ } else { // delete
+ assert(s == GAP_CHAR);
+ for (++i; i < aln_size && saln[i] == GAP_CHAR; ++i) { }
++cns_table[start_soff - 1].del_cnt;
- i = j;
}
}
}
-void
-meap_cns_one_indel(const int sb, const int se, CnsAlns& cns_vec,
- const int min_cov, std::string& aux_qstr,
- std::string& aux_tstr, std::string& cns)
-{
- AlnGraphBoost ag(se - sb + 1);
+static void meap_cns_one_indel(const int sb, const int se, CnsAlns& cns_vec, const int min_cov, std::string& aux_qstr, std::string& aux_tstr, std::string& cns) {
+ ns_meap_cns::AlnGraphBoost ag(se - sb + 1);
int sb_out;
- for (CnsAln* iter = cns_vec.begin(); iter != cns_vec.end(); ++iter)
- {
- if ((*iter).retrieve_aln_subseqs(sb, se, aux_qstr, aux_tstr, sb_out))
- {
+ for (std::vector::iterator a(cns_vec.begin()); a != cns_vec.end(); ++a) {
+ if (a->retrieve_aln_subseqs(sb, se, aux_qstr, aux_tstr, sb_out)) {
ag.addAln(aux_qstr, aux_tstr, sb_out - sb + 1);
}
}
@@ -77,88 +75,100 @@ meap_cns_one_indel(const int sb, const int se, CnsAlns& cns_vec,
ag.consensus(min_cov * 0.4, cns);
}
-void
-meap_consensus_one_segment(CnsTableItem* cns_list, const int cns_list_size,
- uint1* cns_id_vec,
- int start_soff, CnsAlns& cns_vec,
- std::string& aux_qstr, std::string& aux_tstr,
- std::string& target, const int min_cov)
-{
- for (int i = 0; i < cns_list_size; ++i) cns_id_vec[i] = identify_one_consensus_item(cns_list[i], min_cov);
- int i = 0, j;
+static void meap_consensus_one_segment(const std::vector& cns_list, const int cns_list_size, std::vector& cns_id_vec, const int start_soff, CnsAlns& cns_vec, std::string& aux_qstr, std::string& aux_tstr, std::string& target) {
+ if (cns_id_vec.size() < static_cast(cns_list_size)) {
+ cns_id_vec.resize(cns_list_size);
+ }
+ // get types of coverage
+ for (int i(0); i < cns_list_size; ++i) {
+ cns_id_vec[i] = identify_one_consensus_item(cns_list[start_soff + i]);
+ }
std::string cns;
target.clear();
- while (i < cns_list_size && !(cns_id_vec[i] & FMAT)) ++i;
- while (i < cns_list_size)
- {
- target.push_back(cns_list[i].base);
- j = i + 1;
- while (j < cns_list_size && !(cns_id_vec[j] & FMAT)) ++j;
-
- bool need_refinement = false;
- for (int k = i; k < j; ++k)
- if ((cns_id_vec[k] & UNDS) || (cns_id_vec[k] & FDEL)) { need_refinement = true; break; }
- if (need_refinement)
- {
- meap_cns_one_indel(i + start_soff, j + start_soff, cns_vec, cns_list[i].mat_cnt + cns_list[i].ins_cnt, aux_qstr, aux_tstr, cns);
- if (cns.size() > 2) target.append(cns.data() + 1, cns.size() - 2);
+ const uint8_t unds_or_fdel(UNDS | FDEL); // questionable coverage types
+ int i(0);
+ // advance to matching coverage
+ for (; i < cns_list_size && !(cns_id_vec[i] & FMAT); ++i) { }
+ while (i < cns_list_size) {
+ target.push_back(cns_list[start_soff + i].base);
+ const int start(i);
+ // advance to next matching coverage
+ for (++i; i < cns_list_size && !(cns_id_vec[i] & FMAT); ++i) { }
+ int need_refinement(0);
+ // check to see if anything in-between has questionable coverage
+ for (int k(start); k != i; ++k) {
+ if (cns_id_vec[k] & unds_or_fdel) {
+ need_refinement = 1;
+ break;
+ }
+ }
+ if (need_refinement) {
+ meap_cns_one_indel(start_soff + start, start_soff + i, cns_vec, cns_list[start_soff + start].mat_cnt + cns_list[start_soff + start].ins_cnt, aux_qstr, aux_tstr, cns);
+ // trim first and last as they have good coverage
+ if (cns.size() > 2) {
+ target.append(cns.data() + 1, cns.size() - 2);
+ }
}
- i = j;
}
}
-struct CmpMappingRangeBySoff
-{
- bool operator()(const MappingRange& a, const MappingRange& b)
- {
- return (a.start == b.start) ? (a.end > b.end) : (a.start < b.start);
+struct CmpMappingRangeBySoff {
+ bool operator()(const MappingRange& a, const MappingRange& b) {
+ if (a.start != b.start) {
+ return a.start < b.start;
+ } else {
+ return b.end < a.end;
+ }
}
};
-void
-get_effective_ranges(std::vector& mranges, std::vector& eranges, const int read_size, const int min_size)
-{
+static void get_effective_ranges(std::vector& mranges, std::vector& eranges, const int read_size, const int min_size) {
eranges.clear();
- if (mranges.size() == 0) return;
- std::vector::iterator iter;
- for (iter = mranges.begin(); iter != mranges.end(); ++iter)
- if (iter->start <= 500 && read_size - iter->end <= 500)
- {
+ if (mranges.size() == 0) {
+ return;
+ }
+ // see if any ranges are effectively the whole read
+ std::vector::const_iterator a(mranges.begin());
+ const std::vector::const_iterator end_a(mranges.end());
+ for (; a != end_a; ++a) {
+ if (a->start <= 500 && read_size - a->end <= 500) {
eranges.push_back(MappingRange(0, read_size));
return;
}
-
+ }
std::sort(mranges.begin(), mranges.end(), CmpMappingRangeBySoff());
- const int nr = mranges.size();
- int i = 0, j;
- int left = mranges[i].start, right;
- while (i < nr)
- {
- j = i + 1;
- while (j < nr && mranges[j].end <= mranges[i].end) ++j;
- if (j == nr)
- {
- right = mranges[i].end;
- if (right - left >= min_size * 0.95) eranges.push_back(MappingRange(left, right));
+ const int nr(mranges.size());
+ // -1 so we can use > instead of >=
+ const int min_size_95(ceil(min_size * 0.95) - 1);
+ int i(0), left(mranges[0].start);
+ for (;;) {
+ const int right(mranges[i].end);
+ // include ranges we completely overlap
+ for (++i; i < nr && mranges[i].end <= right; ++i) { }
+ if (i == nr) {
+ if (right - left > min_size_95) {
+ eranges.push_back(MappingRange(left, right));
+ }
break;
}
- if (mranges[i].end - mranges[j].start < 1000)
- {
- right = std::min(mranges[i].end, mranges[j].start);
- if (right - left >= min_size * 0.95) eranges.push_back(MappingRange(left, right));
- left = std::max(mranges[i].end, mranges[j].start);
+ // if the next non-contained range overlaps by 1k or more,
+ // treat as part of this range and keep extending
+ if (right - mranges[i].start < 1000) {
+ // truncate current range at overlap start
+ const int eff_right(std::min(right, mranges[i].start));
+ if (eff_right - left > min_size_95) {
+ eranges.push_back(MappingRange(left, eff_right));
+ }
+ // put some space between the effective ranges
+ left = std::max(right, mranges[i].start);
}
- i = j;
}
}
-void
-output_cns_result(std::vector& cns_results,
- CnsResult& cr,
- const index_t beg,
- const index_t end,
- std::string& cns_seq)
-{
+// breaks output sequence into chunks of no more than MaxSeqSize (if needed);
+// split chunks will have an OvlpSize overlap
+
+static void output_cns_result(std::vector& cns_results, CnsResult& cr, const int64_t beg, const int64_t end, const std::string& cns_seq) {
const size_t MaxSeqSize = 60000;
const size_t OvlpSize = 10000;
// BlkSize must be >= OvlpSize
@@ -187,139 +197,156 @@ output_cns_result(std::vector& cns_results,
}
}
-inline bool
-check_ovlp_mapping_range(const int qb, const int qe, const int qs,
- const int sb, const int se, const int ss,
- double ratio)
-{
- const int oq = qe - qb;
- const int qqs = qs * ratio;
- const int os = se - sb;
- const int qss = ss * ratio;
- return oq >= qqs || os >= qss;
+static inline bool check_ovlp_mapping_range(const int qb, const int qe, const int qs, const int sb, const int se, const int ss, double ratio) {
+ return qe - qb >= qs * ratio || se - sb >= ss * ratio;
}
-void
-consensus_worker(CnsTableItem* cns_table,
- uint1* id_list,
- CnsAlns& cns_vec,
- std::string& aux_qstr,
- std::string& aux_tstr,
- std::vector& eranges,
- const int min_cov,
- const int min_size,
- const int read_id,
- std::vector& cns_results)
-{
- index_t beg = 0, end;
+// look for areas of high coverage of about min_size or more,
+// improve them and stick on the results pile
+
+static void consensus_worker(const std::vector& cns_table, std::vector& id_list, CnsAlns& cns_vec, std::string& aux_qstr, std::string& aux_tstr, const std::vector& eranges, const int min_cov, const int min_size, const int read_id, std::vector& cns_results) {
CnsResult cns_result;
- std::string cns_seq;
cns_result.id = read_id;
- std::vector::iterator miter;
- for (miter = eranges.begin(); miter != eranges.end(); ++miter)
- {
- int L = miter->start, R = miter->end;
- beg = L;
- while (beg < R)
- {
- while (beg < R && cns_table[beg].mat_cnt + cns_table[beg].ins_cnt < min_cov) ++beg;
- end = beg + 1;
- while (end < R && cns_table[end].mat_cnt + cns_table[end].ins_cnt >= min_cov) ++end;
- if (end - beg >= 0.95 * min_size)
- {
- meap_consensus_one_segment(cns_table + beg, end - beg, id_list,
- beg, cns_vec, aux_qstr, aux_tstr, cns_seq, min_cov);
-
- if (cns_seq.size() >= min_size) output_cns_result(cns_results, cns_result, beg, end, cns_seq);
+ const int64_t min_size_95(ceil(0.95 * min_size));
+ std::string cns_seq;
+ std::vector::const_iterator a(eranges.begin());
+ const std::vector::const_iterator end_a(eranges.end());
+ for (; a != end_a; ++a) {
+ const int end_i(a->end);
+ for (int64_t i(a->start); i < end_i;) {
+ // find start of next high coverage area
+ for (; i < end_i && cns_table[i].mat_cnt + cns_table[i].ins_cnt < min_cov; ++i) { }
+ const int64_t start(i);
+ // find end of high coverage area
+ for (++i; i < end_i && cns_table[i].mat_cnt + cns_table[i].ins_cnt >= min_cov; ++i) { }
+ if (i - start >= min_size_95) {
+ meap_consensus_one_segment(cns_table, i - start, id_list, start, cns_vec, aux_qstr, aux_tstr, cns_seq);
+ if (cns_seq.size() >= static_cast(min_size)) {
+ output_cns_result(cns_results, cns_result, start, i, cns_seq);
+ }
}
-
- beg = end;
}
}
}
-void
-consensus_one_read_m4_pacbio(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid)
-{
- PackedDB& reads = *ctd->reads;
- ExtensionCandidate* overlaps = ctd->candidates;
- DiffRunningData* drd_s = ctd->drd_s;
- DiffRunningData* drd = NULL;
- M5Record* m5 = ctd->m5;
- CnsAlns& cns_vec = ctd->cns_alns;
- std::vector& cns_results = ctd->cns_results;
- const index_t read_size = overlaps[sid].ssize;
- std::vector& qstr = ctd->query;
- std::vector& tstr = ctd->target;
- tstr.resize(read_size);
- reads.GetSequence(read_id, true, tstr.data(), read_size);
- std::string& nqstr = ctd->qaln;
- std::string& ntstr = ctd->saln;
- const int min_align_size = ctd->rco.min_align_size;
- const int max_added = 60;
+static void decode_and_append_sequence(std::string& s, const std::string& seq, int64_t i, const int64_t end_i) {
+ s.reserve(s.size() + end_i - i);
+ for (; i < end_i; ++i) {
+ s += "ACGT"[static_cast(seq[i])];
+ }
+}
- index_t L, R;
- if (eid - sid <= max_added)
- {
- L = sid;
- R = eid;
+// same as consensus_worker, but produces entire read as one entry;
+// uncorrected sections are just copied as is;
+
+static void consensus_worker_one_read(const std::vector& cns_table, std::vector& id_list, CnsAlns& cns_vec, std::string& aux_qstr, std::string& aux_tstr, const std::vector& eranges, const int min_cov, const int min_size, const int read_id, const std::string& tstr, std::vector& cns_results) {
+ CnsResult cns_result;
+ cns_result.id = read_id;
+ const int64_t min_size_95(std::max(ceil(0.95 * min_size), 1) - 1);
+ std::string cns_seq;
+ std::vector::const_iterator a(eranges.begin());
+ const std::vector::const_iterator end_a(eranges.end());
+ std::vector::const_iterator last_a(eranges.end());
+ for (; a != end_a; last_a = a++) {
+ if (last_a != end_a) { // add in-between range to cns_result
+ decode_and_append_sequence(cns_result.seq, tstr, last_a->end, a->start);
+ } else if (a->start > 0) { // add beginning of read
+ decode_and_append_sequence(cns_result.seq, tstr, 0, a->start);
+ }
+ const int begin_i(a->start - 1);
+ const int end_i(a->end);
+ for (int64_t i(begin_i); i < end_i;) {
+ // find start of next high coverage area
+ const int64_t last_end(i != begin_i ? i : a->start);
+ for (++i; i < end_i && cns_table[i].mat_cnt + cns_table[i].ins_cnt < min_cov; ++i) { }
+ // add low coverage area as-is
+ decode_and_append_sequence(cns_result.seq, tstr, last_end, i);
+ if (i == end_i) {
+ break;
+ }
+ const int64_t start(i);
+ // find end of high coverage area
+ for (++i; i < end_i && cns_table[i].mat_cnt + cns_table[i].ins_cnt >= min_cov; ++i) { }
+ if (i - start > min_size_95) {
+ meap_consensus_one_segment(cns_table, i - start, id_list, start, cns_vec, aux_qstr, aux_tstr, cns_seq);
+ if (cns_seq.size() >= static_cast(min_size)) {
+ // add corrected sequence
+ cns_result.seq += cns_seq;
+ continue;
+ }
+ }
+ // add uncorrected sequence
+ decode_and_append_sequence(cns_result.seq, tstr, start, i);
+ }
}
- else
- {
- L = sid;
- R = L + max_added;
- std::sort(overlaps + sid, overlaps + eid, CompareOverlapByOverlapSize());
+ // add end of read
+ if (last_a != end_a) {
+ decode_and_append_sequence(cns_result.seq, tstr, last_a->end, tstr.size());
}
+ cns_result.range[0] = 0;
+ cns_result.range[1] = cns_result.seq.size();
+ cns_results.push_back(cns_result);
+}
- CnsTableItem* cns_table = ctd->cns_table;
- std::for_each(cns_table, cns_table + read_size, CnsTableItemCleaner());
+void consensus_one_read_m4_pacbio(ConsensusThreadData& ctd, ConsensusPerThreadData &pctd, const int64_t read_id, const int64_t sid, const int64_t eid) {
+ PackedDB& reads(ctd.reads);
+ ExtensionCandidate* overlaps((ExtensionCandidate*)pctd.candidates);
+ DiffRunningData& drd(pctd.drd);
+ M5Record& m5(pctd.m5);
+ CnsAlns& cns_vec(pctd.cns_alns);
+ std::vector& cns_results(pctd.cns_results);
+ const int64_t read_size(overlaps[read_id].ssize);
+ std::string& qstr(pctd.query);
+ std::string& tstr(pctd.target);
+ reads.GetSequence(read_id, 1, tstr);
+ std::string& nqstr(pctd.qaln);
+ std::string& ntstr(pctd.saln);
+ const int min_align_size(ctd.rco.min_align_size);
+ const double error_rate(ctd.rco.error_rate);
+ const int max_added(60);
+ std::vector& cns_table(pctd.cns_table);
+ cns_table.assign(read_size, CnsTableItem()); // reset table
cns_vec.clear();
- for (index_t i = L; i < R; ++i)
- {
- Overlap& ovlp = overlaps[i];
- qstr.resize(ovlp.qsize);
- reads.GetSequence(ovlp.qid, ovlp.qdir == FWD, qstr.data(), ovlp.qsize);
- index_t qext = ovlp.qext;
- index_t sext = ovlp.sext;
- if (ovlp.qdir == REV) qext = ovlp.qsize - 1 - qext;
- drd = drd_s;
- bool r = GetAlignment(qstr.data(), qext, qstr.size(), tstr.data(), sext, tstr.size(), drd, *m5, 0.15, min_align_size);
- if (r)
- {
- normalize_gaps(m5qaln(*m5), m5saln(*m5), strlen(m5qaln(*m5)), nqstr, ntstr, true);
- meap_add_one_aln(nqstr, ntstr, m5soff(*m5), cns_table, tstr.data());
- cns_vec.add_aln(m5soff(*m5), m5send(*m5), nqstr, ntstr);
+ const int64_t L(sid);
+ const int64_t R(eid - sid <= max_added ? eid : L + max_added);
+ if (eid - sid > max_added) { // only use largest max_added overlaps
+ std::sort(overlaps + sid, overlaps + eid, CompareOverlapByOverlapSize());
+ }
+ for (int64_t i(L); i < R; ++i) {
+ Overlap& ovlp(overlaps[i]);
+ reads.GetSequence(ovlp.qid, ovlp.qdir == FWD, qstr);
+ const int64_t qext(ovlp.qdir == FWD ? ovlp.qext : ovlp.qsize - 1 - ovlp.qext);
+ const int r(GetAlignment(qstr, qext, tstr, ovlp.sext, drd, m5, error_rate, min_align_size));
+ if (r) {
+ normalize_gaps(m5.qaln, m5.saln, nqstr, ntstr, 1);
+ meap_add_one_aln(nqstr, ntstr, m5.soff, cns_table);
+ cns_vec.add_aln(m5.soff, m5.send, nqstr, ntstr);
}
}
-
std::vector mranges, eranges;
cns_vec.get_mapping_ranges(mranges);
- get_effective_ranges(mranges, eranges, read_size, ctd->rco.min_size);
-
- consensus_worker(cns_table, ctd->id_list, cns_vec, nqstr, ntstr, eranges, ctd->rco.min_cov, ctd->rco.min_size, read_id, cns_results);
+ get_effective_ranges(mranges, eranges, read_size, ctd.rco.min_size);
+ consensus_worker(cns_table, pctd.id_list, cns_vec, nqstr, ntstr, eranges, ctd.rco.min_cov, ctd.rco.min_size, read_id, cns_results);
}
-void
-consensus_one_read_m4_nanopore(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid)
-{
- PackedDB& reads = *ctd->reads;
- ExtensionCandidate* overlaps = ctd->candidates;
- DiffRunningData* drd_s = ctd->drd_s;
- DiffRunningData* drd = NULL;
- M5Record* m5 = ctd->m5;
- CnsAlns& cns_vec = ctd->cns_alns;
- std::vector& cns_results = ctd->cns_results;
- const index_t read_size = overlaps[sid].ssize;
- std::vector& qstr = ctd->query;
- std::vector& tstr = ctd->target;
- tstr.resize(read_size);
- reads.GetSequence(read_id, true, tstr.data(), read_size);
- std::string& nqstr = ctd->qaln;
- std::string& ntstr = ctd->saln;
- const int min_align_size = ctd->rco.min_align_size;
- const double min_mapping_ratio = ctd->rco.min_mapping_ratio - 0.02;
+void consensus_one_read_m4_nanopore(ConsensusThreadData& ctd, ConsensusPerThreadData &pctd, const int64_t read_id, const int64_t sid, const int64_t eid) {
+ PackedDB& reads = ctd.reads;
+ ExtensionCandidate* overlaps = (ExtensionCandidate*)pctd.candidates;
+ DiffRunningData& drd = pctd.drd;
+ M5Record& m5 = pctd.m5;
+ CnsAlns& cns_vec = pctd.cns_alns;
+ std::vector& cns_results = pctd.cns_results;
+ const int64_t read_size = overlaps[read_id].ssize;
+ std::string& qstr = pctd.query;
+ std::string& tstr = pctd.target;
+ reads.GetSequence(read_id, 1, tstr);
+ std::string& nqstr = pctd.qaln;
+ std::string& ntstr = pctd.saln;
+ const int min_align_size = ctd.rco.min_align_size;
+ const double error_rate(ctd.rco.error_rate);
+ const double min_mapping_ratio = ctd.rco.min_mapping_ratio - 0.02;
- index_t L, R;
+ int64_t L, R;
if (eid - sid <= MAX_CNS_OVLPS)
{
L = sid;
@@ -332,183 +359,152 @@ consensus_one_read_m4_nanopore(ConsensusThreadData* ctd, const index_t read_id,
std::sort(overlaps + sid, overlaps + eid, CompareOverlapByOverlapSize());
}
- CnsTableItem* cns_table = ctd->cns_table;
- std::for_each(cns_table, cns_table + read_size, CnsTableItemCleaner());
+ std::vector& cns_table = pctd.cns_table;
+ cns_table.assign(read_size, CnsTableItem()); // reset table
cns_vec.clear();
- for (index_t i = L; i < R; ++i)
- {
- Overlap& ovlp = overlaps[i];
- qstr.resize(ovlp.qsize);
- reads.GetSequence(ovlp.qid, ovlp.qdir == FWD, qstr.data(), ovlp.qsize);
- index_t qext = ovlp.qext;
- index_t sext = ovlp.sext;
- if (ovlp.qdir == REV) qext = ovlp.qsize - 1 - qext;
- drd = drd_s;
- bool r = GetAlignment(qstr.data(), qext, qstr.size(), tstr.data(), sext, tstr.size(), drd, *m5, 0.20, min_align_size);
- if (r && check_ovlp_mapping_range(m5qoff(*m5), m5qend(*m5), ovlp.qsize, m5soff(*m5), m5send(*m5), ovlp.ssize, min_mapping_ratio))
- {
- normalize_gaps(m5qaln(*m5), m5saln(*m5), strlen(m5qaln(*m5)), nqstr, ntstr, true);
- meap_add_one_aln(nqstr, ntstr, m5soff(*m5), cns_table, tstr.data());
- cns_vec.add_aln(m5soff(*m5), m5send(*m5), nqstr, ntstr);
+ for (int64_t i(L); i < R; ++i) {
+ Overlap& ovlp(overlaps[i]);
+ reads.GetSequence(ovlp.qid, ovlp.qdir == FWD, qstr);
+ const int64_t qext(ovlp.qdir == FWD ? ovlp.qext : ovlp.qsize - 1 - qext);
+ const int r(GetAlignment(qstr, qext, tstr, ovlp.sext, drd, m5, error_rate, min_align_size));
+ if (r && check_ovlp_mapping_range(m5.qoff, m5.qend, ovlp.qsize, m5.soff, m5.send, ovlp.ssize, min_mapping_ratio)) {
+ normalize_gaps(m5.qaln, m5.saln, nqstr, ntstr, 1);
+ meap_add_one_aln(nqstr, ntstr, m5.soff, cns_table);
+ cns_vec.add_aln(m5.soff, m5.send, nqstr, ntstr);
}
}
std::vector mranges, eranges;
eranges.push_back(MappingRange(0, read_size));
- consensus_worker(cns_table, ctd->id_list, cns_vec, nqstr, ntstr, eranges, ctd->rco.min_cov, ctd->rco.min_size, read_id, cns_results);
+ consensus_worker(cns_table, pctd.id_list, cns_vec, nqstr, ntstr, eranges, ctd.rco.min_cov, ctd.rco.min_size, read_id, cns_results);
}
-struct CmpExtensionCandidateByScore
-{
- bool operator()(const ExtensionCandidate& a, const ExtensionCandidate& b)
- {
- if (a.score != b.score) return a.score > b.score;
- if (a.qid != b.qid) return a.qid < b.qid;
- return a.qext < b.qext;
- }
-};
+// limit coverage - if there's not yet sufficient coverage, add coverage and return true;
+// once coverage gets high enough, stop adding coverage and return false
-inline bool
-check_cov_stats(u1_t* cov_stats, int soff, int send)
-{
- const int max_cov = 20;
- int n = 0;
- for (int i = soff; i < send; ++i)
- if (cov_stats[i] >= max_cov)
+static inline int check_cov_stats(std::vector& cov_stats, const int soff, const int send) {
+ if (cov_stats.size() < static_cast(send)) {
+ cov_stats.resize(send, 0);
+ }
+ int n(0);
+ for (int i(soff); i < send; ++i) {
+ if (cov_stats[i] >= 20) { // max coverage
++n;
- if (send - soff >= n + 200)
- {
- for (int i = soff; i < send; ++i) ++cov_stats[i];
- return true;
+ }
}
- return false;
+ if (send - soff >= n + 200) {
+ for (int i(soff); i < send; ++i) {
+ // don't let small redundant region cause overflow
+ if (cov_stats[i] != std::numeric_limits::max()) {
+ ++cov_stats[i];
+ }
+ }
+ return 1;
+ }
+ return 0;
}
-void
-consensus_one_read_can_pacbio(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid)
-{
- PackedDB& reads = *ctd->reads;
- ExtensionCandidate* candidates = ctd->candidates;
- DiffRunningData* drd_s = ctd->drd_s;
- DiffRunningData* drd = NULL;
- M5Record* m5 = ctd->m5;
- CnsAlns& cns_vec = ctd->cns_alns;
- std::vector& cns_results = ctd->cns_results;
- const index_t read_size = candidates[sid].ssize;
- std::vector& qstr = ctd->query;
- std::vector& tstr = ctd->target;
- tstr.resize(read_size);
- reads.GetSequence(read_id, true, tstr.data(), read_size);
- std::string& nqstr = ctd->qaln;
- std::string& ntstr = ctd->saln;
- const int min_align_size = ctd->rco.min_align_size;
- const double min_mapping_ratio = ctd->rco.min_mapping_ratio - 0.02;
- const int max_added = 60;
-
- std::sort(candidates + sid, candidates + eid, CmpExtensionCandidateByScore());
- int num_added = 0;
- int num_ext = 0;
- const int max_ext = 200;
- CnsTableItem* cns_table = ctd->cns_table;
- std::for_each(cns_table, cns_table + read_size, CnsTableItemCleaner());
+void consensus_one_read_can_pacbio(ConsensusThreadData& ctd, ConsensusPerThreadData& pctd, const int64_t read_id, const int64_t sid, int64_t eid) {
+ const PackedDB& reads(ctd.reads);
+ ExtensionCandidateCompressed* candidates((ExtensionCandidateCompressed*)pctd.candidates);
+ DiffRunningData& drd(pctd.drd);
+ M5Record& m5(pctd.m5);
+ CnsAlns& cns_vec(pctd.cns_alns);
+ std::vector& cns_results(pctd.cns_results);
+ const int64_t ssize(reads.read_size(read_id));
+ std::string& qstr(pctd.query);
+ std::string& tstr(pctd.target);
+ reads.GetSequence(read_id, 1, tstr);
+ std::string& nqstr(pctd.qaln);
+ std::string& ntstr(pctd.saln);
+ const int min_align_size(ctd.rco.min_align_size);
+ const double error_rate(ctd.rco.error_rate);
+ const double min_mapping_ratio(ctd.rco.min_mapping_ratio - 0.02);
+ int num_added(0);
+ const int max_added(60);
+ eid = std::min(eid, sid + 200); // max of 200 extents
+ std::vector& cns_table(pctd.cns_table);
+ cns_table.assign(ssize, CnsTableItem()); // reset table
cns_vec.clear();
- std::set used_ids;
- u1_t* cov_stats = ctd->id_list;
- std::fill(cov_stats, cov_stats + MAX_SEQ_SIZE, 0);
- for (idx_t i = sid; i < eid && num_added < max_added && num_ext < max_ext; ++i)
- {
- ++num_ext;
- ExtensionCandidate& ec = candidates[i];
- r_assert(ec.sdir == FWD);
- if (used_ids.find(ec.qid) != used_ids.end()) continue;
- qstr.resize(ec.qsize);
- reads.GetSequence(ec.qid, ec.qdir == FWD, qstr.data(), ec.qsize);
- index_t qext = ec.qext;
- index_t sext = ec.sext;
- if (ec.qdir == REV) qext = ec.qsize - 1 - qext;
- drd = drd_s;
- bool r = GetAlignment(qstr.data(), qext, qstr.size(), tstr.data(), sext, tstr.size(), drd, *m5, 0.15, min_align_size);
- if (r && check_ovlp_mapping_range(m5qoff(*m5), m5qend(*m5), ec.qsize, m5soff(*m5), m5send(*m5), ec.ssize, min_mapping_ratio))
- {
- if (check_cov_stats(cov_stats, m5soff(*m5), m5send(*m5)))
- {
+ std::set used_ids;
+ std::vector& id_list(pctd.id_list); // used to be called cov_stats
+ id_list.clear();
+ for (int64_t i(sid); i < eid && num_added < max_added; ++i) {
+ const ExtensionCandidateCompressed& ec(candidates[i]);
+ if (used_ids.find(ec.qid) != used_ids.end()) {
+ continue;
+ }
+ const int64_t qsize(reads.read_size(ec.qid));
+ reads.GetSequence(ec.qid, ec.qdir() == FWD, qstr);
+ const int64_t qext(ec.qdir() == FWD ? ec.qext() : qsize - 1 - ec.qext());
+ const int r(GetAlignment(qstr, qext, tstr, ec.sext, drd, m5, error_rate, min_align_size));
+ if (r && check_ovlp_mapping_range(m5.qoff, m5.qend, qsize, m5.soff, m5.send, ssize, min_mapping_ratio)) {
+ if (check_cov_stats(id_list, m5.soff, m5.send)) {
++num_added;
used_ids.insert(ec.qid);
- normalize_gaps(m5qaln(*m5), m5saln(*m5), strlen(m5qaln(*m5)), nqstr, ntstr, true);
- meap_add_one_aln(nqstr, ntstr, m5soff(*m5), cns_table, tstr.data());
- cns_vec.add_aln(m5soff(*m5), m5send(*m5), nqstr, ntstr);
+ normalize_gaps(m5.qaln, m5.saln, nqstr, ntstr, 1);
+ meap_add_one_aln(nqstr, ntstr, m5.soff, cns_table);
+ cns_vec.add_aln(m5.soff, m5.send, nqstr, ntstr);
}
}
}
-
std::vector mranges, eranges;
cns_vec.get_mapping_ranges(mranges);
- get_effective_ranges(mranges, eranges, read_size, ctd->rco.min_size);
-
- consensus_worker(cns_table, ctd->id_list, cns_vec, nqstr, ntstr, eranges, ctd->rco.min_cov, ctd->rco.min_size, read_id, cns_results);
+ get_effective_ranges(mranges, eranges, ssize, ctd.rco.min_size);
+ if (ctd.rco.full_reads) {
+ consensus_worker_one_read(cns_table, pctd.id_list, cns_vec, nqstr, ntstr, eranges, ctd.rco.min_cov, ctd.rco.min_size, read_id, tstr, cns_results);
+ } else {
+ consensus_worker(cns_table, pctd.id_list, cns_vec, nqstr, ntstr, eranges, ctd.rco.min_cov, ctd.rco.min_size, read_id, cns_results);
+ }
}
-void
-consensus_one_read_can_nanopore(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid)
-{
- PackedDB& reads = *ctd->reads;
- ExtensionCandidate* candidates = ctd->candidates;
- DiffRunningData* drd_s = ctd->drd_s;
- DiffRunningData* drd = NULL;
- M5Record* m5 = ctd->m5;
- CnsAlns& cns_vec = ctd->cns_alns;
- std::vector& cns_results = ctd->cns_results;
- const index_t read_size = candidates[sid].ssize;
- std::vector& qstr = ctd->query;
- std::vector& tstr = ctd->target;
- tstr.resize(read_size);
- reads.GetSequence(read_id, true, tstr.data(), read_size);
- std::string& nqstr = ctd->qaln;
- std::string& ntstr = ctd->saln;
- const int min_align_size = ctd->rco.min_align_size;
- const double min_mapping_ratio = ctd->rco.min_mapping_ratio - 0.02;
-
- std::sort(candidates + sid, candidates + eid, CmpExtensionCandidateByScore());
- int num_added = 0;
- int num_ext = 0;
- const int max_ext = 200;
- CnsTableItem* cns_table = ctd->cns_table;
- std::for_each(cns_table, cns_table + read_size, CnsTableItemCleaner());
+void consensus_one_read_can_nanopore(ConsensusThreadData& ctd, ConsensusPerThreadData &pctd, const int64_t read_id, const int64_t sid, const int64_t eid) {
+ PackedDB& reads(ctd.reads);
+ ExtensionCandidateCompressed* candidates((ExtensionCandidateCompressed*)pctd.candidates);
+ DiffRunningData& drd(pctd.drd);
+ M5Record& m5(pctd.m5);
+ CnsAlns& cns_vec(pctd.cns_alns);
+ std::vector& cns_results(pctd.cns_results);
+ const int64_t ssize(reads.read_size(read_id));
+ std::string& qstr(pctd.query);
+ std::string& tstr(pctd.target);
+ reads.GetSequence(read_id, 1, tstr);
+ std::string& nqstr(pctd.qaln);
+ std::string& ntstr(pctd.saln);
+ const int min_align_size(ctd.rco.min_align_size);
+ const double error_rate(ctd.rco.error_rate);
+ const double min_mapping_ratio(ctd.rco.min_mapping_ratio - 0.02);
+ int num_added(0);
+ int num_ext(0);
+ const int max_ext(200);
+ std::vector& cns_table(pctd.cns_table);
+ cns_table.assign(ssize, CnsTableItem()); // reset table
cns_vec.clear();
std::set used_ids;
- u1_t* cov_stats = ctd->id_list;
- std::fill(cov_stats, cov_stats + MAX_SEQ_SIZE, 0);
- for (idx_t i = sid; i < eid && num_added < MAX_CNS_OVLPS && num_ext < max_ext; ++i)
- {
- ++num_ext;
- ExtensionCandidate& ec = candidates[i];
- r_assert(ec.sdir == FWD);
- if (used_ids.find(ec.qid) != used_ids.end()) continue;
- qstr.resize(ec.qsize);
- reads.GetSequence(ec.qid, ec.qdir == FWD, qstr.data(), ec.qsize);
- index_t qext = ec.qext;
- index_t sext = ec.sext;
- if (ec.qdir == REV) qext = ec.qsize - 1 - qext;
- drd = drd_s;
- bool r = GetAlignment(qstr.data(), qext, qstr.size(), tstr.data(), sext, tstr.size(), drd, *m5, 0.20, min_align_size);
- if (r && check_ovlp_mapping_range(m5qoff(*m5), m5qend(*m5), ec.qsize, m5soff(*m5), m5send(*m5), ec.ssize, min_mapping_ratio))
- {
- if (check_cov_stats(cov_stats, m5soff(*m5), m5send(*m5)))
- {
+ std::vector& id_list(pctd.id_list); // used to be called cov_stats
+ id_list.clear();
+ for (int64_t i(sid); i < eid && num_added < MAX_CNS_OVLPS && num_ext < max_ext; ++i) {
+ ++num_ext;
+ const ExtensionCandidateCompressed& ec(candidates[i]);
+ if (used_ids.find(ec.qid) != used_ids.end()) {
+ continue;
+ }
+ const int64_t qsize(reads.read_size(ec.qid));
+ reads.GetSequence(ec.qid, ec.qdir() == FWD, qstr);
+ const int64_t qext(ec.qdir() == FWD ? ec.qext() : qsize - 1 - ec.qext());
+ const int r(GetAlignment(qstr, qext, tstr, ec.sext, drd, m5, error_rate, min_align_size));
+ if (r && check_ovlp_mapping_range(m5.qoff, m5.qend, qsize, m5.soff, m5.send, ssize, min_mapping_ratio)) {
+ if (check_cov_stats(id_list, m5.soff, m5.send)) {
++num_added;
used_ids.insert(ec.qid);
- normalize_gaps(m5qaln(*m5), m5saln(*m5), strlen(m5qaln(*m5)), nqstr, ntstr, true);
- meap_add_one_aln(nqstr, ntstr, m5soff(*m5), cns_table, tstr.data());
- cns_vec.add_aln(m5soff(*m5), m5send(*m5), nqstr, ntstr);
+ normalize_gaps(m5.qaln, m5.saln, nqstr, ntstr, 1);
+ meap_add_one_aln(nqstr, ntstr, m5.soff, cns_table);
+ cns_vec.add_aln(m5.soff, m5.send, nqstr, ntstr);
}
}
}
-
std::vector mranges, eranges;
- eranges.push_back(MappingRange(0, read_size));
-
- consensus_worker(cns_table, ctd->id_list, cns_vec, nqstr, ntstr, eranges, ctd->rco.min_cov, ctd->rco.min_size, read_id, cns_results);
+ eranges.push_back(MappingRange(0, ssize));
+ consensus_worker(cns_table, pctd.id_list, cns_vec, nqstr, ntstr, eranges, ctd.rco.min_cov, ctd.rco.min_size, read_id, cns_results);
}
-
-} // namespace ns_meap_cns {
diff --git a/src/mecat2cns/mecat_correction.h b/src/mecat2cns/mecat_correction.h
index c00d463..176fcea 100644
--- a/src/mecat2cns/mecat_correction.h
+++ b/src/mecat2cns/mecat_correction.h
@@ -1,22 +1,16 @@
#ifndef MEAP_CORRECTION_H
#define MEAP_CORRECTION_H
-#include "reads_correction_aux.h"
+#include // int64_t
-namespace ns_meap_cns {
+#include "reads_correction_aux.h" // ConsensusPerThreadData, ConsensusThreadData
-void
-consensus_one_read_m4_pacbio(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid);
+void consensus_one_read_m4_pacbio(ConsensusThreadData& ctd, ConsensusPerThreadData &cptd, const int64_t read_id, const int64_t sid, const int64_t eid);
-void
-consensus_one_read_m4_nanopore(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid);
+void consensus_one_read_m4_nanopore(ConsensusThreadData& ctd, ConsensusPerThreadData &cptd, const int64_t read_id, const int64_t sid, const int64_t eid);
-void
-consensus_one_read_can_pacbio(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid);
+void consensus_one_read_can_pacbio(ConsensusThreadData& ctd, ConsensusPerThreadData &cptd, const int64_t read_id, const int64_t sid, const int64_t eid);
-void
-consensus_one_read_can_nanopore(ConsensusThreadData* ctd, const index_t read_id, const index_t sid, const index_t eid);
-
-} // namespace ns_meap_cns
+void consensus_one_read_can_nanopore(ConsensusThreadData& ctd, ConsensusPerThreadData &cptd, const int64_t read_id, const int64_t sid, const int64_t eid);
#endif // MEAP_CORRECTION_H
diff --git a/src/mecat2cns/options.cpp b/src/mecat2cns/options.cpp
index 1ddaee1..e7da475 100644
--- a/src/mecat2cns/options.cpp
+++ b/src/mecat2cns/options.cpp
@@ -4,31 +4,33 @@
#include
#include
+#include
+#include
-using namespace std;
+#include "../common/defs.h" // TECH_NANOPORE, TECH_PACBIO
-static int input_type_pacbio = 1;
-static int num_threads_pacbio = 1;
-static index_t batch_size_pacbio = 100000;
-static double mapping_ratio_pacbio = 0.9;
-static int align_size_pacbio = 2000;
-static int cov_pacbio = 6;
-static int min_size_pacbio = 5000;
-static bool print_usage_pacbio = false;
-static int tech_pacbio = TECH_PACBIO;
-static int num_partition_files = 10;
+static const double mapping_ratio_pacbio = 0.6;
+static const int align_size_pacbio = 1000;
+static const int cov_pacbio = 4;
+static const int min_size_pacbio = 2000;
+static const int tech_pacbio = TECH_PACBIO;
+static const double error_rate_pacbio = .15;
-static int input_type_nanopore = 1;
-static int num_threads_nanopore = 1;
-static index_t batch_size_nanopore = 100000;
-static double mapping_ratio_nanopore = 0.4;
-static int align_size_nanopore = 400;
-static int cov_nanopore = 6;
-static int min_size_nanopore = 2000;
-static bool print_usage_nanopore = false;
-static int tech_nanopore = TECH_NANOPORE;
+static const double mapping_ratio_nanopore = 0.4;
+static const int align_size_nanopore = 400;
+static const int cov_nanopore = 6;
+static const int min_size_nanopore = 2000;
+static const int tech_nanopore = TECH_NANOPORE;
+static const double error_rate_nanopore = .2;
-static int default_tech = TECH_PACBIO;
+static const int default_tech = TECH_PACBIO;
+static const int num_partition_files = 0;
+static const int full_reads = 0;
+static const size_t read_buffer_size = 0;
+static const size_t batch_size = size_t(1) << 33; // 8 GB
+static const int print_usage_info = 0;
+static const int input_type = 1;
+static const int num_threads = 1;
static const char input_type_n = 'i';
static const char num_threads_n = 't';
@@ -40,264 +42,306 @@ static const char min_size_n = 'l';
static const char usage_n = 'h';
static const char tech_n = 'x';
static const char num_partition_files_n = 'k';
+static const char grid_options_n = 'G';
+static const char grid_options_split_n = 'S';
+static const char job_index_n = 'I';
+static const char reads_to_correct_n = 'R';
+static const char grid_start_delay_n = 'D';
+static const char full_reads_n = 'F';
+static const char read_buffer_size_n = 'b';
-void
-print_pacbio_default_options()
-{
- cerr << '-' << input_type_n << ' ' << input_type_pacbio
- << ' '
- << '-' << num_threads_n << ' ' << num_threads_pacbio
- << ' '
- << '-' << batch_size_n << ' ' << batch_size_pacbio
- << ' '
- << '-' << mapping_ratio_n << ' ' << mapping_ratio_pacbio
- << ' '
- << '-' << align_size_n << ' ' << align_size_pacbio
- << ' '
- << '-' << cov_n << ' ' << cov_pacbio
- << ' '
- << '-' << min_size_n << ' ' << min_size_pacbio
- << ' '
- << '-' << num_partition_files_n << ' ' << num_partition_files
- << "\n";
+// convert number with potential suffix (k, m, g)
+static size_t convert_integer(const std::string& s) {
+ std::istringstream x(s);
+ size_t value;
+ x >> value;
+ const size_t i(s.find_first_not_of("0123456789"));
+ if (i != std::string::npos) {
+ switch (s[i]) {
+ case 'g':
+ value *= 1024;
+ case 'm':
+ value *= 1024;
+ case 'k':
+ value *= 1024;
+ }
+ }
+ return value;
}
-void print_nanopore_default_options()
-{
- cerr << '-' << input_type_n << ' ' << input_type_nanopore
- << ' '
- << '-' << num_threads_n << ' ' << num_threads_nanopore
- << ' '
- << '-' << batch_size_n << ' ' << batch_size_nanopore
- << ' '
- << '-' << mapping_ratio_n << ' ' << mapping_ratio_nanopore
- << ' '
- << '-' << align_size_n << ' ' << align_size_nanopore
- << ' '
- << '-' << cov_n << ' ' << cov_nanopore
- << ' '
- << '-' << min_size_n << ' ' << min_size_nanopore
- << ' '
- << '-' << num_partition_files_n << ' ' << num_partition_files
+static void print_pacbio_default_options() {
+ std::cerr << " -" << mapping_ratio_n << " " << mapping_ratio_pacbio
+ << " -" << align_size_n << " " << align_size_pacbio
+ << " -" << cov_n << " " << cov_pacbio << " "
+ << " -" << min_size_n << " " << min_size_pacbio
<< "\n";
}
-void
-print_usage(const char* prog)
-{
- cerr << "USAGE:\n"
- << prog << ' ' << "[options]" << ' ' << "input" << ' ' << "reads" << ' ' << "output" << "\n";
- cerr << "\n" << "OPTIONS:" << "\n";
-
- cerr << "-" << tech_n << " <0/1>\t" << "sequencing platform: 0 = PACBIO, 1 = NANOPORE" << "\n"
- << "\t\t" << "default: 0" << "\n";
-
- cerr << "-" << input_type_n << " <0/1>\t" << "input type: 0 = candidte, 1 = m4" << "\n";
-
- cerr << "-" << num_threads_n << " \t" << "number of threads (CPUs)" << "\n";
-
- cerr << "-" << batch_size_n << " \t" << "batch size that the reads will be partitioned" << "\n";
-
- cerr << "-" << mapping_ratio_n << " \t" << "minimum mapping ratio" << "\n";
-
- cerr << "-" << align_size_n << " \t" << "minimum overlap size" << "\n";
-
- cerr << "-" << cov_n << " \t" << "minimum coverage under consideration" << "\n";
-
- cerr << "-" << min_size_n << " \t" << "minimum length of corrected sequence" << "\n";
-
- cerr << "-" << num_partition_files_n << " \t"
- << "number of partition files when partitioning overlap results"
- << " (if < 0, then it will be set to system limit value)"
+static void print_nanopore_default_options() {
+ std::cerr << " -" << mapping_ratio_n << " " << mapping_ratio_nanopore
+ << " -" << align_size_n << " " << align_size_nanopore
+ << " -" << cov_n << " " << cov_nanopore
+ << " -" << min_size_n << " " << min_size_nanopore
<< "\n";
-
- cerr << "-" << usage_n << "\t\t" << "print usage info." << "\n";
-
- cerr << "\n"
- << "If 'x' is set to be '0' (pacbio), then the other options have the following default values: \n";
+}
+
+// given options, recreate arguments from the command line
+void make_options(const ReadsCorrectionOptions& options, std::string& new_options) {
+ std::ostringstream cmd;
+ cmd << " -" << input_type_n << " " << (options.input_type == INPUT_TYPE_CAN ? 0 : 1);
+ if (options.num_threads > -1) {
+ cmd << " -" << num_threads_n << " " << options.num_threads;
+ }
+ if (options.batch_size && options.batch_size != batch_size) {
+ cmd << " -" << batch_size_n << " " << options.batch_size;
+ }
+ if (options.min_mapping_ratio >= 0 && ((options.tech == TECH_PACBIO && options.min_mapping_ratio != mapping_ratio_pacbio) || (options.tech != TECH_PACBIO && options.min_mapping_ratio != mapping_ratio_nanopore))) {
+ cmd << " -" << mapping_ratio_n << " " << options.min_mapping_ratio;
+ }
+ if (options.min_align_size >= 0 && ((options.tech == TECH_PACBIO && options.min_align_size != align_size_pacbio) || (options.tech != TECH_PACBIO && options.min_align_size != align_size_nanopore))) {
+ cmd << " -" << align_size_n << " " << options.min_align_size;
+ }
+ if (options.min_cov >= 0 && ((options.tech == TECH_PACBIO && options.min_cov != cov_pacbio) || (options.tech != TECH_PACBIO && options.min_cov != cov_nanopore))) {
+ cmd << " -" << cov_n << " " << options.min_cov;
+ }
+ if (options.min_size >= 0 && ((options.tech == TECH_PACBIO && options.min_size != min_size_pacbio) || (options.tech != TECH_PACBIO && options.min_size != min_size_nanopore))) {
+ cmd << " -" << min_size_n << " " << options.min_size;
+ }
+ if (options.num_partition_files > 0 && options.num_partition_files != num_partition_files) {
+ cmd << " -" << num_partition_files_n << " " << options.num_partition_files;
+ }
+ if (options.grid_options != NULL) {
+ cmd << " -" << grid_options_n << " \"" << options.grid_options << "\"";
+ }
+ if (options.grid_options_split != NULL) {
+ cmd << " -" << grid_options_split_n << " \"" << options.grid_options_split << "\"";
+ }
+ if (options.job_index != -1) {
+ cmd << " -" << job_index_n << " " << options.job_index;
+ }
+ if (options.reads_to_correct) {
+ cmd << " -" << reads_to_correct_n << " " << options.reads_to_correct;
+ }
+ if (options.grid_start_delay) {
+ cmd << " -" << grid_start_delay_n << " " << options.grid_start_delay;
+ }
+ if (options.full_reads) {
+ cmd << " -" << full_reads_n;
+ }
+ if (options.read_buffer_size) {
+ cmd << " -" << read_buffer_size_n << " " << options.read_buffer_size;
+ }
+ cmd << " " << options.m4;
+ cmd << " " << options.reads;
+ cmd << " " << options.corrected_reads;
+ new_options = cmd.str();
+}
+
+void print_usage(const char* prog) {
+ std::cerr << "USAGE:\n"
+ << prog << " [options] input reads output\n"
+ << "\n"
+ << "OPTIONS:\n"
+ << "-" << tech_n << " <0/1>\tsequencing platform: 0 = PACBIO, 1 = NANOPORE [0]\n"
+ << "-" << input_type_n << " <0/1>\tinput type: 0 = candidate, 1 = m4\n"
+ << "-" << num_threads_n << " \tnumber of threads (CPUs)\n"
+ << "-" << batch_size_n << " \tmemory used for holding candidates [8 GB]\n"
+ << "-" << mapping_ratio_n << " \tminimum mapping ratio\n"
+ << "-" << align_size_n << " \tminimum overlap size\n"
+ << "-" << cov_n << " \tminimum coverage under consideration\n"
+ << "-" << min_size_n << " \tminimum length of corrected sequence\n"
+ << "-" << num_partition_files_n << " \tnumber of partition files per pass when partitioning overlap results (if 0, then use system limit)\n"
+ << "-" << grid_options_n << " \toptions for grid submission\n"
+ << "-" << grid_options_split_n << " \toptions for split grid submission\n"
+ << "-" << reads_to_correct_n << " \tnumber of reads to correct [all]\n"
+ << "-" << grid_start_delay_n << " \tseconds to delay between starting grid jobs\n"
+ << "-" << full_reads_n << "\t\toutput full reads, not just the corrected parts\n"
+ << "-" << read_buffer_size_n << " \tbytes of memory to buffer reads [no buffer]\n"
+ << "-" << usage_n << "\t\tprint usage info\n"
+ << "\n"
+ << "If 'x' is set to be '0' (pacbio), then the other options have the following default values: \n";
print_pacbio_default_options();
-
- cerr << "\n"
- << "If 'x' is set to be '1' (nanopore), then the other options have the following default values: \n";
+ std::cerr << "\n"
+ << "If 'x' is set to be '1' (nanopore), then the other options have the following default values: \n";
print_nanopore_default_options();
}
-ConsensusOptions
-init_consensus_options(int tech)
-{
- ConsensusOptions t;
+ReadsCorrectionOptions init_consensus_options(const int tech) {
+ ReadsCorrectionOptions t;
+ t.m4 = NULL;
+ t.reads = NULL;
+ t.corrected_reads = NULL;
+ t.grid_options = NULL;
+ t.grid_options_split = NULL;
+ t.num_partition_files = num_partition_files;
+ t.job_index = -1;
+ t.reads_to_correct = 0;
+ t.grid_start_delay = 0;
+ t.full_reads = full_reads;
+ t.read_buffer_size = read_buffer_size;
+ t.batch_size = batch_size;
+ t.print_usage_info = print_usage_info;
+ t.input_type = input_type;
+ t.num_threads = num_threads;
if (tech == TECH_PACBIO) {
- t.input_type = input_type_pacbio;
- t.m4 = NULL;
- t.reads = NULL;
- t.corrected_reads = NULL;
- t.num_threads = num_threads_pacbio;
- t.batch_size = batch_size_pacbio;
t.min_mapping_ratio = mapping_ratio_pacbio;
t.min_align_size = align_size_pacbio;
t.min_cov = cov_pacbio;
t.min_size = min_size_pacbio;
- t.print_usage_info = print_usage_pacbio;
- t.num_partition_files = num_partition_files;
t.tech = tech_pacbio;
+ t.error_rate = error_rate_pacbio;
} else {
- t.input_type = input_type_nanopore;
- t.m4 = NULL;
- t.reads = NULL;
- t.corrected_reads = NULL;
- t.num_threads = num_threads_nanopore;
- t.batch_size = batch_size_nanopore;
t.min_mapping_ratio = mapping_ratio_nanopore;
t.min_align_size = align_size_nanopore;
t.min_cov = cov_nanopore;
t.min_size = min_size_nanopore;
- t.print_usage_info = print_usage_nanopore;
- t.num_partition_files = num_partition_files;
t.tech = tech_nanopore;
+ t.error_rate = error_rate_nanopore;
}
- return t;
+ return t;
}
-int detect_tech(int argc, char* argv[])
-{
- int t = default_tech;
- char tech_nstr[64]; tech_nstr[0] = '-'; tech_nstr[1] = tech_n; tech_nstr[2] = '\0';
-
- for (int i = 0; i < argc; ++i) {
- if (strcmp(tech_nstr, argv[i]) == 0) {
- if (i + 1 == argc) {
- fprintf(stderr, "argument to option '%c' is missing.\n", tech_n);
- t = -1;
- }
- cout << tech_nstr << "\n";
- if (argv[i + 1][0] == '0') {
- t = TECH_PACBIO;
- } else if (argv[i + 1][0] == '1') {
- t = TECH_NANOPORE;
- } else {
- fprintf(stderr, "invalid argument to option '%c': %s\n", tech_n, argv[i + 1]);
- t = -1;
- }
- break;
- }
- }
-
- return t;
+int detect_tech(int argc, char* argv[]) {
+ int t = default_tech;
+ const char tech_nstr[] = {'-', tech_n, 0};
+ for (int i = 0; i < argc; ++i) {
+ if (strcmp(tech_nstr, argv[i]) == 0) {
+ if (i + 1 == argc) {
+ std::cerr << "argument to option '" << tech_n << "' is missing.\n";
+ t = -1;
+ }
+ std::cout << tech_nstr << "\n";
+ if (argv[i + 1][0] == '0') {
+ t = TECH_PACBIO;
+ } else if (argv[i + 1][0] == '1') {
+ t = TECH_NANOPORE;
+ } else {
+ std::cerr << "invalid argument to option '" << tech_n << "': " << argv[i + 1] << "\n";
+ t = -1;
+ }
+ break;
+ }
+ }
+ return t;
}
-int
-parse_arguments(int argc, char* argv[], ConsensusOptions& t)
-{
- bool parse_success = true;
- int tech = detect_tech(argc, argv);
+int parse_arguments(int argc, char* argv[], ReadsCorrectionOptions& t) {
+ int parse_success(1);
+ const int tech(detect_tech(argc, argv));
if (tech == -1) {
return 1;
}
t = init_consensus_options(tech);
-
int opt_char;
- char err_char;
- opterr = 0;
- while((opt_char = getopt(argc, argv, "i:t:p:r:a:c:l:x:k:h")) != -1) {
+ opterr = 0;
+ while ((opt_char = getopt(argc, argv, "i:t:p:r:a:c:l:x:hG:I:n:R:S:D:k:Fb:")) != -1) {
switch (opt_char) {
case input_type_n:
- if (optarg[0] == '0')
+ if (optarg[0] == '0') {
t.input_type = INPUT_TYPE_CAN;
- else if (optarg[0] == '1')
+ } else if (optarg[0] == '1') {
t.input_type = INPUT_TYPE_M4;
- else {
- fprintf(stderr, "invalid argument to option '%c': %s\n", input_type_n, optarg);
+ } else {
+ std::cerr << "invalid argument to option '" << input_type_n << "': " << optarg << "\n";
return 1;
}
break;
case num_threads_n:
- t.num_threads = atoi(optarg);
+ t.num_threads = convert_integer(optarg);
break;
case batch_size_n:
- t.batch_size = atoll(optarg);
+ t.batch_size = convert_integer(optarg);
break;
case mapping_ratio_n:
t.min_mapping_ratio = atof(optarg);
break;
case align_size_n:
- t.min_align_size = atoi(optarg);
+ t.min_align_size = convert_integer(optarg);
break;
case cov_n:
- t.min_cov = atoi(optarg);
+ t.min_cov = convert_integer(optarg);
break;
case min_size_n:
- t.min_size = atoll(optarg);
+ t.min_size = convert_integer(optarg);
break;
case usage_n:
- t.print_usage_info = true;
+ t.print_usage_info = 1;
+ break;
+ case grid_options_n:
+ t.grid_options = optarg;
+ break;
+ case grid_options_split_n:
+ t.grid_options_split = optarg;
+ break;
+ case job_index_n:
+ t.job_index = convert_integer(optarg);
+ break;
+ case reads_to_correct_n:
+ t.reads_to_correct = convert_integer(optarg);
+ break;
+ case grid_start_delay_n:
+ t.grid_start_delay = convert_integer(optarg);
break;
case tech_n:
break;
case num_partition_files_n:
- t.num_partition_files = atoi(optarg);
+ t.num_partition_files = convert_integer(optarg);
+ break;
+ case full_reads_n:
+ t.full_reads = 1;
+ break;
+ case read_buffer_size_n:
+ t.read_buffer_size = convert_integer(optarg);
break;
case '?':
- err_char = (char)optopt;
- fprintf(stderr, "unrecognised option '%c'\n", err_char);
- return 1;
- break;
- case ':':
- err_char = (char)optopt;
- fprintf(stderr, "argument to option '%c' is missing.\n", err_char);
- return 1;
- break;
+ std::cerr << "unrecognised option '" << char(optopt) << "'\n";
+ return 1;
+ break;
+ case ':':
+ std::cerr << "argument to option '" << char(optopt) << "' is missing.\n";
+ return 1;
+ break;
}
}
-
- if (t.num_threads <= 0)
- {
- std::cerr << "cpu threads must be greater than 0\n";
- parse_success = false;
- }
- if (t.batch_size <= 0)
- {
+ if (t.batch_size == 0) {
std::cerr << "batch size must be greater than 0\n";
- parse_success = false;
+ parse_success = 0;
}
- if (t.min_mapping_ratio < 0.0)
- {
+ if (t.min_mapping_ratio < 0.0) {
std::cerr << "mapping ratio must be >= 0.0\n";
- parse_success = false;
- }
- if (t.min_cov < 0)
- {
- std::cerr << "coverage must be >= 0\n";
- parse_success = false;
+ parse_success = 0;
}
- if (t.min_cov < 0)
- {
+ if (t.min_size < 0) {
std::cerr << "sequence size must be >= 0\n";
- parse_success = false;
+ parse_success = 0;
+ }
+ if (argc - optind < 3) {
+ return 1;
}
-
- if (argc < 3) return 1;
-
t.m4 = argv[argc - 3];
t.reads = argv[argc - 2];
t.corrected_reads = argv[argc - 1];
-
- if (parse_success) return 0;
- return 1;
+ return parse_success;
}
-void
-print_options(ConsensusOptions& t)
-{
- cout << "input_type:\t" << t.input_type << "\n";
- if (t.m4) cout << "reads\t" << t.m4 << "\n";
- if (t.reads) cout << "output\t" << t.reads << "\n";
- if (t.corrected_reads) cout << "m4\t" << t.corrected_reads << "\n";
- cout << "number of threads:\t" << t.num_threads << "\n";
- cout << "batch size:\t" << t.batch_size << "\n";
- cout << "mapping ratio:\t" << t.min_mapping_ratio << "\n";
- cout << "align size:\t" << t.min_align_size << "\n";
- cout << "cov:\t" << t.min_cov << "\n";
- cout << "min size:\t" << t.min_size << "\n";
- cout << "partition files:\t" << t.num_partition_files << "\n";
- cout << "tech:\t" << t.tech << "\n";
+// not currently used
+#if 0
+static void print_options(ReadsCorrectionOptions& t) {
+ std::cout << "input_type:\t" << t.input_type << "\n";
+ if (t.m4) std::cout << "reads\t" << t.m4 << "\n";
+ if (t.reads) std::cout << "output\t" << t.reads << "\n";
+ if (t.corrected_reads) std::cout << "m4\t" << t.corrected_reads << "\n";
+ if (t.grid_options) std::cout << "grid\t" << t.grid_options << "\n";
+ if (t.grid_options_split) std::cout << "grid_split\t" << t.grid_options_split << "\n";
+ if (t.full_reads) std::cout << "full reads\n";
+ if (t.read_buffer_size) std::cout << "read_buffer_size\t" << t.read_buffer_size << "\n";
+ std::cout << "number of threads:\t" << t.num_threads << "\n";
+ std::cout << "batch size:\t" << t.batch_size << "\n";
+ std::cout << "mapping ratio:\t" << t.min_mapping_ratio << "\n";
+ std::cout << "align size:\t" << t.min_align_size << "\n";
+ std::cout << "cov:\t" << t.min_cov << "\n";
+ std::cout << "min size:\t" << t.min_size << "\n";
+ std::cout << "partition files:\t" << t.num_partition_files << "\n";
+ std::cout << "tech:\t" << t.tech << "\n";
}
+#endif
diff --git a/src/mecat2cns/options.h b/src/mecat2cns/options.h
index 6dac71c..489b91b 100644
--- a/src/mecat2cns/options.h
+++ b/src/mecat2cns/options.h
@@ -1,37 +1,39 @@
#ifndef OPTIONS_H
#define OPTIONS_H
-#include "../common/defs.h"
+#include // string
#define INPUT_TYPE_CAN 0
#define INPUT_TYPE_M4 1
-struct ConsensusOptions
-{
+struct ReadsCorrectionOptions {
int input_type;
const char* m4;
const char* reads;
const char* corrected_reads;
+ const char* grid_options;
+ const char* grid_options_split;
int num_threads;
- index_t batch_size;
+ size_t batch_size;
double min_mapping_ratio;
int min_align_size;
int min_cov;
- index_t min_size;
- bool print_usage_info;
+ int min_size;
+ int print_usage_info;
int tech;
- int num_partition_files;
+ int num_partition_files;
+ int job_index;
+ int reads_to_correct;
+ int grid_start_delay;
+ int full_reads;
+ size_t read_buffer_size;
+ double error_rate; // .15 for pacbio, .2 for nanopore
};
-void
-print_usage(const char* prog);
+void print_usage(const char* prog);
-int
-parse_arguments(int argc, char* argv[], ConsensusOptions& t);
+int parse_arguments(int argc, char* argv[], ReadsCorrectionOptions& t);
-void
-print_options(ConsensusOptions& t);
-
-typedef ConsensusOptions ReadsCorrectionOptions;
+void make_options(const ReadsCorrectionOptions& t, std::string& new_options_str);
#endif // OPTIONS_H
diff --git a/src/mecat2cns/overlaps_partition.cpp b/src/mecat2cns/overlaps_partition.cpp
index 9698265..ec560a5 100644
--- a/src/mecat2cns/overlaps_partition.cpp
+++ b/src/mecat2cns/overlaps_partition.cpp
@@ -1,424 +1,247 @@
#include "overlaps_partition.h"
+#include // assert()
#include
#include
#include
+#include // stat(), struct stat
#include
-
-#include
+#include // make_pair(), pair<>
+#include // ceil()
#include "overlaps_store.h"
#include "reads_correction_aux.h"
+#include "packed_db.h" // PackedDB
-using namespace std;
-
-#define error_and_exit(msg) { std::cerr << msg << "\n"; abort(); }
-
-inline bool
-check_m4record_mapping_range(const M4Record& m4, const double min_cov_ratio)
-{
- const index_t qm = m4qend(m4) - m4qoff(m4);
- const index_t qs = m4qsize(m4) * min_cov_ratio;
- const index_t sm = m4send(m4) - m4soff(m4);
- const index_t ss = m4ssize(m4) * min_cov_ratio;
- return qm >= qs || sm >= ss;
+inline static bool query_is_contained(const M4Record& m4, const double min_cov_ratio) {
+ return m4qend(m4) - m4qoff(m4) >= m4qsize(m4) * min_cov_ratio;
}
-inline bool
-query_is_contained(const M4Record& m4, const double min_cov_ratio)
-{
- const index_t qm = m4qend(m4) - m4qoff(m4);
- const index_t qs = m4qsize(m4) * min_cov_ratio;
- return qm >= qs;
+inline static bool subject_is_contained(const M4Record& m4, const double min_cov_ratio) {
+ return m4send(m4) - m4soff(m4) >= m4ssize(m4) * min_cov_ratio;
}
-inline bool
-subject_is_contained(const M4Record& m4, const double min_cov_ratio)
-{
- const index_t sm = m4send(m4) - m4soff(m4);
- const index_t ss = m4ssize(m4) * min_cov_ratio;
- return sm >= ss;
+inline static bool check_m4record_mapping_range(const M4Record& m4, const double min_cov_ratio) {
+ return query_is_contained(m4, min_cov_ratio) || subject_is_contained(m4, min_cov_ratio);
}
-void
-get_qualified_m4record_counts(const char* m4_file_name, const double min_cov_ratio, index_t& num_qualified_records, index_t& num_reads)
-{
- std::ifstream in;
- open_fstream(in, m4_file_name, std::ios::in);
- num_qualified_records = 0;
- num_reads = -1;
- index_t num_records = 0;
- M4Record m4;
+static int64_t get_qualified_m4record_counts(const char* const m4_file_name, const double min_cov_ratio) {
+ std::ifstream in;
+ open_fstream(in, m4_file_name, std::ios::in);
+ int64_t num_reads(-1), num_records(0), num_qualified_records(0);
+ M4Record m4;
m4qext(m4) = m4sext(m4) = INVALID_IDX;
- while (in >> m4)
- {
- if (m4qext(m4) == INVALID_IDX || m4sext(m4) == INVALID_IDX)
- {
+ while (in >> m4) {
+ if (m4qext(m4) == INVALID_IDX || m4sext(m4) == INVALID_IDX) {
ERROR("no gapped start position is provided, please make sure that you have run \'meap_pairwise\' with option \'-g 1\'");
}
- ++num_records;
- if (check_m4record_mapping_range(m4, min_cov_ratio)) ++num_qualified_records;
- num_reads = std::max(num_reads, m4qid(m4));
- num_reads = std::max(num_reads, m4sid(m4));
- }
- close_fstream(in);
-
- LOG(stderr, "there are %d overlaps, %d are qualified.", (int)num_records, (int)num_qualified_records);
- ++num_reads;
+ ++num_records;
+ if (check_m4record_mapping_range(m4, min_cov_ratio)) {
+ ++num_qualified_records;
+ }
+ num_reads = std::max(num_reads, m4qid(m4));
+ num_reads = std::max(num_reads, m4sid(m4));
+ }
+ close_fstream(in);
+ LOG(stderr, "there are %ld overlaps, %ld are qualified.", num_records, num_qualified_records);
+ return num_reads + 1;
}
-void
-get_repeat_reads(const char* m4_file_name, const double min_cov_ratio, const index_t num_reads, set& repeat_reads)
-{
- const char MaxContained = 100;
- char cnt_table[MaxContained + 1];
- for (int i = 0; i < MaxContained; ++i) cnt_table[i] = i + 1;
- cnt_table[MaxContained] = MaxContained;
- vector cnts(num_reads, 0);
+// not in use at the moment
+#if 0
+static void get_repeat_reads(const char* const m4_file_name, const double min_cov_ratio, const int64_t num_reads, std::set& repeat_reads) {
+ const int max_contained(100);
+ // used to increment to a max of max_contained
+ char count_table[max_contained + 1];
+ for (int i(0); i < max_contained; ++i) {
+ count_table[i] = i + 1;
+ }
+ count_table[max_contained] = max_contained;
+ std::vector counts(num_reads, 0);
std::ifstream in;
- open_fstream(in, m4_file_name, std::ios::in);
+ open_fstream(in, m4_file_name, std::ios::in);
M4Record m4;
m4qext(m4) = m4sext(m4) = INVALID_IDX;
- while (in >> m4)
- {
- if (query_is_contained(m4, min_cov_ratio))
- {
- index_t qid = m4qid(m4);
- cnts[qid] = cnt_table[cnts[qid]];
+ while (in >> m4) {
+ if (query_is_contained(m4, min_cov_ratio)) {
+ const int64_t qid = m4qid(m4);
+ // increments count up to max_contained
+ counts[qid] = count_table[static_cast(counts[qid])];
}
- if (subject_is_contained(m4, min_cov_ratio))
- {
- index_t sid = m4sid(m4);
- cnts[sid] = cnt_table[cnts[sid]];
+ if (subject_is_contained(m4, min_cov_ratio)) {
+ const int64_t sid = m4sid(m4);
+ // increments count up to max_contained
+ counts[sid] = count_table[static_cast(counts[sid])];
}
- }
- close_fstream(in);
-
- for(index_t i = 0; i < num_reads; ++i)
- if (cnts[i] >= MaxContained)
- {
- cerr << "repeat read " << i << "\n";
+ }
+ close_fstream(in);
+ for (int64_t i(0); i < num_reads; ++i) {
+ if (counts[i] == max_contained) {
+ std::cerr << "repeat read " << i << "\n";
repeat_reads.insert(i);
}
-
- LOG(stderr, "number of repeat reads: %d", (int)repeat_reads.size());
-}
-
-void
-generate_partition_index_file_name(const char* m4_file_name, std::string& ret)
-{
- ret = m4_file_name;
- ret += ".partition_files";
-}
-
-void
-generate_partition_file_name(const char* m4_file_name, const index_t part, std::string& ret)
-{
- ret = m4_file_name;
- ret += ".part";
- std::ostringstream os;
- os << part;
- ret += os.str();
-}
-
-idx_t
-get_num_reads(const char* candidates_file)
-{
- ifstream in;
- open_fstream(in, candidates_file, ios::in);
- ExtensionCandidate ec;
- int max_id = -1;
- while (in >> ec)
- {
- max_id = std::max(ec.qid, max_id);
- max_id = std::max(ec.sid, max_id);
}
- close_fstream(in);
- return max_id + 1;
+ LOG(stderr, "number of repeat reads: %lu", repeat_reads.size());
}
+#endif
-void
-normalise_candidate(ExtensionCandidate& src, ExtensionCandidate& dst, const bool subject_is_target)
-{
- if (subject_is_target)
- {
- dst = src;
- }
- else
- {
- dst.qdir = src.sdir;
- dst.qid = src.sid;
- dst.qext = src.sext;
- dst.qsize = src.ssize;
- dst.sdir = src.qdir;
- dst.sid = src.qid;
- dst.sext = src.qext;
- dst.ssize = src.qsize;
- dst.score = src.score;
- }
-
- if (dst.sdir == REV)
- {
- dst.qdir = REVERSE_STRAND(dst.qdir);
- dst.sdir = REVERSE_STRAND(dst.sdir);
- }
+void generate_partition_index_file_name(const std::string& input_file_name, std::string& ret) {
+ ret = input_file_name + ".partition_files";
}
-int
-fix_file_counts(int num_files) {
- if (num_files < 0) {
- num_files = sysconf(_SC_OPEN_MAX) - 10;
- }
- return num_files;
+static void generate_partition_file_name(const std::string& input_file_name, const int part, std::string& ret) {
+ std::ostringstream os;
+ os << part;
+ ret = input_file_name + ".part" + os.str();
}
-void
-partition_candidates(const char* input, const idx_t batch_size, const int min_read_size, int num_files)
-{
- DynamicTimer dt(__func__);
-
- num_files = fix_file_counts(num_files);
- const idx_t num_reads = get_num_reads(input);
- const idx_t num_batches = (num_reads + batch_size - 1) / batch_size;
- string idx_file_name;
+void partition_candidates(const std::string& input, const std::string& pac_prefix, const size_t file_size, const int max_files_per_batch, const int64_t num_reads) {
+ DynamicTimer dtimer(__func__);
+ struct stat buf;
+ if (stat(input.c_str(), &buf) == -1) {
+ ERROR("Could not get file size: %s", input.c_str());
+ }
+ // each candidate line takes on average 44 characters, but go with 32;
+ // each one produces two candidates (forward and reverse)
+ const int num_files(ceil(double(buf.st_size / 32 * 2) * sizeof(ExtensionCandidateCompressed) / file_size));
+ // separate them by read id as we don't know how many candidates each
+ // read is part of; we're assuming an even distribution on average
+ const int64_t reads_per_file((num_reads + num_files - 1) / num_files);
+ std::vector read_sizes;
+ PackedDB::read_sizes(pac_prefix, read_sizes);
+ PartitionResultsWriter prw(max_files_per_batch);
+ int i(0);
+ off_t input_pos;
+ // if restarting, both i and input_pos can be changed
+ int is_restart(prw.restart(input, generate_partition_file_name, "partition.done", i, input_pos));
+ std::string idx_file_name;
generate_partition_index_file_name(input, idx_file_name);
- ofstream idx_file;
- open_fstream(idx_file, idx_file_name.c_str(), ios::out);
-
- ExtensionCandidate ec, nec;
- PartitionResultsWriter prw(num_files);
- for (idx_t i = 0; i < num_batches; i += num_files) {
- const idx_t sfid = i;
- const idx_t efid = min(sfid + num_files, num_batches);
- const int nf = efid - sfid;
- const idx_t Lid = batch_size * sfid;
- const idx_t Rid = batch_size * efid;
- cout << "Lid = " << Lid
- << ", Rid = " << Rid
- << "\n";
- ifstream in;
- open_fstream(in, input, ios::in);
- prw.OpenFiles(sfid, efid, input, generate_partition_file_name);
-
+ std::ofstream idx_file;
+ open_fstream(idx_file, idx_file_name.c_str(), std::ios::out);
+ ExtensionCandidate ec;
+ ExtensionCandidateCompressed nec;
+ // and here we go through the input file to write num_files,
+ // being limited by the number of open output files we can have
+ for (; i < num_files; i += prw.kNumFiles) {
+ const int sfid(i);
+ const int efid(std::min(sfid + prw.kNumFiles, num_files));
+ const int nf(efid - sfid);
+ const int64_t L(sfid * reads_per_file);
+ const int64_t R(std::min(efid * reads_per_file, num_reads));
+ std::ifstream in;
+ open_fstream(in, input.c_str(), std::ios::in);
+ if (is_restart) {
+ if (!in.seekg(input_pos)) {
+ ERROR("Input seek failed while restoring checkpoint: %s", input.c_str());
+ }
+ is_restart = 0;
+ } else {
+ prw.OpenFiles(sfid, efid, input, generate_partition_file_name, "partition.done");
+ }
while (in >> ec) {
- if (ec.qsize < min_read_size || ec.ssize < min_read_size) continue;
- if (ec.qid >= Lid && ec.qid < Rid) {
- normalise_candidate(ec, nec, false);
- prw.WriteOneResult((ec.qid - Lid) / batch_size, ec.qid, nec);
+ // size has been set to zero for too-small reads
+ if (!read_sizes[ec.sid] || !read_sizes[ec.qid]) {
+ continue;
}
- if (ec.sid >= Lid && ec.sid < Rid)
- {
- normalise_candidate(ec, nec, true);
- prw.WriteOneResult((ec.sid - Lid) / batch_size, ec.sid, nec);
+ // make sure values match before tossing the values
+ assert(ec.ssize == read_sizes[ec.sid] && ec.qsize == read_sizes[ec.qid]);
+ // as variable sizes may be different, make sure the ones
+ // we read in can be safely stored
+ assert(ec.sid <= ExtensionCandidateCompressed::max_value);
+ assert(ec.qid <= ExtensionCandidateCompressed::max_value);
+ assert(ec.sext <= ExtensionCandidateCompressed::max_value);
+ assert(ec.score <= ExtensionCandidateCompressed::max_value);
+ // qext is one bit smaller than the others
+ assert(ec.qext <= ExtensionCandidateCompressed::max_qext);
+ if (L <= ec.sid && ec.sid < R) {
+ nec.set(ec);
+ if (prw.WriteOneResult((nec.sid - L) / reads_per_file, nec.sid, nec)) {
+ prw.checkpoint(in.tellg());
+ }
+ }
+ if (L <= ec.qid && ec.qid < R) {
+ nec.set_swap(ec);
+ if (prw.WriteOneResult((nec.sid - L) / reads_per_file, nec.sid, nec)) {
+ prw.checkpoint(in.tellg());
+ }
}
}
- for (int k = 0; k < nf; ++k)
- {
- if (prw.max_seq_ids[k] == std::numeric_limits::min()) continue;
- idx_file << prw.file_names[k] << "\t" << prw.min_seq_ids[k] << "\t" << prw.max_seq_ids[k] << "\n";
- fprintf(stderr, "%s contains reads %d --- %d\n", prw.file_names[k].c_str(), (int)prw.min_seq_ids[k], (int)prw.max_seq_ids[k]);
+ for (int k(0); k < nf; ++k) {
+ std::cerr << prw.file_names[k] << " contains " << prw.counts[k] << " overlaps\n";
+ if (prw.counts[k] != 0) {
+ idx_file << prw.file_names[k] << "\n";
+ }
}
prw.CloseFiles();
}
close_fstream(idx_file);
+ prw.finalize();
}
-/*
-void
-partition_candidates(const char* input, const idx_t batch_size, const int min_read_size, const int num_files)
-{
+void partition_m4records(const char* const m4_file_name, const double min_cov_ratio, const size_t file_size, const int min_read_size, const int max_files_per_batch) {
DynamicTimer dtimer(__func__);
-
- idx_t num_reads = get_num_reads(input);
- const index_t num_batches = (num_reads + batch_size - 1) / batch_size;
+ int64_t num_reads(get_qualified_m4record_counts(m4_file_name, min_cov_ratio));
+ std::set repeat_reads;
+ //get_repeat_reads(m4_file_name, min_cov_ratio, num_reads, repeat_reads);
+ struct stat buf;
+ if (stat(m4_file_name, &buf) == -1) {
+ ERROR("Could not get file size: %s", m4_file_name);
+ }
+ const int num_files(ceil(double(buf.st_size / 32 * 2) * sizeof(ExtensionCandidate) / file_size));
+ const int64_t reads_per_file((num_reads + num_files - 1) / num_files);
std::string idx_file_name;
- generate_partition_index_file_name(input, idx_file_name);
- std::ofstream idx_file;
- open_fstream(idx_file, idx_file_name.c_str(), std::ios::out);
-
- ExtensionCandidate ec, nec;
- PartitionResultsWriter prw;
- for (index_t i = 0; i < num_batches; i += PartitionResultsWriter::kNumFiles)
- {
- const index_t sfid = i;
- const index_t efid = std::min(sfid + PartitionResultsWriter::kNumFiles, num_batches);
- const int nf = efid - sfid;
- const index_t L = batch_size * sfid;
- const index_t R = batch_size * efid;
- std::ifstream in;
- open_fstream(in, input, std::ios::in);
- prw.OpenFiles(sfid, efid, input, generate_partition_file_name);
-
- while (in >> ec)
- {
- if (ec.qsize < min_read_size || ec.ssize < min_read_size) continue;
- if (ec.qid >= L && ec.qid < R)
- {
- normalise_candidate(ec, nec, false);
- prw.WriteOneResult((ec.qid - L) / batch_size, ec.qid, nec);
+ generate_partition_index_file_name(m4_file_name, idx_file_name);
+ std::ofstream idx_file;
+ open_fstream(idx_file, idx_file_name.c_str(), std::ios::out);
+ M4Record m4, nm4;
+ ExtensionCandidate ec;
+ PartitionResultsWriter prw(max_files_per_batch);
+ for (int i(0); i < num_files; i += prw.kNumFiles) {
+ const int sfid(i);
+ const int efid(std::min(sfid + prw.kNumFiles, num_files));
+ const int nf(efid - sfid);
+ const int64_t L(reads_per_file * sfid);
+ const int64_t R(efid < num_files ? reads_per_file * efid : num_reads);
+ std::ifstream in;
+ open_fstream(in, m4_file_name, std::ios::in);
+ prw.OpenFiles(sfid, efid, m4_file_name, generate_partition_file_name, "partition.done");
+ while (in >> m4) {
+ if (m4qsize(m4) < min_read_size || m4ssize(m4) < min_read_size) {
+ continue;
+ } else if (!check_m4record_mapping_range(m4, min_cov_ratio)) {
+ continue;
+ } else if (repeat_reads.find(m4qid(m4)) != repeat_reads.end() || repeat_reads.find(m4sid(m4)) != repeat_reads.end()) {
+ continue;
}
- if (ec.sid >= L && ec.sid < R)
- {
- normalise_candidate(ec, nec, true);
- prw.WriteOneResult((ec.sid - L) / batch_size, ec.sid, nec);
+ if (m4qid(m4) >= L && m4qid(m4) < R) {
+ normalize_m4record(m4, false, nm4);
+ m4_to_candidate(nm4, ec);
+ prw.WriteOneResult((m4qid(m4) - L) / reads_per_file, m4qid(m4), ec);
+ }
+ if (m4sid(m4) >= L && m4sid(m4) < R) {
+ normalize_m4record(m4, true, nm4);
+ m4_to_candidate(nm4, ec);
+ prw.WriteOneResult((m4sid(m4) - L) / reads_per_file, m4sid(m4), ec);
+ }
+ }
+ for (int k(0); k < nf; ++k) {
+ std::cerr << prw.file_names[k] << " contains " << prw.counts[k] << " overlaps\n";
+ if (prw.counts[k] != 0) {
+ idx_file << prw.file_names[k] << "\n";
}
}
-
- for (int k = 0; k < nf; ++k)
- {
- if (prw.max_seq_ids[k] == std::numeric_limits::min()) continue;
- idx_file << prw.file_names[k] << "\t" << prw.min_seq_ids[k] << "\t" << prw.max_seq_ids[k] << "\n";
- fprintf(stderr, "%s contains reads %d --- %d\n", prw.file_names[k].c_str(), (int)prw.min_seq_ids[k], (int)prw.max_seq_ids[k]);
- }
+ prw.CloseFiles();
}
- prw.CloseFiles();
close_fstream(idx_file);
}
-*/
-
-/*
-void
-partition_m4records(const char* m4_file_name, const double min_cov_ratio, const index_t batch_size, const int min_read_size)
-{
- DynamicTimer dtimer(__func__);
-
- index_t num_reads, num_qualified_records;
- get_qualified_m4record_counts(m4_file_name, min_cov_ratio, num_qualified_records, num_reads);
- set repeat_reads;
- //get_repeat_reads(m4_file_name, min_cov_ratio, num_reads, repeat_reads);
- const index_t num_batches = (num_reads + batch_size - 1) / batch_size;
- std::string idx_file_name;
- generate_partition_index_file_name(m4_file_name, idx_file_name);
- std::ofstream idx_file;
- open_fstream(idx_file, idx_file_name.c_str(), std::ios::out);
-
- M4Record m4, nm4;
- ExtensionCandidate ec;
- PartitionResultsWriter prw;
- for (index_t i = 0; i < num_batches; i += PartitionResultsWriter::kNumFiles)
- {
- const index_t sfid = i;
- const index_t efid = std::min(sfid + PartitionResultsWriter::kNumFiles, num_batches);
- const int nf = efid - sfid;
- const index_t L = batch_size * sfid;
- const index_t R = batch_size * efid;
- std::ifstream in;
- open_fstream(in, m4_file_name, std::ios::in);
- prw.OpenFiles(sfid, efid, m4_file_name, generate_partition_file_name);
-
- while (in >> m4)
- {
- if (m4qsize(m4) < min_read_size || m4ssize(m4) < min_read_size) continue;
- if (!check_m4record_mapping_range(m4, min_cov_ratio)) continue;
- if (repeat_reads.find(m4qid(m4)) != repeat_reads.end()
- ||
- repeat_reads.find(m4sid(m4)) != repeat_reads.end()) continue;
-
- if (m4qid(m4) >= L && m4qid(m4) < R)
- {
- normalize_m4record(m4, false, nm4);
- m4_to_candidate(nm4, ec);
- prw.WriteOneResult((m4qid(m4) - L) / batch_size, m4qid(m4), ec);
- }
- if (m4sid(m4) >= L && m4sid(m4) < R)
- {
- normalize_m4record(m4, true, nm4);
- m4_to_candidate(nm4, ec);
- prw.WriteOneResult((m4sid(m4) - L) / batch_size, m4sid(m4), ec);
- }
- }
- for (int k = 0; k < nf; ++k)
- {
- if (prw.max_seq_ids[k] == std::numeric_limits::min()) continue;
- idx_file << prw.file_names[k] << "\t" << prw.min_seq_ids[k] << "\t" << prw.max_seq_ids[k] << "\n";
- fprintf(stderr, "%s contains reads %d --- %d\n", prw.file_names[k].c_str(), (int)prw.min_seq_ids[k], (int)prw.max_seq_ids[k]);
- }
-
- prw.CloseFiles();
- }
- close_fstream(idx_file);
-}
-*/
-
-void
-partition_m4records(const char* m4_file_name,
- const double min_cov_ratio,
- const index_t batch_size,
- const int min_read_size,
- int num_files)
-{
- DynamicTimer dtimer(__func__);
-
- num_files = fix_file_counts(num_files);
- index_t num_reads, num_qualified_records;
- get_qualified_m4record_counts(m4_file_name, min_cov_ratio, num_qualified_records, num_reads);
- set repeat_reads;
- const index_t num_batches = (num_reads + batch_size - 1) / batch_size;
- std::string idx_file_name;
- generate_partition_index_file_name(m4_file_name, idx_file_name);
- std::ofstream idx_file;
- open_fstream(idx_file, idx_file_name.c_str(), std::ios::out);
-
- M4Record m4, nm4;
- ExtensionCandidate ec;
- PartitionResultsWriter prw(num_files);
- for (index_t i = 0; i < num_batches; i += num_files)
- {
- const index_t sfid = i;
- const index_t efid = std::min(sfid + num_files, num_batches);
- const int nf = efid - sfid;
- const index_t L = batch_size * sfid;
- const index_t R = batch_size * efid;
- std::ifstream in;
- open_fstream(in, m4_file_name, std::ios::in);
- prw.OpenFiles(sfid, efid, m4_file_name, generate_partition_file_name);
-
- while (in >> m4)
- {
- if (m4qsize(m4) < min_read_size || m4ssize(m4) < min_read_size) continue;
- if (!check_m4record_mapping_range(m4, min_cov_ratio)) continue;
- if (repeat_reads.find(m4qid(m4)) != repeat_reads.end()
- ||
- repeat_reads.find(m4sid(m4)) != repeat_reads.end()) continue;
-
- if (m4qid(m4) >= L && m4qid(m4) < R)
- {
- normalize_m4record(m4, false, nm4);
- m4_to_candidate(nm4, ec);
- prw.WriteOneResult((m4qid(m4) - L) / batch_size, m4qid(m4), ec);
- }
- if (m4sid(m4) >= L && m4sid(m4) < R)
- {
- normalize_m4record(m4, true, nm4);
- m4_to_candidate(nm4, ec);
- prw.WriteOneResult((m4sid(m4) - L) / batch_size, m4sid(m4), ec);
- }
- }
-
- for (int k = 0; k < nf; ++k)
- {
- if (prw.max_seq_ids[k] == std::numeric_limits::min()) continue;
- idx_file << prw.file_names[k] << "\t" << prw.min_seq_ids[k] << "\t" << prw.max_seq_ids[k] << "\n";
- fprintf(stderr, "%s contains reads %d --- %d\n", prw.file_names[k].c_str(), (int)prw.min_seq_ids[k], (int)prw.max_seq_ids[k]);
- }
-
- prw.CloseFiles();
- }
- close_fstream(idx_file);
-}
-
-void
-load_partition_files_info(const char* idx_file_name, std::vector& file_info_vec)
-{
- file_info_vec.clear();
- std::ifstream in;
- open_fstream(in, idx_file_name, std::ios::in);
- PartitionFileInfo pfi;
- while (in >> pfi.file_name)
- {
- in >> pfi.min_seq_id >> pfi.max_seq_id;
- file_info_vec.push_back(pfi);
- }
- close_fstream(in);
+void load_partition_files_info(const char* const idx_file_name, std::vector& file_info_vec) {
+ file_info_vec.clear();
+ std::ifstream in;
+ open_fstream(in, idx_file_name, std::ios::in);
+ std::string line;
+ while (in >> line) {
+ file_info_vec.push_back(line);
+ }
+ close_fstream(in);
}
diff --git a/src/mecat2cns/overlaps_partition.h b/src/mecat2cns/overlaps_partition.h
index 275a972..f8f7cc3 100644
--- a/src/mecat2cns/overlaps_partition.h
+++ b/src/mecat2cns/overlaps_partition.h
@@ -1,37 +1,16 @@
#ifndef OVERLAPS_PARTITION_H
#define OVERLAPS_PARTITION_H
-#include
+#include // string
+#include // vector<>
+#include // int64_t
-#include "../common/alignment.h"
+void generate_partition_index_file_name(const std::string& input_file_name, std::string& ret);
-void
-generate_partition_index_file_name(const char* m4_file_name, std::string& ret);
+void partition_m4records(const char* m4_file_name, double min_cov_ratio, size_t batch_size, int min_read_size, int num_files);
-void
-generate_partition_file_name(const char* m4_file_name, const index_t part, std::string& ret);
+void partition_candidates(const std::string& input, const std::string& pac_prefix, size_t batch_size, int num_files, int64_t num_reads);
-void
-partition_m4records(const char* m4_file_name,
- const double min_cov_ratio,
- const index_t batch_size,
- const int min_read_size,
- int num_files);
-
-void
-partition_candidates(const char* input,
- const idx_t batch_size,
- const int min_read_size,
- int num_files);
-
-struct PartitionFileInfo
-{
- std::string file_name;
- index_t min_seq_id;
- index_t max_seq_id;
-};
-
-void
-load_partition_files_info(const char* idx_file_name, std::vector& file_info_vec);
+void load_partition_files_info(const char* idx_file_name, std::vector& file_info_vec);
#endif // OVERLAPS_PARTITION_H
diff --git a/src/mecat2cns/overlaps_store.h b/src/mecat2cns/overlaps_store.h
index a545982..e015ae1 100644
--- a/src/mecat2cns/overlaps_store.h
+++ b/src/mecat2cns/overlaps_store.h
@@ -1,121 +1,219 @@
#ifndef _OVERLAPS_STORE_H
#define _OVERLAPS_STORE_H
-#include
-#include
-#include
+#include // ifstream, ofstream, streampos, streamsize
+#include // string
+#include // int64_t
+#include // time(), time_t
+#include // _SC_OPEN_MAX, F_OK, access(), off_t, rename(), sysconf(), unlink()
+#include // vector<>
-#include "../common/defs.h"
-#include "../common/pod_darr.h"
+#include "../common/defs.h" // ERROR(), LOG()
+#include "../common/pod_darr.h" // PODArray<>
-template
-class PartitionResultsWriter
-{
-public:
- typedef void (*file_name_generator)(const char* prefix, const idx_t id, std::string& name);
-
-public:
- PartitionResultsWriter(const int max_num_files)
- {
- file_is_open = 0;
- num_open_files = 0;
- MaxNumFiles = max_num_files;
- results = new PODArray[MaxNumFiles];
- files = new std::ofstream[MaxNumFiles];
- file_names = new std::string[MaxNumFiles];
- min_seq_ids = new idx_t[MaxNumFiles];
- max_seq_ids = new idx_t[MaxNumFiles];
- for (int i = 0; i < MaxNumFiles; ++i) results[i].reserve(kStoreSize);
+template class PartitionResultsWriter {
+ public:
+ typedef void (*file_name_generator)(const std::string& prefix, int id, std::string& name);
+ public:
+ const int kNumFiles; // effective open file limit
+ int kStoreSize;
+ int num_open_files;
+ PODArray* results; // can't use vector<>, causes memory corruption
+ std::ofstream* files; // can't use vector<>, non-copyable
+ std::vector file_names;
+ std::vector counts;
+ public:
+ // can't make kNumFiles static, as sysconf() is run-time only;
+ // leave room for stdin, stdout, stderr, a few others
+ explicit PartitionResultsWriter(const int max_files_per_batch) : kNumFiles(max_files_per_batch > 0 ? max_files_per_batch : sysconf(_SC_OPEN_MAX) - 10), kStoreSize(0), num_open_files(0), results(0), files(0) { }
+ ~PartitionResultsWriter() {
+ CloseFiles();
+ }
+ void OpenFiles(const int sfid, const int efid, const std::string& prefix, file_name_generator fng, const std::string& done_file) {
+ CloseFiles();
+ if (efid <= sfid) {
+ return;
+ }
+ done_file_ = done_file;
+ ckpt_file_ = done_file_ + ".ckpt";
+ ckpt_file_tmp_ = ckpt_file_ + ".tmp";
+ num_open_files = efid - sfid;
+ batch_start_ = sfid;
+ allocate_data(prefix, fng, 0);
}
-
- ~PartitionResultsWriter()
- {
+ void CloseFiles() {
+ if (num_open_files == 0) {
+ return;
+ }
+ for (int i(0); i != num_open_files; ++i) { // flush buffers
+ if (results[i].size()) {
+ write_buffer_to_disk(i);
+ }
+ files[i].close();
+ // don't attempt to re-finish finished files
+ if (!file_names[i].empty()) {
+ const std::string tmp_file(file_names[i] + ".tmp");
+ if (rename(tmp_file.c_str(), file_names[i].c_str()) == -1) {
+ ERROR("Could not rename %s", tmp_file.c_str());
+ }
+ }
+ }
delete[] results;
delete[] files;
- delete[] file_names;
- delete[] min_seq_ids;
- delete[] max_seq_ids;
+ file_names.clear();
+ counts.clear();
+ kStoreSize = num_open_files = 0;
+ results = 0;
+ files = 0;
}
-
- void OpenFiles(const idx_t sfid, const idx_t efid, const std::string& prefix, file_name_generator fng)
- {
- if (file_is_open) CloseFiles();
-
- const int nf = efid - sfid;
- if (nf == 0) return;
- for (int i = 0; i < nf; ++i)
- {
- fng(prefix.data(), i + sfid, file_names[i]);
- open_fstream(files[i], file_names[i].c_str(), std::ios::binary);
- min_seq_ids[i] = std::numeric_limits::max();
- max_seq_ids[i] = std::numeric_limits::min();
+ void finalize() {
+ // touch done file
+ std::ofstream done(done_file_.c_str());
+ if (!done) {
+ ERROR("Could not create done file %s", done_file_.c_str());
+ }
+ done.close();
+ // remove checkpoint file
+ unlink(std::string(done_file_ + ".ckpt").c_str());
+ }
+ int WriteOneResult(const int i, const int64_t seq_id, const T& r) {
+ ++counts[i];
+ results[i].push_back(r);
+ if (results[i].size() == kStoreSize) {
+ write_buffer_to_disk(i);
results[i].clear();
- }
- num_open_files = nf;
- file_is_open = true;
+ if (time(0) >= next_checkpoint_time_) {
+ return 1;
+ }
+ }
+ return 0;
+ }
+ int restart(const std::string& prefix, file_name_generator fng, const std::string& done_file, int& sfid, off_t& input_pos) {
+ CloseFiles();
+ done_file_ = done_file;
+ ckpt_file_ = done_file_ + ".ckpt";
+ ckpt_file_tmp_ = ckpt_file_ + ".tmp";
+ std::ifstream in(ckpt_file_.c_str());
+ if (!in) {
+ return 0;
+ }
+ in >> batch_start_ >> num_open_files >> input_pos;
+ if (!in) {
+ ERROR("Read error while restoring checkpoint from %s", ckpt_file_.c_str());
+ }
+ allocate_data(prefix, fng, 1);
+ for (int i(0); i != num_open_files; ++i) {
+ off_t file_pos;
+ in >> file_pos >> counts[i];
+ if (!in) {
+ ERROR("Read error while restore checkpoint from %s (%d)", ckpt_file_.c_str(), i);
+ }
+ // don't need to restart finished files
+ if (!file_names[i].empty() && !files[i].seekp(file_pos)) {
+ ERROR("Seek failed while restoring checkpoint from %s (%d)", ckpt_file_.c_str(), i);
+ }
+ }
+ sfid = batch_start_;
+ return 1;
}
-
- void CloseFiles()
- {
- if (!file_is_open) return;
- for (int i = 0; i < num_open_files; ++i)
- {
- if (results[i].size())
- {
- char* buf = (char*)results[i].data();
- std::streamsize s = sizeof(T) * results[i].size();
- files[i].write(buf, s);
+ void checkpoint(const off_t input_pos) {
+ std::ofstream out(ckpt_file_tmp_.c_str());
+ if (!out) {
+ LOG(stderr, "Checkpoint failed: couldn't open %s", ckpt_file_tmp_.c_str());
+ return;
+ }
+ out << batch_start_ << " " << num_open_files << " " << input_pos << "\n";
+ if (!out) {
+ LOG(stderr, "Checkpoint failed: write failed: %s", ckpt_file_tmp_.c_str());
+ return;
+ }
+ // flush buffers
+ for (int i(0); i != num_open_files; ++i) {
+ if (results[i].size()) {
+ write_buffer_to_disk(i);
+ results[i].clear();
+ }
+ files[i].flush();
+ out << off_t(files[i].tellp()) << " " << counts[i] << "\n";
+ if (!out) {
+ LOG(stderr, "Checkpoint failed: write failed: %s", ckpt_file_tmp_.c_str());
+ return;
}
}
- for (int i = 0; i < num_open_files; ++i) close_fstream(files[i]);
- file_is_open = false;
- num_open_files = 0;
+ out.close();
+ if (rename(ckpt_file_tmp_.c_str(), ckpt_file_.c_str()) == -1) {
+ LOG(stderr, "Checkpoint failed: rename failed: %s", ckpt_file_.c_str());
+ }
+ next_checkpoint_time_ = time(0) + 300;
}
-
- void WriteOneResult(const int fid, const idx_t seq_id, const T& r)
- {
- if (fid >= num_open_files) {
- std::cout << "fid = " << fid
- << ", num_open_files = "
- << num_open_files
- << "\n";
+ int64_t total_count() const {
+ int64_t total(0);
+ std::vector::const_iterator a(counts.begin());
+ const std::vector::const_iterator end_a(counts.end());
+ for (; a != end_a; ++a) {
+ total += *a;
}
- r_assert(fid < num_open_files);
- min_seq_ids[fid] = std::min(min_seq_ids[fid], seq_id);
- max_seq_ids[fid] = std::max(max_seq_ids[fid], seq_id);
- results[fid].push_back(r);
- if (results[fid].size() == kStoreSize)
- {
- char* buf = (char*)results[fid].data();
- std::streamsize s = sizeof(T) * results[fid].size();
- files[fid].write(buf, s);
- results[fid].clear();
+ return total;
+ }
+ private:
+ int batch_start_ ;
+ time_t next_checkpoint_time_;
+ std::string done_file_;
+ std::string ckpt_file_;
+ std::string ckpt_file_tmp_;
+ private:
+ void allocate_data(const std::string& prefix, file_name_generator fng, const int is_restart) {
+ // allocate about a gb of memory as buffer, split among num_open_files
+ kStoreSize = (1 << 30) / sizeof(T) / num_open_files;
+ results = new PODArray[num_open_files];
+ file_names.assign(num_open_files, "");
+ counts.assign(num_open_files, 0);
+ files = new std::ofstream[num_open_files];
+ for (int i(0); i != num_open_files; ++i) {
+ fng(prefix, i + batch_start_, file_names[i]);
+ const std::string tmp_file(file_names[i] + ".tmp");
+ if (is_restart) {
+ if (access(file_names[i].c_str(), F_OK) == 0) {
+ // mark as already finished
+ file_names[i].clear();
+ // use /dev/null to prevent write errors
+ files[i].open("/dev/null", std::ios::binary);
+ } else {
+ // don't truncate on restart
+ files[i].open(tmp_file.c_str(), std::ios::binary | std::ios::in);
+ }
+ } else {
+ files[i].open(tmp_file.c_str(), std::ios::binary);
+ }
+ if (!files[i]) {
+ ERROR("Open failed on %s", tmp_file.c_str());
+ }
+ results[i].reserve(kStoreSize);
+ }
+ next_checkpoint_time_ = time(0) + 300;
+ }
+ void write_buffer_to_disk(const int i) {
+ // can't use static_cast<>
+ const char* const buf((char*)results[i].data());
+ const std::streamsize s(sizeof(T) * results[i].size());
+ if (!files[i].write(buf, s)) {
+ ERROR("Error writing to %s", file_names[i].c_str());
}
}
-
-public:
- int MaxNumFiles;
- static const int kStoreSize = 500000;
- PODArray* results;
- bool file_is_open;
- int num_open_files;
- std::ofstream* files;
- std::string* file_names;
- idx_t* min_seq_ids;
- idx_t* max_seq_ids;
};
-template
-T* load_partition_data(const char* path, idx_t& num_results)
-{
+template T* load_partition_data(const char* const path, int64_t& num_results) {
std::ifstream in;
open_fstream(in, path, std::ios::binary);
in.seekg(0, std::ios::end);
- std::streampos fs = in.tellg();
+ const std::streampos fs(in.tellg());
in.seekg(0, std::ios::beg);
num_results = fs / sizeof(T);
- T* arr = new T[num_results];
- in.read((char*)arr, fs);
+ T* const arr(new T[num_results]);
+ in.read((char*)arr, fs); // can't use static_cast<>
+ if (!in) {
+ ERROR("Error reading partition data: %s", path);
+ }
close_fstream(in);
return arr;
}
diff --git a/src/mecat2cns/packed_db.cpp b/src/mecat2cns/packed_db.cpp
new file mode 100644
index 0000000..b4682db
--- /dev/null
+++ b/src/mecat2cns/packed_db.cpp
@@ -0,0 +1,318 @@
+#include "packed_db.h"
+
+#include // ifstream, ofstream
+#include // set<>
+#include // uint8_t
+#include // rename()
+#include
+#include