diff --git a/graphUtils.cpp b/graphUtils.cpp index 19909a9..02e7700 100644 --- a/graphUtils.cpp +++ b/graphUtils.cpp @@ -13,8 +13,7 @@ void graphUtils::read_graph() n_vtx = gfa_n_vtx(g); // Resize node_len node_len.resize(n_vtx, 0); - // Array of vectors - adj_ = new std::vector[n_vtx]; + adj_.resize(n_vtx); /* Node_len */ for (int v = 0; v < n_vtx/2; v++) { @@ -811,14 +810,11 @@ std::vector graphUtils::Chaining(std::vector anchors) int N = M[cid].size(); // #Anchors int K = path_cover[cid].size(); // #Paths /* Initialise Search Trees */ - std::pair , int64_t> default_value = {{std::numeric_limits::min(), -1}, (int64_t)-1}; - std::vector I(K, SearchTree(default_value)); // Search Tree /* Initialise T */ std::vector T; // Tuples of Anchors - std::vector> C; // Array of Cost of each Anchor + std::vector> C; // Array of Cost of each Anchor C.resize(N); - int64_t sf = scale_factor; // Scale Factor - int64_t cost =(int64_t) (M[cid][0].d - M[cid][0].c + 1)*sf; // Cost of each Anchor + int64_t cost =(int64_t) (M[cid][0].d - M[cid][0].c + 1); // Cost of each Anchor for (int j = 0; j < N; j++) // Add Tuple of Anchors { int node = M[cid][j].v; @@ -885,13 +881,25 @@ std::vector graphUtils::Chaining(std::vector anchors) int64_t _infty_int64 = std::numeric_limits::min(); std::pair , int64_t> rmq; // rmq std::pair key; + + float c_1 = expf(-div*(float)kmer_len), c_2 = 0.05*c_1; // penalities + if(!is_ggen) // minimap2 gap cost + { + c_1 = 0.01f * kmer_len; // minimap2 + c_2 = 0.0f; + } + + std::vector, std::pair>> MAP(K); + std::map, std::pair>::iterator start_it, end_it; + for (auto t:T) // Process Tuples in the Lexicographic order of (rank(v), pos, task) { if(t.task == 0) // update the score { - int64_t val_1 = ( M[cid][t.anchor].c - 1 + M[cid][t.anchor].x - 1 + dist2begin[cid][t.path][t.v] + Distance[cid][t.path][t.w]); - int64_t val_2 = sf*(M[cid][t.anchor].d - M[cid][t.anchor].c + 1); - int64_t range = (int64_t)(dist2begin[cid][t.path][t.v] + Distance[cid][t.path][t.w] + M[cid][t.anchor].x - G - 1); + int64_t weight = (M[cid][t.anchor].d - M[cid][t.anchor].c + 1); + int64_t range = (dist2begin[cid][t.path][t.v] + Distance[cid][t.path][t.w] + M[cid][t.anchor].x - G - 1); + int64_t Rd_1 = (dist2begin[cid][t.path][t.v] + Distance[cid][t.path][t.w] + M[cid][t.anchor].x); + int64_t Qd_1 = (M[cid][t.anchor].c); if (index[cid][t.path][t.v] != -1) // same path { for (int l = x[t.path]; l < path_distance[t.path].size(); l++) @@ -907,35 +915,43 @@ std::vector graphUtils::Chaining(std::vector anchors) for (int l = x[t.path]; l < rmq_coor[t.path]; l++) { key = path_distance[t.path][l].second; - I[t.path].remove(key); + MAP[t.path].erase(key); } // Update the pointer x[t.path] = rmq_coor[t.path]; } - // Extend the chain by adding the current anchor - rmq = I[t.path].RMQ_2({M[cid][t.anchor].c_, infty_int}, {M[cid][t.anchor].c - 1, infty_int}, range); - if (rmq.first.first > _infty_int64) - { - C[t.anchor] = std::max(C[t.anchor], { rmq.first.first - val_1 + val_2, rmq.first.second}); - } - if (param_z) // debug + // Minigraph gap cost + start_it = MAP[t.path].lower_bound({M[cid][t.anchor].c_, infty_int}); + end_it = MAP[t.path].upper_bound({M[cid][t.anchor].d - 1, infty_int}); // overlap + // // Maximize the score function from window query + std::pair max_cost = std::make_pair(_infty_int64, -1); + int h = 0; + for (auto i = std::reverse_iterator(end_it); i != std::reverse_iterator(start_it) && h < max_itr; i++) { - std::cerr << " cid : " << cid << " idx : " << t.anchor << " top_v :" << t.top_v << " pos : " << t.pos << " task : " << t.task << " path : " << t.path << " parent : " << C[t.anchor].second << " node : " << M[cid][t.anchor].v << " index : " << index[cid][t.path][t.w] << " C[j] : " << C[t.anchor].first << " update_C : " << (rmq.first.first - val_1 + val_2) << " rmq.first : " << rmq.first.first << " val_1 : " << val_1 << " dist2begin : " << dist2begin[cid][t.path][t.v] << " Distance : " << Distance[cid][t.path][t.w] << " M[i].d : " << t.d << "\n"; + int idx = i->second.first; // anchor_idx + int v = i->second.second; // vertex + int64_t Rd_2 = dist2begin[cid][t.path][v] + M[cid][idx].y + 1; // count of bases hence +1 + int64_t Qd_2 = M[cid][idx].d + 1; + if (Rd_2 - 1 > range) + { + int64_t g = abs((Rd_1 - Rd_2) - (Qd_1 - Qd_2)); + int64_t d = std::min((Rd_1 - Rd_2), (Qd_1 - Qd_2)); + float log_pen = g >= 2 ? mg_log2(g) : 0.0f; // mg_log2() only works for dd>=2 + int64_t gap = (c_1*(float)g + c_2*(float)d + log_pen); + weight = std::min(d, kmer_len); + max_cost = max_cost.first > (C[idx].first - gap) ? max_cost : std::make_pair(C[idx].first - gap, idx); + C[t.anchor] = std::max(C[t.anchor], { weight + max_cost.first, max_cost.second}); + h++; + } } }else // update the priority of the anchor { - int64_t val_3 = (M[cid][t.anchor].d + M[cid][t.anchor].y + dist2begin[cid][t.path][t.v]); // M[i].d + M[i].y + dist2begin[v] - I[t.path].add({M[cid][t.anchor].d, t.anchor}, std::make_pair(std::make_pair(C[t.anchor].first + val_3 , t.anchor), (int64_t)(M[cid][t.anchor].y + dist2begin[cid][t.path][t.v]))); // C[j].first + M[i].d + M[i].y + dist2begin[v] (priority, anchor) + MAP[t.path][{M[cid][t.anchor].d, t.anchor}] = {t.anchor, t.v}; path_distance[t.path].push_back({(int64_t)(M[cid][t.anchor].y + dist2begin[cid][t.path][t.v]), {M[cid][t.anchor].d, t.anchor}}); // (value, key) pair - if (param_z) // debug - { - std::cerr << " cid : " << cid << " idx : " << t.anchor << " top_v :" << t.top_v << " pos : " << t.pos << " task : " << t.task << " path : " << t.path << " val_3 : " << val_3 << " M.y : " << M[cid][t.anchor].y << " dist2begin : " << dist2begin[cid][t.path][t.v] << " M[i].d : " << t.d << "\n"; - } } } - // Backtracking for read mapping and graph generation if (N!=0) { @@ -955,12 +971,12 @@ std::vector graphUtils::Chaining(std::vector anchors) return std::tie(a.first.first, a.first.second) > std::tie(b.first.first, b.first.second); // sort by score and index }); // Find the minimum score - int min_score = scale_factor * (float)(M[cid][0].d - M[cid][0].c + 1); // minimum score for a contig + int min_score = (M[cid][0].d - M[cid][0].c + 1); // minimum score for a contig // Find the maximum score std::pair chain_rmq = D[0].first; // Max value of (score, index) int prev_idx = 0; // index of the maximum score max_score = chain_rmq.first; // maximum score - int64_t threshold_score = is_ggen > 0 ? tau_1*(float)min_score : tau_1*(float)max_score; // threshold score + int64_t threshold_score = is_ggen > 0 ? 1000 : tau_1*(float)max_score; // threshold score int64_t best_cid_score = max_score; // best score for a contig std::pair, int64_t> chain_pair; // (chain, score) bool flag; // flag for disjoint chain @@ -988,7 +1004,7 @@ std::vector graphUtils::Chaining(std::vector anchors) } } // Store the anchors if chain is disjoint and clear the keys - if (flag == true && temp_chain.size() > 3 ) // minimap2 min_cnt + if (flag == true && temp_chain.size()) // minimap2 min_cnt { for (int i = 0; i < temp_chain.size(); i++) { @@ -1068,24 +1084,6 @@ std::vector graphUtils::Chaining(std::vector anchors) } // Reverse the order of the elements in best std::reverse(best.begin(),best.end()); - std::vector red_idx; - std::vector count(node_len.size(), 0); // Initialize the count array with all 0s - /* Count the frequency of each node */ - for (int i = 0; i < best.size(); i++) { - int node = (int)(best[i].x >> 32); - count[node]++; - } - /* Find the indices of anchors with frequency <= 5 */ - for (int i = 0; i < best.size(); i++) { - int node = (int)(best[i].x >> 32); - if (count[node] <= 5) { - red_idx.push_back(i); - } - } - /* Remove anchors from collected indices */ - for (int i = red_idx.size() - 1; i >= 0; i--) { - best.erase(best.begin()+red_idx[i]); - } } if (param_z) { diff --git a/graphUtils.h b/graphUtils.h index 066862c..7e12f7b 100644 --- a/graphUtils.h +++ b/graphUtils.h @@ -15,226 +15,7 @@ #include "minigraph.h" #include "mgpriv.h" #include -#include - -// AVL Tree -template -class AVLTree { - public: - AVLTree(const V &default_value = V()) : default_value(default_value) , root_(-1) {} - - void add(T key, V value) { root_ = add(root_, key, value); } - - void update(T key, V value) { - auto node = find(root_, key); - if (node != -1) { - nodes_[node].value = value; - } - } - - V RMQ(T key1, T key2) { - return RMQ(root_, key1, key2); - } - - V RMQ_1(T key1, T key2) { - return RMQ_1(root_, key1, key2); - } - - V RMQ_2(T key1, T key2, int64_t range) { - return RMQ_2(root_, key1, key2, range); - } - - void remove(T key) { root_ = remove(root_, key); } - - V get(T key) { - int node = find(root_, key); - if (node != -1) { - return nodes_[node].value; - } - return default_value; - } - - private: - struct Node { - T key; - V value; - int height; - int left; - int right; - - Node(T k, V v) : key(k), value(v), height(1), left(-1), right(-1) {} - }; - int root_; - std::vector nodes_; - V default_value; - - int height(int node) { return node != -1 ? nodes_[node].height : 0; } - - int balanceFactor(int node) { return height(nodes_[node].left) - height(nodes_[node].right); } - - void updateHeight(int node) { nodes_[node].height = std::max(height(nodes_[node].left), height(nodes_[node].right)) + 1; } - - int rotateRight(int node) { - int left = nodes_[node].left; - nodes_[node].left = nodes_[left].right; - nodes_[left].right = node; - updateHeight(node); - updateHeight(left); - return left; - } - - int rotateLeft(int node) { - int right = nodes_[node].right; - nodes_[node].right = nodes_[right].left; - nodes_[right].left = node; - updateHeight(node); - updateHeight(right); - return right; - } - - int balance(int node) { - updateHeight(node); - if (balanceFactor(node) == 2) { - if (balanceFactor(nodes_[node].left) < 0) { - nodes_[node].left = rotateLeft(nodes_[node].left); - } - return rotateRight(node); - } else if (balanceFactor(node) == -2) { - if (balanceFactor(nodes_[node].right) > 0) { - nodes_[node].right = rotateRight(nodes_[node].right); - } - return rotateLeft(node); - } - return node; - } - - int add(int node, T key, V value) { - if (node == -1) { - nodes_.emplace_back(key, value); - return (int)nodes_.size() - 1; - } - if (key < nodes_[node].key) { - nodes_[node].left = add(nodes_[node].left, key, value); - } else if (key > nodes_[node].key) { - nodes_[node].right = add(nodes_[node].right, key, value); - } else { - nodes_[node].value = std::max(nodes_[node].value, value); - } - return balance(node); - } - - int find(int node, T key) { - if (node == -1) { - return -1; - } - if (key < nodes_[node].key) { - return find(nodes_[node].left, key); - } else if (key > nodes_[node].key) { - return find(nodes_[node].right, key); - } else { - return node; - } - } - - int findMin(int node) { - if (nodes_[node].left == -1) { - return node; - } - return findMin(nodes_[node].left); - } - - - int remove(int node, T key) { - if (node == -1) { - return -1; - } - if (key < nodes_[node].key) { - nodes_[node].left = remove(nodes_[node].left, key); - } else if (key > nodes_[node].key) { - nodes_[node].right = remove(nodes_[node].right, key); - } else { - if (nodes_[node].left == -1 && nodes_[node].right == -1) { - return -1; - } - if (nodes_[node].left == -1) { - return nodes_[node].right; - } - if (nodes_[node].right == -1) { - return nodes_[node].left; - } - int next = findMin(nodes_[node].right); - nodes_[node].key = nodes_[next].key; - nodes_[node].value = nodes_[next].value; - nodes_[node].right = remove(nodes_[node].right, nodes_[next].key); - } - return balance(node); - } - - V RMQ(int node, T key1, T key2) { - if (node == -1) { - return default_value; - } - if (key1 <= nodes_[node].key && key2 >= nodes_[node].key) { - V leftMax = RMQ(nodes_[node].left, key1, key2); - V rightMax = RMQ(nodes_[node].right, key1, key2); - return std::max(nodes_[node].value, std::max(leftMax, rightMax)); - } else if (key1 > nodes_[node].key) { - return RMQ(nodes_[node].right, key1, key2); - } else { - return RMQ(nodes_[node].left, key1, key2); - } - } - - V RMQ_1(int node, T key1, T key2) { - if (node == -1) { - return default_value; - } - V maxValue = default_value; - std::stack s; - s.push(node); - while (!s.empty()) { - node = s.top(); - s.pop(); - if (key1 <= nodes_[node].key && key2 >= nodes_[node].key) { - maxValue = std::max(maxValue, nodes_[node].value); - } - if (nodes_[node].left != -1 && key1 <= nodes_[node].key) { - s.push(nodes_[node].left); - } - if (nodes_[node].right != -1 && key2 > nodes_[node].key) { - s.push(nodes_[node].right); - } - } - return maxValue; - } - - V RMQ_2(int node, T key1, T key2, int64_t range) { - if (node == -1) { - return default_value; - } - V maxValue = default_value; - std::stack s; - s.push(node); - while (!s.empty()) { - node = s.top(); - std::pair, int64_t> value = nodes_[node].value; - s.pop(); - if (key1 <= nodes_[node].key && key2 >= nodes_[node].key && value.second > range) { - maxValue = std::max(maxValue, nodes_[node].value); - } - if (nodes_[node].left != -1 && key1 <= nodes_[node].key) { - s.push(nodes_[node].left); - } - if (nodes_[node].right != -1 && key2 > nodes_[node].key) { - s.push(nodes_[node].right); - } - } - return maxValue; - } - -}; - -typedef AVLTree, std::pair, int64_t>> SearchTree; +#include // Anchors @@ -293,7 +74,7 @@ class graphUtils { public: gfa_t * g; // This is Graph - std::vector *adj_; // This is the adjacency list + std::vector> adj_; // This is the adjacency list std::vector> conn_comp; // connected components std::vector component; // component id int num_comp; // number of connected components @@ -327,6 +108,9 @@ class graphUtils float tau_1; float tau_2; bool is_ggen; + int64_t kmer_len; + float div; + int max_itr; graphUtils(gfa_t *g); // This is constructor diff --git a/index.c b/index.c index bcb5e8f..59cffe5 100644 --- a/index.c +++ b/index.c @@ -211,11 +211,10 @@ mg_idx_t *mg_index_core(gfa_t *g, int k, int w, int b, int n_threads) /* Pass parameters */ params* par; -void pass_par(bool ¶m_z, float &scale_factor, int &G) +void pass_par(bool ¶m_z, int &G) { par = new params(); par->param_z = param_z; - par->scale_factor = scale_factor; par->G = G; } @@ -240,7 +239,6 @@ mg_idx_t *mg_index(gfa_t *g, const mg_idxopt_t *io, int n_threads, mg_mapopt_t * graphOp->read_graph(); omp_set_dynamic(1); omp_set_num_threads(0); - graphOp->scale_factor = par->scale_factor; // graphOp->print_graph(); graphOp->Connected_components(); int cycle_count = graphOp->is_cyclic(); diff --git a/lchain.c b/lchain.c index e135932..fd858f0 100644 --- a/lchain.c +++ b/lchain.c @@ -126,7 +126,7 @@ static inline int32_t comput_sc(const mg128_t *ai, const mg128_t *aj, int32_t ma dg = dr < dq? dr : dq; q_span = aj->y>>32&0xff; sc = q_span < dg? q_span : dg; - if (dd || dg > q_span) { + if (dd || dg > q_span && false) { float lin_pen, log_pen; lin_pen = chn_pen_gap * (float)dd + chn_pen_skip * (float)dg; log_pen = dd >= 1? mg_log2(dd + 1) : 0.0f; // mg_log2() only works for dd>=2 @@ -241,7 +241,7 @@ static inline int32_t comput_sc_simple(const mg128_t *ai, const mg128_t *aj, flo q_span = aj->y>>32&0xff; sc = q_span < dg? q_span : dg; if (exact) *exact = (dd == 0 && dg <= q_span); - if (dd || dq > q_span) { + if (dd || dq > q_span && false) { float lin_pen, log_pen; lin_pen = chn_pen_gap * (float)dd + chn_pen_skip * (float)dg; log_pen = dd >= 1? mg_log2(dd + 1) : 0.0f; // mg_log2() only works for dd>=2 diff --git a/main.c b/main.c index 4a352f0..47dd928 100644 --- a/main.c +++ b/main.c @@ -7,6 +7,7 @@ #include "ketopt.h" #include "graphUtils.h" #include +#include #ifdef __linux__ #include @@ -102,9 +103,8 @@ int main(int argc, char *argv[]) mg_mapopt_t opt; mg_idxopt_t ipt; mg_ggopt_t gpt; - int i, c, ret, n_threads = 4; - float scale_factor = 200.0f; - int G = 10000; // read mapping + int i, c, ret, n_threads = sysconf(_SC_NPROCESSORS_ONLN); // get maximal threads as default + int G = 5000; // read mapping bool z = false; char *s; FILE *fp_help = stderr; @@ -135,7 +135,6 @@ int main(int argc, char *argv[]) if (c == 'w') ipt.w = atoi(o.arg); else if (c == 'k') ipt.k = atoi(o.arg); else if (c == 't') n_threads = atoi(o.arg); - else if (c == 's') scale_factor = atof(o.arg); else if (c == 'z') z = atoi(o.arg); else if (c == 'G') G = atoi(o.arg); else if (c == 'f') opt.occ_max1_frac = atof(o.arg); @@ -253,7 +252,6 @@ int main(int argc, char *argv[]) fprintf(fp_help, " -N INT retain at most INT secondary mappings [%d]\n", opt.best_n); fprintf(fp_help, " -D skip self diagonal matches\n"); fprintf(fp_help, " -z BOOL print chain information [%d]\n", z); - fprintf(fp_help, " -s FLOAT factor to scale anchor weights during chaining [%f]\n", scale_factor); fprintf(fp_help, " Graph generation:\n"); fprintf(fp_help, " --ggen perform incremental graph generation\n"); fprintf(fp_help, " -q INT min mapping quality [%d]\n", gpt.min_mapq); @@ -276,7 +274,7 @@ int main(int argc, char *argv[]) return fp_help == stdout? 0 : 1; } - pass_par(z, scale_factor, G); + pass_par(z, G); g = gfa_read(argv[o.ind]); if (g == 0) { diff --git a/map-algo.c b/map-algo.c index 8aeb9fa..0c2dc51 100644 --- a/map-algo.c +++ b/map-algo.c @@ -404,7 +404,10 @@ void mg_map_frag(const mg_idx_t *gi, int n_segs, const int *qlens, const char ** graphOp->tau_1 = opt->tau_1; graphOp->tau_2 = opt->tau_2; graphOp->is_ggen = opt->is_ggen; - if (graphOp->G == 10000) // No user input + graphOp->kmer_len = gi->k; + graphOp->div = opt->div; + graphOp->max_itr = opt->max_itr; + if (graphOp->G == 5000) // No user input { graphOp->G = opt->G; } // else G will be passed from user input diff --git a/mgpriv.h b/mgpriv.h index f25d577..efbc8dc 100644 --- a/mgpriv.h +++ b/mgpriv.h @@ -121,10 +121,9 @@ void radix_sort_128x(mg128_t *beg, mg128_t *end); void radix_sort_gfa64(uint64_t *beg, uint64_t *end); uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); -void pass_par(bool ¶m_z, float &scale_factor, int &G); +void pass_par(bool ¶m_z, int &G); struct params { - float scale_factor; bool param_z; int G; }; diff --git a/minigraph.h b/minigraph.h index 42919c2..dfed8d2 100644 --- a/minigraph.h +++ b/minigraph.h @@ -5,7 +5,7 @@ #include "gfa.h" #define MG_VERSION "0.20-r559" -#define MC_VERSION "1.1" +#define MC_VERSION "1.2" #define MG_M_SPLICE 0x10 #define MG_M_SR 0x20 @@ -79,6 +79,7 @@ typedef struct { int min_cov_mapq, min_cov_blen; bool is_ggen; int G; + int max_itr; } mg_mapopt_t; typedef struct { @@ -178,4 +179,4 @@ int mg_ggen(gfa_t *g, int32_t n_fn, const char **fn, const mg_idxopt_t *ipt, con } #endif -#endif \ No newline at end of file +#endif diff --git a/options.c b/options.c index c11a6d0..0617bc3 100644 --- a/options.c +++ b/options.c @@ -28,8 +28,8 @@ void mg_mapopt_init(mg_mapopt_t *mo) mo->mini_batch_size = 500000000; mo->div = 0.1f; mo->l_chn_pen_gap = 0.0f, mo->l_chn_pen_skip = 0.0f; // Added for linear chaining - mo->chn_pen_gap = 1.0f, mo->chn_pen_skip = 0.05f; - mo->min_lc_cnt = 2, mo->min_lc_score = 0; // mo->min_lc_cnt = 5, mo->min_lc_score = 40; + mo->chn_pen_gap = 0.0f, mo->chn_pen_skip = 0.00f; + mo->min_lc_cnt = 6, mo->min_lc_score = 0; mo->min_gc_cnt = 5, mo->min_gc_score = 50; mo->gdp_max_ed = 10000; mo->lc_max_trim = 50; @@ -44,11 +44,13 @@ void mg_mapopt_init(mg_mapopt_t *mo) mo->min_cov_blen = 1000; mo->cap_kalloc = 1000000000; // tau_1 : threshold to compute "intra cid" disjoint set of chains - mo->tau_1 = 0.99; + mo->tau_1 = 0.99f; // tau_2 : threshold to pick "inter cid" set of chains - mo->tau_2 = 0.95; + mo->tau_2 = 0.95f; mo->is_ggen = 0; - mo->G = 10000; // for long read mapping + mo->G = 5000; // for long read mapping + mo->max_itr = 100; + } void mg_ggopt_init(mg_ggopt_t *go) @@ -77,24 +79,25 @@ int mg_opt_set(const char *preset, mg_idxopt_t *io, mg_mapopt_t *mo, mg_ggopt_t mg_ggopt_init(go); } else if (strcmp(preset, "lr") == 0) { // this is the default } else if (strcmp(preset, "asm") == 0 || strcmp(preset, "ggs") == 0) { - io->k = 17, io->w = 11; + io->k = 19, io->w = 10; mo->flag |= MG_M_RMQ; - mo->occ_max1 = 10, mo->occ_max1_cap = 100; + mo->occ_max1 = 10, mo->occ_max1_cap = 100; mo->bw = 1000, mo->bw_long = 150000; - mo->max_gap = 10000, mo->max_gap_pre = 1000; - mo->min_lc_cnt = 1, mo->min_lc_score = 0; // minigraph default: mo->min_lc_cnt = 5, mo->min_lc_score = 40; - mo->min_gc_cnt = 5, mo->min_gc_score = 1000; + mo->max_gap = 150000, mo->max_gap_pre = 1000; + mo->min_lc_cnt = 1, mo->min_lc_score = 0; + mo->min_gc_cnt = 1, mo->min_gc_score = 1000; mo->min_cov_mapq = 5; mo->min_cov_blen = 100000; mo->max_lc_skip = mo->max_gc_skip = 50; mo->div = 0.01f; mo->mini_batch_size = 4000000000LL; // tau_1 : threshold to compute "intra cid" disjoint set of chains - mo->tau_1 = 5.00f; // to match minimap2 defualt threshold + mo->tau_1 = 0.0f; // tau_2 : threshold to pick "inter cid" set of chains - mo->tau_2 = 0.95f; + mo->tau_2 = 0.0f; mo->is_ggen = 1; - mo->G = 5000; // for graph generation + mo->G = 150000; // for graph generation + mo->max_itr = 500; // works reasonably well if (strcmp(preset, "ggs") == 0) go->algo = MG_G_GGSIMPLE, mo->best_n = 0; } else if (strcmp(preset, "se") == 0 || strcmp(preset, "sr") == 0) {