xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-01 03:08:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

ont simplex support

Browse Source
This commit is contained in:
chhylp123
2024-11-27 13:27:13 -05:00
parent 80fa5ed436
commit 676385cf8e
10 changed files with 896 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+38 -11
View File
@@ -7,6 +7,9 @@
#include <sys/time.h>
#include "CommandLines.h"
#include "ketopt.h"
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)
#define DEFAULT_OUTPUT "hifiasm.asm"
@@ -72,7 +75,9 @@ static ko_longopt_t long_options[] = {
{ "telo-s", ko_required_argument, 357},
{ "ctg-n", ko_required_argument, 358},
{ "ont", ko_no_argument, 359},
{ "sc-n", ko_no_argument, 360},
// { "sc-n", ko_no_argument, 360},
{ "chem-c", ko_required_argument, 361},
{ "chem-f", ko_required_argument, 362},
// { "path-round", ko_required_argument, 348},
{ 0, 0, 0 }
};
@@ -210,6 +215,14 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " --telo-s INT\n");
fprintf(stderr, " min score for telomere reads [%ld]\n", asm_opt->telo_mic_sc);
fprintf(stderr, " ONT simplex assembly (beta):\n");
fprintf(stderr, " --ont assemble ONT simplex reads in fastq format\n");
// fprintf(stderr, " --sc-n consider base qual value for assembly\n");
fprintf(stderr, " --chem-c INT\n");
fprintf(stderr, " detect chemical reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov);
fprintf(stderr, " --chem-f INT\n");
fprintf(stderr, " length of flanking regions for chemical read detection [%lu]\n", asm_opt->chemical_flank);
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n");
@@ -342,6 +355,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->is_ont = 0;
asm_opt->is_sc = 0;
asm_opt->chemical_cov = 1;
asm_opt->chemical_flank = 256;
}
void destory_enzyme(enzyme* f)
@@ -692,25 +707,33 @@ int check_option(hifiasm_opt_t* asm_opt)
void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
{
if(opt->ind == argc)
{
if(opt->ind == argc) {
return;
}
asm_opt->num_reads = argc - opt->ind;
asm_opt->read_file_names = (char**)malloc(sizeof(char*)*asm_opt->num_reads);
long long i;
gzFile dfp;
for (i = 0; i < asm_opt->num_reads; i++)
{
long long i; int ret;
gzFile dfp; kseq_t *ks = NULL;
for (i = 0; i < asm_opt->num_reads; i++) {
asm_opt->read_file_names[i] = argv[i + opt->ind];
dfp = gzopen(asm_opt->read_file_names[i], "r");
if (dfp == 0)
{
if (dfp == 0) {
fprintf(stderr, "[ERROR] Cannot find the input read file: %s\n",
asm_opt->read_file_names[i]);
exit(0);
} else if(asm_opt->is_sc){
ks = kseq_init(dfp);
while (((ret = kseq_read(ks)) >= 0)) {
if((ks->qual.l == 0) || (ks->qual.s == NULL)) {
fprintf(stderr, "[ERROR] %s is in fasta format rather than fastq format\n", asm_opt->read_file_names[i]);
asm_opt->is_sc = 0;
exit(0);
}
break;
}
kseq_destroy(ks); ks = NULL;
}
gzclose(dfp);
}
@@ -914,9 +937,13 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 357) asm_opt->telo_mic_sc = atol(opt.arg);
else if (c == 358) asm_opt->max_contig_tip = atol(opt.arg);
else if (c == 359) {
asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37;
} else if (c == 360) {
asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; asm_opt->is_sc = 1; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37;
} /**else if (c == 360) {
asm_opt->is_sc = 1;
}**/ else if (c == 361) {
asm_opt->chemical_cov = atol(opt.arg);
} else if (c == 362) {
asm_opt->chemical_flank = atol(opt.arg);
} else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
}
CommandLines.h
+3 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.21.0-r666"
#define HA_VERSION "0.21.0-r686"
#define VERBOSE 0
@@ -163,6 +163,8 @@ typedef struct {
uint64_t is_ont;
uint64_t is_sc;
uint64_t chemical_cov;
uint64_t chemical_flank;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Correct.cpp
+108 -2
View File
@@ -24531,7 +24531,108 @@ void reassign_gaps(overlap_region *z, overlap_region *aux, char* qstr, int64_t q
// }
}
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf)
uint32_t is_ovlp_debug(int64_t s, int64_t e, int64_t ws, int64_t we, int64_t op)
{
int64_t os, oe, ovlp;
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0);
if(op != 2) {
return (!!ovlp);
} else {
if(ws >= s && we <= e) return 1;
}
return 0;
}
uint32_t inline ff_tend(overlap_region *z, int64_t wn, int64_t dn, double dr, double er, int64_t min_err)
{
int64_t k, zwn = z->w_list.n, err, mm, qi, ci, cn, ql, ws, we, zs, ze, s, e, os, oe, ovlp; bit_extz_t ez; uint32_t cl; uint16_t c;
zs = z->x_pos_s; ze = z->x_pos_e + 1; ql = ze - zs;
if(ql < wn) return 0; if(dn > (ql*dr)) dn = ql*dr; if(dn < wn) return 0; if(ql < dn) return 0;
s = zs; e = zs + dn;
// if(z->y_id == 27) fprintf(stderr, "-a-[M::%s] tid::%u\t%.*s\tz::[%ld,%ld)\ti::[%ld,%ld)\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), zs, ze, s, e);
for (k = err = mm = 0, qi = zs; (k < zwn) && (z->w_list.a[k].x_start < e); k++) {
// if(z->y_id == 27) fprintf(stderr, "-a-[M::%s] tid::%u\t%.*s\tw::[%d,%d)\terr::%d\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].error);
if(!(is_ualn_win(z->w_list.a[k]))) {
set_bit_extz_t(ez, (*z), k);
ci = 0; cn = ez.cigar.n; qi = ez.ts; //ti = ez.ps;
while (ci < cn && qi < e) {
ws = qi;
ci = pop_trace(&(ez.cigar), ci, &c, &cl);
if(c!=2) qi += cl;
// if(c!=3) ti += cl;
we = qi;
if(c == 0) {
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0); mm += ovlp;
} else {
err += cl;
}
// assert(is_ovlp_debug(s, e, ws, we, c));
if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) {
// fprintf(stderr, "-0-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er);
return 1;
}
}
} else {
ws = z->w_list.a[k].x_start; we = z->w_list.a[k].x_end + 1;
err += we - ws;
// assert(is_ovlp_debug(s, e, ws, we, -1));
if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) {
// fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er);
return 1;
}
}
}
s = ze - dn; e = ze;
for (k = zwn - 1, err = mm = 0, qi = ze; (k >= 0) && ((z->w_list.a[k].x_end + 1) > s); k--) {
if(!(is_ualn_win(z->w_list.a[k]))) {
set_bit_extz_t(ez, (*z), k);
ci = ((int64_t)ez.cigar.n) - 1; qi = ez.te + 1; //ti = ez.pe + 1;
while (ci >= 0 && qi > s) {
we = qi;
ci = pop_trace_back(&(ez.cigar), ci, &c, &cl);
if(c!=2) qi -= cl;
// if(c!=3) ti += cl;
ws = qi;
if(c == 0) {
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0); mm += ovlp;
} else {
err += cl;
}
// assert(is_ovlp_debug(s, e, ws, we, c));
if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) {
// fprintf(stderr, "-2-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er);
return 1;
}
}
} else {
ws = z->w_list.a[k].x_start; we = z->w_list.a[k].x_end + 1;
err += we - ws;
// assert(is_ovlp_debug(s, e, ws, we, -1));
if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) {
// fprintf(stderr, "-3-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er);
return 1;
}
}
}
return 0;
}
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop)
{
uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re;
overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1;
@@ -24565,6 +24666,9 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre
if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue;
if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue;
// if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re);
reassign_gaps(z, aux_o, qu->seq, ql, NULL, -1, rref, tu, buf);
@@ -24586,7 +24690,7 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre
}
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf)
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop)
{
uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re;
overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1;
@@ -24628,6 +24732,8 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads
if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue;
if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue;
// if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re);
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
Correct.h
+2 -2
View File
@@ -1391,8 +1391,8 @@ int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *t
const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o,
bit_extz_t *exz, double e_rate, int64_t qs);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop);
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop);
uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8);
void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
Hash_Table.h
+3 -2
View File
@@ -9,8 +9,9 @@
#define WINDOW 375
#define WINDOW_BOUNDARY 375
#define WINDOW_HC 775
#define WINDOW_OHC 475
// #define WINDOW_OHC 375 ///ONT high error
///ONT high error
// #define WINDOW_OHC 475
#define WINDOW_OHC 375
#define WINDOW_HC_FAST 512
///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
#define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
Overlaps.cpp
+58 -9
View File
@@ -19,6 +19,7 @@
#include "gfa_ut.h"
#include "assert.h"
#include "khash.h"
#include "ecovlp.h"
KHASH_SET_INIT_INT64(64)
@@ -63,7 +64,6 @@ KRADIX_SORT_INIT(ha_mzl_t_srt1, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x)
#define Uc_beg(z) ((uint32_t)((z).a[0]>>32))
#define Uc_end(z) ((uint32_t)((z).a[(z).n-1]>>32)^1)
#define UL_COV_THRES 2
#define PHASE_SEP 64
#define PHASE_SEF 2
#define PHASE_SEP_RATE 0.04
@@ -2969,10 +2969,18 @@ int max_hang, int min_ovlp)
void prt_specific_overlap(ma_hit_t_alloc *src, uint64_t qn, uint64_t tn, const char *cmd)
{
int64_t idx = get_specific_overlap(&(src[qn]), qn, tn);
const ma_hit_t *h = &(src[qn].buffer[idx]);
fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__,
Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml);
ma_hit_t *h = NULL; int64_t idx = -1, k;
if(tn != ((uint64_t)-1)) {
idx = get_specific_overlap(&(src[qn]), qn, tn); h = &(src[qn].buffer[idx]);
fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__,
Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml);
} else {
for (k = 0; k < src[qn].length; k++) {
h = &(src[qn].buffer[k]);
fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__,
Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml);
}
}
}
asg_t *ma_sg_gen_ul(ma_hit_t_alloc* sources, int64_t n_read, const ma_sub_t *coverage_cut,
@@ -5667,8 +5675,12 @@ long long weakID, uint32_t w_qs, uint32_t w_qe)
{
strongID = Get_tn(aim_paf->buffer[i]);
index = get_specific_overlap(&(reverse_paf_list[strongID]), strongID, weakID);
if(index != -1)
{
if(index != -1) {
// if((Get_qn(aim_paf->buffer[i]) == 27087 && weakID == 27128) || (Get_qn(aim_paf->buffer[i]) == 27128 && weakID == 27087)) {
// ma_hit_t *h = &(aim_paf->buffer[i]);
// fprintf(stderr, "[M::%s]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", __func__, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev],
// (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl);
// }
return 0;
}
}
@@ -11193,6 +11205,11 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
index = get_specific_overlap(&(sources[tn]), tn, qn);
// if(index < 0 || index >= sources[tn].length) fprintf(stderr, "sb, tn: %u, qn: %u, index: %ld, length: %u\n", tn, qn, index, sources[tn].length);
sources[tn].buffer[index].bl |= ((uint32_t)0x40000000);
// ma_hit_t *h = &(sources[i].buffer[j]);
// if((Get_qn(*h) == 27087 && Get_tn(*h) == 27128) || (Get_tn(*h) == 27087 && Get_qn(*h) == 27128)) {
// fprintf(stderr, "[M::%s]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", __func__, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev],
// (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl);
// }
}
}
}
@@ -38528,17 +38545,43 @@ void prt_dbg_gfa(asg_t *sg, const char *suffix, ma_sub_t *cov, ma_hit_t_alloc* s
free(gfa_name); free(o_file);
}
void prt_dbg_rid_ovlp(ma_hit_t_alloc *ov, int64_t rid, char *rn, const char *cmd)
{
uint64_t k; ma_hit_t *h = NULL;
if(rid < 0) {
for (k = 0; k < R_INF.total_reads; k++) {
if (memcmp(rn, Get_NAME((R_INF), k), Get_NAME_LENGTH((R_INF), k)) == 0) break;
}
if(k >= R_INF.total_reads) return;
rid = k;
}
fprintf(stderr, "\n[M::%s::%s::id::%ld]\n", __func__, cmd, rid);
for (k = 0; k < ov[rid].length; k++) {
h = &(ov[rid].buffer[k]);
if(h->del) continue;
fprintf(stderr, "%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev],
(int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl);
}
}
asg_t *gen_init_sg(int32_t min_dp, uint64_t n_read, int64_t mini_overlap_length, int64_t max_hang_length, int64_t gap_fuzz,
ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t, ma_sub_t** cov, all_ul_t *ul, telo_end_t *te)
{
asg_t *sg = NULL;
// prt_specific_overlap(src, 22233, 22235, "1");
// prt_dbg_rid_ovlp(src, 3700, NULL, "1");
if(ul) rescue_src_ul(src, n_read, UL_COV_THRES);
// prt_dbg_rid_ovlp(src, 3700, NULL, "2");
ma_hit_sub(min_dp, src, n_read, readLen, mini_overlap_length, cov);
// prt_dbg_rid_ovlp(src, 3700, NULL, "3");
detect_chimeric_reads(src, n_read, readLen, *cov, asm_opt.max_ov_diff_final*2.0, ul, UL_COV_THRES);
// prt_dbg_rid_ovlp(src, 3700, NULL, "4");
ma_hit_cut(src, n_read, readLen, mini_overlap_length, cov);
// prt_dbg_rid_ovlp(src, 3700, NULL, "5");
ma_hit_flt(src, n_read, *cov, max_hang_length, mini_overlap_length);
// prt_dbg_rid_ovlp(src, 3700, NULL, "6");
ma_hit_contained_advance(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length);
// prt_dbg_rid_ovlp(src, 3700, NULL, "7");
if(!ul) {
sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length);
@@ -38650,6 +38693,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
}
///just for debug
renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read);
// if(asm_opt.is_ont) handle_chemical_arc(asm_opt.thread_num, R_INF.total_reads);
// if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads);
///it's hard to say which function is better
///normalize_ma_hit_t_single_side(sources, n_read);
@@ -38681,7 +38726,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
}
// prt_specific_overlap(sources, 22233, 22235, "0-b");
// prt_specific_overlap(sources, 22235, 22233, "0-b");
clean_weak_ma_hit_t(sources, reverse_sources, n_read, asm_opt.ar?UL_COV_THRES:(uint32_t)-1);
// prt_dbg_rid_ovlp(sources, 27087, NULL, "0-a");
if(!(asm_opt.is_ont)) clean_weak_ma_hit_t(sources, reverse_sources, n_read, asm_opt.ar?UL_COV_THRES:(uint32_t)-1);
// prt_dbg_rid_ovlp(sources, 27087, NULL, "0-b");
// prt_specific_overlap(sources, 22233, 22235, "0-c");
// prt_specific_overlap(sources, 22235, 22233, "0-c");
sg = gen_init_sg(min_dp, n_read, mini_overlap_length, max_hang_length, gap_fuzz, sources, readLen, ruIndex,
@@ -39006,7 +39053,9 @@ long long bubble_dist, int read_graph, int write)
if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY))
{
if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads);
try_rescue_overlaps(sources, reverse_sources, n_read, 4);
clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
Overlaps.h
+1
View File
@@ -35,6 +35,7 @@
// #define ALTER_LABLE 2
// #define HAP_LABLE 4
#define ug_ext_len 75000
#define UL_COV_THRES 2
#define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
#define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
ecovlp.cpp
+501 -12
View File
@@ -14,6 +14,56 @@
#define del_cns_nn(z, nn_i) ((z).a[(nn_i)].sc == CNS_DEL_V)
#define REFRESH_N 128
KDQ_INIT(uint32_t)
typedef struct {
uint32_t v:31, f:1;
uint32_t sc;
} cns_arc;
typedef struct {size_t n, m, nou; cns_arc *a; } cns_arc_v;
typedef struct {
// uint16_t c:2, t:2, f:1, sc:3;
uint32_t c:2, f:1, sc:29;
cns_arc_v arc;
}cns_t;
typedef struct {
size_t n, m;
cns_t *a;
uint32_t si, ei, off, bn, bb0, bb1, cns_g_wl;
kdq_t(uint32_t) *q;
}cns_gfa;
typedef struct {
// chaining and overlapping related buffers
UC_Read self_read, ovlp_read;
Candidates_list clist;
overlap_region_alloc olist;
ha_abuf_t *ab;
// int64_t num_read_base, num_correct_base, num_recorrect_base;
uint64_t cnt[6], rr;
haplotype_evdience_alloc hap;
bit_extz_t exz;
kv_ul_ov_t pidx;
asg64_v v64;
asg32_v v32;
asg16_v v16;
asg8_v v8q, v8t;
kvec_t_u8_warp k_flag;
st_mt_t sp;
cns_gfa cns;
} ec_ovec_buf_t0;
typedef struct {
ec_ovec_buf_t0 *a;
uint32_t n, rev;
} ec_ovec_buf_t;
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n);
void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p);
#define generic_key(x) (x)
KRADIX_SORT_INIT(ec16, uint16_t, generic_key, 2)
@@ -2608,6 +2658,7 @@ void push_ff_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag,
if(z->el == 1) cnt[4]++;
if(z->no_l_indel) cnt[5]++;
}
z->del = 0;
}
}
@@ -2729,7 +2780,7 @@ inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl,
return 0;
}
void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in)
void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop)
{
if(ol->length <= 0) return;
@@ -2744,7 +2795,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *
}
if(!(srt->n)) {
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop);
} else {
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
@@ -2778,7 +2829,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
if(on > nec) gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf);
if(on > nec) gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop);
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
@@ -2986,6 +3037,158 @@ void debug_retrive_bqual(asg8_v *vq, asg8_v *vt, uint64_t id, uint64_t rn)
}
}
uint32_t is_uncorrected_read(overlap_region_alloc* ov, asg64_v *idx, int64_t len, int64_t min_len)
{
uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed;
for (k = idx->n = 0; k < ov->length; k++) {
s = ov->list[k].x_pos_s; e = ov->list[k].x_pos_e + 1;
kv_push(uint64_t, (*idx), (s<<1));
kv_push(uint64_t, (*idx), (e<<1)|1);
}
radix_sort_ec64(idx->a, idx->a + idx->n);
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
if(old_dp == 0) {
if((ed - st) >= min_len) return 1;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
if(old_dp == 0) {
if((ed - st) >= min_len) return 1;
if((ed - st) >= len) return 1;
}
}
return 0;
}
uint32_t is_chemical_r_qual(overlap_region_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t frank_len, asg8_v *qv, uint64_t rid)
{
uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, s1, e0, e1, rr, qk;
if((frank_len) > (len *0.01)) frank_len = len *0.01;
for (k = idx->n = 0; k < ov->length; k++) {
s = ov->list[k].x_pos_s; e = ov->list[k].x_pos_e + 1;
kv_push(uint64_t, (*idx), (s<<1));
kv_push(uint64_t, (*idx), (e<<1)|1);
// fprintf(stderr, "[M::%s]\trid::%lu\ts::%lu\te::%lu\n", __func__, rid, s, e);
}
radix_sort_ec64(idx->a, idx->a + idx->n); s0 = s1 = e0 = e1 = rr = -1;
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
if(old_dp <= cov) {
rr = 1;
} else {
if(s0 < 0) {
s0 = st; s1 = ed;
}
e0 = st; e1 = ed;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
if(old_dp <= cov) {
rr = 1;
} else {
if(s0 < 0) {
s0 = st; s1 = ed;
}
e0 = st; e1 = ed;
}
}
if((s0 != e0) && (s1 != e1) && (s0 <= frank_len) && ((len - e1) <= frank_len) && (rr > 0)) {
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
if((old_dp <= cov) && (st >= s0) && (ed <= e1)) {
retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn);
fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]);
fprintf(stderr, "\n");
return 1;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
if((old_dp <= cov) && (st >= s0) && (ed <= e1)) {
retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn);
fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]);
fprintf(stderr, "\n");
return 1;
}
}
}
// for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
// old_dp = dp;
// ///if a[j] is qe
// if (idx->a[k]&1) --dp;
// else ++dp;
// ed = idx->a[k]>>1;
// if(ed > st) {
// if((old_dp <= cov)) {
// retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn);
// fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
// for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]);
// fprintf(stderr, "\n");
// return 1;
// }
// }
// st = ed;
// }
// ed = len; old_dp = dp;
// if(ed > st) {
// if((old_dp <= cov)) {
// retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn);
// fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
// for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]);
// fprintf(stderr, "\n");
// return 1;
// }
// }
return 0;
}
static void worker_hap_ec(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
@@ -3006,7 +3209,7 @@ static void worker_hap_ec(void *data, long i, int tid)
// if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
// if (memcmp("c42804f3-0e13-43a0-8a71-b91b40accf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// if (memcmp("b2e68ecf-381a-439c-b676-c1e6831d6acf", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// if (memcmp("0aec8c4f-c849-4c31-85ba-4ffb297eeb28", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// if (memcmp("64b2c27d-86b8-451e-9330-6ba62be2ffcc", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
@@ -3023,7 +3226,7 @@ static void worker_hap_ec(void *data, long i, int tid)
recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length;
h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**((asm_opt.is_ont)?(0.05):(0.02))**/0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error
h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error
// b->num_read_base += b->olist.length;
b->cnt[0] += b->self_read.length;
@@ -3033,7 +3236,7 @@ static void worker_hap_ec(void *data, long i, int tid)
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]));
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont);
// prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length);
@@ -3061,7 +3264,10 @@ static void worker_hap_ec(void *data, long i, int tid)
push_nec_re(aux_o, &(scc.a[i]));
push_nec_re(aux_o, &(scb.a[i]));
// if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) {
// // b->olist.length = 0;
// fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
// }
push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])/**, i, &b->self_read, &b->ovlp_read**/);
push_ne_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, NULL/**, i, NULL, NULL**/);
@@ -3605,9 +3811,9 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid)
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; uint32_t qlen = 0;
recover_UC_Read(&b->self_read, &R_INF, i);
recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length;
h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**0.02**/0.001, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);
@@ -3624,6 +3830,11 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid)
h_ec_lchain_fast_new(b->ab, i, &b->self_read, &b->ovlp_read, &R_INF, &b->olist, &b->clist, &b->exz, &b->v16, &b->v64, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), 0.866666);
if((asm_opt.is_ont) && (is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))) {
b->olist.length = 0;
// fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
push_ff_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, b->cnt);
push_ff_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, b->cnt);
@@ -3641,6 +3852,236 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid)
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
uint32_t is_chemical_r(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t frank_len)
{
uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, s1, e0, e1, rr;
if((frank_len) > (len *0.01)) frank_len = len *0.01;
for (k = idx->n = 0; k < ov->length; k++) {
s = (uint32_t)ov->buffer[k].qns; e = ov->buffer[k].qe;
kv_push(uint64_t, (*idx), (s<<1));
kv_push(uint64_t, (*idx), (e<<1)|1);
}
radix_sort_ec64(idx->a, idx->a + idx->n); s0 = s1 = e0 = e1 = rr = -1;
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
if(old_dp <= cov) {
rr = 1;
} else {
if(s0 < 0) {
s0 = st; s1 = ed;
}
e0 = st; e1 = ed;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
if(old_dp <= cov) {
rr = 1;
} else {
if(s0 < 0) {
s0 = st; s1 = ed;
}
e0 = st; e1 = ed;
}
}
if((s0 != e0) && (s1 != e1) && (s0 <= frank_len) && ((len - e1) <= frank_len) && (rr > 0)) {
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
if((old_dp <= cov) && (st >= s0) && (ed <= e1)) {
// if((ov->buffer[0].qns>>32) == 3364) fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
return 1;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
if((old_dp <= cov) && (st >= s0) && (ed <= e1)) {
// if((ov->buffer[0].qns>>32) == 3364) fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp);
return 1;
}
}
}
return 0;
}
uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t cut_len, double dup_rate)
{
uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, e0, rr, lt;
for (k = idx->n = 0; k < ov->length; k++) {
s0 = (uint32_t)ov->buffer[k].qns; e0 = ov->buffer[k].qe;
if(s0 > 0) s0 += cut_len;
if(e0 < len) e0 -= cut_len;
if(e0 <= s0) continue;
s = s0; e = e0;
lt = Get_READ_LENGTH((R_INF), ov->buffer[k].tn);
rr = (lt >= len)?(lt - len):(len - lt);
if((rr <= (len*dup_rate)) && (rr <= (lt*dup_rate)) && (ov->buffer[k].rev)) {
dp = (ov->buffer[k].qe) - ((uint32_t)ov->buffer[k].qns); dp = len - dp;
old_dp = ov->buffer[k].te - ov->buffer[k].ts; old_dp = lt - old_dp;
if((dp <= (len*dup_rate)) && (old_dp <= (lt*dup_rate))) continue;
}
kv_push(uint64_t, (*idx), (s<<1));
kv_push(uint64_t, (*idx), (e<<1)|1);
}
radix_sort_ec64(idx->a, idx->a + idx->n); s0 = e0 = rr = -1;
for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
old_dp = dp;
///if a[j] is qe
if (idx->a[k]&1) --dp;
else ++dp;
ed = idx->a[k]>>1;
if(ed > st) {
// if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) {
// fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32);
// }
if(old_dp <= cov) {
// if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len);
return 1;
}
}
st = ed;
}
ed = len; old_dp = dp;
if(ed > st) {
// if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) {
// fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32);
// }
if(old_dp <= cov) {
// if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len);
return 1;
}
}
return 0;
}
void prt_dbg_rid_paf(ma_hit_t_alloc *ov, UC_Read *ra, asg8_v *qa)
{
if(!(ov->length)) return;
uint64_t k, qn = (ov->buffer[0].qns>>32), qn_n, i, m; char *nn = NULL; FILE *fp = NULL; ma_hit_t *h = NULL;
qn_n = Get_NAME_LENGTH((R_INF), qn) + 64; MALLOC(nn, qn_n);
sprintf(nn, "%.*s.qry.fq", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w");
for (k = 0; k < ov->length; k++) {
i = ov->buffer[k].tn;
recover_UC_Read(ra, &R_INF, i);
fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq);
fprintf(fp, "+\n");
retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn);
for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1));
fprintf(fp, "\n");
}
fclose(fp);
sprintf(nn, "%.*s.ref.fq", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w");
i = qn;
recover_UC_Read(ra, &R_INF, i);
fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq);
fprintf(fp, "+\n");
retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn);
for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1));
fprintf(fp, "\n");
fclose(fp);
sprintf(nn, "%.*s.ref.fa", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w");
i = qn;
recover_UC_Read(ra, &R_INF, i);
fprintf(fp, ">%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq);
// fprintf(fp, "+\n");
// retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn);
// for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1));
// fprintf(fp, "\n");
fclose(fp);
sprintf(nn, "%.*s.ov.paf", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w");
for (k = 0; k < ov->length; k++) {
h = &(ov->buffer[k]);
fprintf(fp, "%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev],
(int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl);
}
fclose(fp);
free(nn);
}
static void worker_hap_dc_ec_chemical_r(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
ma_hit_t_alloc *paf = &(R_INF.paf[i]); uint64_t k, m;
// if (memcmp("3ed80bc4-1169-4948-a9ff-9c2463b7f7a2", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-a-[M::%s-beg] rid->%ld, b->rr->%lu\n", __func__, i, b->rr);
// }
if(b->cnt[1] == 0) {
// if(i == 6204620) prt_dbg_rid_paf(&(R_INF.paf[i]), &(b->self_read), &(b->v8q));
// if(is_chemical_r(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), 3, 16)) {
if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) {
// fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));
R_INF.paf[i].length = 0; b->cnt[0]++;
}
} else if(b->cnt[1] == 1) {
for (k = 0; k < paf->length; k++) {
if(R_INF.paf[paf->buffer[k].tn].length == 0) {
paf->buffer[k].tn = (uint32_t)-1; b->cnt[0]++;
}
}
} else {
for (k = m = 0; k < paf->length; k++) {
if(paf->buffer[k].tn == ((uint32_t)-1)) continue;
paf->buffer[m++] = paf->buffer[k];
}
paf->length = m;
}
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
static void worker_hap_dc_ec_chemical_arc(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
ma_hit_t_alloc *paf = &(R_INF.paf[i]); uint64_t k;
if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) {
// fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));
for (k = 0; k < paf->length; k++) paf->buffer[k].del = 1; b->cnt[0]++;
}
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
void gen_ovlst_paf(ma_hit_t_alloc *in_e, ma_hit_t_alloc *in_r, asg64_v *ou)
{
uint32_t n = 0, k;
@@ -3843,7 +4284,7 @@ overlap_region* h_ec_lchain_re1(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea
// fprintf(stderr, "-0-[M::%s]\n", __func__);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf); rs = qu->seq;
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq;
// fprintf(stderr, "-1-[M::%s]\n", __func__);
@@ -4507,7 +4948,7 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea
// fprintf(stderr, "-0-[M::%s]\n", __func__);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf); rs = qu->seq;
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq;
// fprintf(stderr, "-1-[M::%s]\n", __func__);
@@ -5588,7 +6029,55 @@ void sl_ec_r(uint64_t n_thre, uint64_t n_a)
kt_for(n_thre, worker_sl_ec, b, n_a);///debug_for_fix
for (k = 0; k < n_thre; k++) {
free(b[k].a); destory_UC_Read(&b[k].z);
free(b[k].a); destory_UC_Read(&b[k].z); kv_destroy(b[k].q);
}
free(b);
}
void handle_chemical_r(uint64_t n_thre, uint64_t n_a)
{
ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0, dedup = 0;
b = gen_ec_ovec_buf_t(n_thre);
for (k = 0; k < n_thre; ++k) {
b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 0;
}
kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a);
for (k = 0; k < n_thre; ++k) {
chem_n += b->a[k].cnt[0];
b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 1;
}
kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a);
for (k = 0; k < n_thre; ++k) {
dedup += b->a[k].cnt[0];
b->a[k].cnt[1] = 2;
}
kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a);
fprintf(stderr, "[M::%s] # chemical reads: %lu, # arcs:: %lu\n", __func__, chem_n, dedup);
destroy_ec_ovec_buf_t(b);
}
void handle_chemical_arc(uint64_t n_thre, uint64_t n_a)
{
ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0;
b = gen_ec_ovec_buf_t(n_thre);
for (k = 0; k < n_thre; ++k) {
b->a[k].cnt[0] = 0;
}
kt_for(n_thre, worker_hap_dc_ec_chemical_arc, b, n_a);
for (k = 0; k < n_thre; ++k) {
chem_n += b->a[k].cnt[0];
}
fprintf(stderr, "[M::%s] # chemical reads: %lu\n", __func__, chem_n);
destroy_ec_ovec_buf_t(b);
}
ecovlp.h
+2 -48
View File
@@ -7,58 +7,12 @@
#include "Process_Read.h"
#include "kdq.h"
KDQ_INIT(uint32_t)
typedef struct {
uint32_t v:31, f:1;
uint32_t sc;
} cns_arc;
typedef struct {size_t n, m, nou; cns_arc *a; } cns_arc_v;
typedef struct {
// uint16_t c:2, t:2, f:1, sc:3;
uint32_t c:2, f:1, sc:29;
cns_arc_v arc;
}cns_t;
typedef struct {
size_t n, m;
cns_t *a;
uint32_t si, ei, off, bn, bb0, bb1, cns_g_wl;
kdq_t(uint32_t) *q;
}cns_gfa;
typedef struct {
// chaining and overlapping related buffers
UC_Read self_read, ovlp_read;
Candidates_list clist;
overlap_region_alloc olist;
ha_abuf_t *ab;
// int64_t num_read_base, num_correct_base, num_recorrect_base;
uint64_t cnt[6], rr;
haplotype_evdience_alloc hap;
bit_extz_t exz;
kv_ul_ov_t pidx;
asg64_v v64;
asg32_v v32;
asg16_v v16;
asg8_v v8q, v8t;
kvec_t_u8_warp k_flag;
st_mt_t sp;
cns_gfa cns;
} ec_ovec_buf_t0;
typedef struct {
ec_ovec_buf_t0 *a;
uint32_t n, rev;
} ec_ovec_buf_t;
ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n);
void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p);
void prt_chain(overlap_region_alloc *o);
void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, uint64_t is_sv, uint64_t *tot_b, uint64_t *tot_e);
void sl_ec_r(uint64_t n_thre, uint64_t n_a);
void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx);
void handle_chemical_r(uint64_t n_thre, uint64_t n_a);
void handle_chemical_arc(uint64_t n_thre, uint64_t n_a);
#endif
gfa_ut.cpp
+180 -1
View File
@@ -1407,6 +1407,173 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_
if (cnt > 0) asg_cleanup(g);
}
uint32_t is_dedup_weak_arc(asg_arc_t *av, uint32_t an, uint32_t ak, ma_hit_t_alloc *rev)
{
uint32_t k, w, m; ma_hit_t_alloc *z = NULL;
for (k = 0; k < an; k++) {
if((k == ak) || (av[k].del) || (!av[k].strong)) continue;
if(av[k].ol > av[ak].ol) {
z = &(rev[av[k].v>>1]); w = (av[ak].v^(av[k].v&1));
for (m = 0; m < z->length; m++) {
if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == w) break;
}
if(m < z->length) return 1;
}
}
return 0;
}
void asg_arc_cut_weak(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio,
uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, uint32_t test_bub, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t i, k, v, w, wz, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou, olw[2], m;
asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max; ma_hit_t_alloc *z = NULL;
if(in) b = in;
else b = &tx;
b->n = 0;
for (v = 0; v < n_vtx; ++v) {
// if((v>>1)==17078) fprintf(stderr, "[M::%s::] v:%u, del:%u, seq_vis:%u\n", __func__, v, g->seq[v>>1].del, g->seq_vis[v]);
if (g->seq[v>>1].del) continue;
if((test_bub == 0) || (g->seq_vis[v] == 0)) {
av = asg_arc_a(g, v); nv = asg_arc_n(g, v);
if (nv < 2) continue;
for (i = kv = 0; i < nv; ++i) {
if(av[i].del) continue;
kv++;
}
if(kv < 2) continue;
for (i = olw[1] = 0, olw[0] = (uint32_t)-1; i < nv; ++i) {
if(av[i].del) continue;
if(av[i].strong) {
if(av[i].ol > olw[av[i].strong]) olw[av[i].strong] = av[i].ol;
} else {
if(av[i].ol < olw[av[i].strong]) olw[av[i].strong] = av[i].ol;
}
}
if(olw[1] <= olw[0]) continue;
for (i = 0; i < nv; ++i) {
if(av[i].del || av[i].strong) continue;
if(av[i].ol >= olw[1]) continue;
if(max_drop_len && av[i].ol >= (*max_drop_len)) continue;
if(is_dedup_weak_arc(av, nv, i, rev)) {
kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i)));
}
}
}
}
if(rev && rI) memset(g->seq_vis, 0, g->n_seq*2*sizeof(uint8_t));
radix_sort_srt64(b->a, b->a + b->n);
for (k = 0; k < b->n; k++) {
if(g->arc[(uint32_t)b->a[k]].del) continue;
v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1;
if(g->seq[v>>1].del || g->seq[w>>1].del) continue;
nv = asg_arc_n(g, v); nw = asg_arc_n(g, w);
av = asg_arc_a(g, v); aw = asg_arc_a(g, w);
if(nv<=1 && nw <= 1) continue;
if(is_trio) {
if(get_arcs(g, v, NULL, 0)<=1 && get_arcs(g, w, NULL, 0)<=1) continue;///speedup
trioF = get_tip_trio_infor(g, v^1);
ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1));
}
ve = &(g->arc[(uint32_t)b->a[k]]);
for (i = 0; i < nw; ++i) {
if (aw[i].v == (v^1)) {
we = &(aw[i]);
break;
}
}
///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted
mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou);
for (i = kv = ol_max = ou_max = 0, vl_max = NULL; i < nv; ++i) {
if(av[i].del) continue;
kv++;
if((av[i].v == ve->v) || (!av[i].strong) || (av[i].ol <= ve->ol)) continue;
if(is_trio && get_tip_trio_infor(g, av[i].v) == ntrioF) continue;
z = &(rev[av[i].v>>1]); wz = (ve->v^(av[i].v&1));
for (m = 0; m < z->length; m++) {
if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == wz) break;
}
if(m >= z->length) continue;
if(ol_max < av[i].ol) ol_max = av[i].ol, vl_max = &(av[i]);
if(ou_max < av[i].ou) ou_max = av[i].ou;
}
if (kv < 1) continue;
if (kv >= 2) {
if (mm_ol > ol_max*len_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
if ((mm_ol + min_diff) > ol_max) continue;
}
for (i = kw = ol_max = ou_max = 0, wl_max = NULL; i < nw; ++i) {
if(aw[i].del) continue;
kw++;
if((aw[i].v == we->v) || (!aw[i].strong) || (aw[i].ol <= we->ol)) continue;
if(is_trio && get_tip_trio_infor(g, aw[i].v) == ntrioF) continue;
z = &(rev[aw[i].v>>1]); wz = (we->v^(aw[i].v&1));
for (m = 0; m < z->length; m++) {
if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == wz) break;
}
if(m >= z->length) continue;
if(ol_max < aw[i].ol) ol_max = aw[i].ol, wl_max = &(aw[i]);
if(ou_max < aw[i].ou) ou_max = aw[i].ou;
}
if (kw < 1) continue;
if (kw >= 2) {
if (mm_ol > ol_max*len_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
if ((mm_ol + min_diff) > ol_max) continue;
}
if (kv <= 1 && kw <= 1) continue;
to_del = 0;
if(is_topo) {
if (kv > 1 && kw > 1) {
to_del = 1;
} else if (kw == 1) {
if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1;
} else if (kv == 1) {
if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1;
}
}
if(rev && rI) {
if((to_del == 0) && vl_max && (ve->v!=vl_max->v) && (trans_path_check(ve->v, vl_max->v, g, rev, rI, max_ext, b)==0)) {
to_del = 1;
}
if((to_del == 0) && wl_max && (we->v!=wl_max->v) && (trans_path_check(we->v, wl_max->v, g, rev, rI, max_ext, b)==0)) {
to_del = 1;
}
if(vl_max && wl_max) assert(ve->v!=vl_max->v||we->v!=wl_max->v);
}
if (to_del) {
ve->del = we->del = 1, ++cnt;
}
}
// stats_sysm(g);
if(!in) free(tx.a);
if (cnt > 0) asg_cleanup(g);
}
void asg_arc_cut_length_adv(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio,
uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
@@ -2784,7 +2951,9 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819));
// debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0");
// debug_info_of_specfic_node("m64012_190920_173625/163644465/ccs", sg, rI, "beg-0");
// debug_info_of_specfic_node("bcb40bcc-d9cf-48e6-88ee-47ac3dde22ff", sg, rI, "beg-0");
if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo
// fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847));
@@ -2795,6 +2964,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(drop <= 0.500001) min_diff = step_diff>>1;
else min_diff = step_diff;
}
if(asm_opt.is_ont) {
asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
}
// fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop);
// prt_specfic_sge(sg, 10531, 10519, "--0--");
@@ -2843,6 +3017,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
}
}
if(asm_opt.is_ont) {
asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
}
if(is_ou) min_diff = step_diff;
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty4.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "end");