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

has bug multiple threads

Browse Source
This commit is contained in:
chhylp123
2021-12-28 21:48:41 -05:00
parent 8a0ec3d575
commit dbc72b8c2d
5 changed files with 447 additions and 160 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+7 -25
View File
@@ -31102,9 +31102,14 @@ char *get_outfile_name(char* output_file_name)
return buf; return buf;
} }
void create_ul_info() void create_ul_info(ma_hit_t_alloc* sources, int max_hang, int min_ovlp, long long gap_fuzz)
{ {
ug_opt_t opt; memset(&opt, 0, sizeof(opt)); ug_opt_t opt; memset(&opt, 0, sizeof(opt));
opt.sources = sources;
opt.max_hang = max_hang;
opt.min_ovlp = min_ovlp;
opt.gap_fuzz = gap_fuzz;
ul_load(&opt); ul_load(&opt);
exit(1); exit(1);
} }
@@ -31117,30 +31122,6 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
ma_sub_t **coverage_cut_ptr, int debug_g) ma_sub_t **coverage_cut_ptr, int debug_g)
{ {
if(asm_opt.ar) create_ul_info();
char *o_file = get_outfile_name(output_file_name); char *o_file = get_outfile_name(output_file_name);
ma_sub_t *coverage_cut = *coverage_cut_ptr; ma_sub_t *coverage_cut = *coverage_cut_ptr;
asg_t *sg = *sg_ptr; asg_t *sg = *sg_ptr;
@@ -31166,6 +31147,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
{ {
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
} }
if(asm_opt.ar) create_ul_info(sources, max_hang_length, mini_overlap_length, gap_fuzz);
///print_binned_reads(sources, n_read, coverage_cut); ///print_binned_reads(sources, n_read, coverage_cut);
clean_weak_ma_hit_t(sources, reverse_sources, n_read); clean_weak_ma_hit_t(sources, reverse_sources, n_read);
Process_Read.cpp
+404 -102
View File
@@ -28,7 +28,7 @@ uint8_t seq_nt6_table[256] = {
char bit_t_seq_table[256][4] = {{0}}; char bit_t_seq_table[256][4] = {{0}};
char bit_t_seq_table_rc[256][4] = {{0}}; char bit_t_seq_table_rc[256][4] = {{0}};
char s_H[5] = {'A', 'C', 'G', 'T', 'N'}; char s_H[5] = {'A', 'C', 'G', 'T', 'N'};
char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'}; char rc_Table[6] = {'T', 'G', 'C', 'A', 'N', 'N'};
void init_All_reads(All_reads* r) void init_All_reads(All_reads* r)
{ {
@@ -341,87 +341,92 @@ void init_UC_Read(UC_Read* r)
} }
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) void recover_UC_Read_sub_region(char* r, int64_t start_pos, int64_t length, uint8_t strand, All_reads* R_INF, int64_t ID)
{ {
long long readLen = Get_READ_LENGTH((*R_INF), ID); int64_t readLen = Get_READ_LENGTH((*R_INF), ID);
int64_t end_pos = start_pos + length - 1, begLen, tailLen, offset, mn, src_i, des_i, i;
uint8_t* src = Get_READ((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID);
long long i; if(strand == 0) {
long long copyLen; offset = start_pos&3;
long long end_pos = start_pos + length - 1; begLen = 4-offset;
if(begLen > length) begLen = length;
tailLen = (length-begLen)&3;
mn = (length - begLen - tailLen)>>2;
src_i = start_pos; des_i = 0; i = 0;
if (strand == 0) if(begLen > 0) {
{ memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]]+offset, begLen);
i = start_pos; des_i += begLen; src_i += begLen;
copyLen = 0;
long long initLen = start_pos % 4;
if (initLen != 0)
{
memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
} }
while (copyLen < length) for (i = 0; i < mn; i++) {
{ memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], 4);
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); des_i += 4; src_i += 4;
copyLen = copyLen + 4;
i = i + 4;
} }
if (R_INF->N_site[ID]) if(tailLen > 0) {
{ memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], tailLen);
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) des_i += tailLen; src_i += tailLen;
{ }
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
{ if (R_INF->N_site[ID]) {
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) {
if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos) {
r[R_INF->N_site[ID][i] - start_pos] = 'N'; r[R_INF->N_site[ID][i] - start_pos] = 'N';
} }
else if((long long)R_INF->N_site[ID][i] > end_pos) else if((long long)R_INF->N_site[ID][i] > end_pos) {
{
break; break;
} }
} }
} }
} }
else else {
{
start_pos = readLen - start_pos - 1; start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1; end_pos = readLen - end_pos - 1;
///start_pos > end_pos begLen = (start_pos+1)&3;
i = start_pos; offset = 4 - begLen;
copyLen = 0; if(begLen > length) begLen = length;
long long initLen = (start_pos + 1) % 4; tailLen = (length-begLen)&3;
mn = (length - begLen - tailLen)>>2;
src_i = start_pos; des_i = 0; i = 0;
if (initLen != 0) if(begLen > 0) {
{ memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]]+offset, begLen);
memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen); des_i += begLen; src_i -= begLen;
copyLen = copyLen + initLen;
i = i - initLen;
} }
while (copyLen < length) for (i = 0; i < mn; i++) {
{ memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], 4);
memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4); des_i += 4; src_i -= 4;
copyLen = copyLen + 4;
i = i - 4;
} }
if (R_INF->N_site[ID]) if(tailLen > 0) {
{ memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], tailLen);
long long offset = readLen - start_pos - 1; des_i += tailLen; src_i -= tailLen;
}
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) /**
{ if (R_INF->N_site[ID]) {
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) offset = readLen - start_pos - 1;
{ for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) {
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) {
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
} }
else if((long long)R_INF->N_site[ID][i] > start_pos) else if((long long)R_INF->N_site[ID][i] > start_pos) {
{ break;
}
}
}
**/
if (R_INF->N_site[ID]) {
start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1;
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) {
offset = readLen - R_INF->N_site[ID][i] - 1;
if(offset >= start_pos && offset <= end_pos) {
r[offset - start_pos] = 'N';
} else if(offset < start_pos) {
break; break;
} }
} }
@@ -676,6 +681,7 @@ void reverse_complement(char* pattern, uint64_t length)
for (i = 0; i < end; i++) for (i = 0; i < end; i++)
{ {
index = length - i - 1; index = length - i - 1;
k = pattern[index]; k = pattern[index];
pattern[index] = RC_CHAR(pattern[i]); pattern[index] = RC_CHAR(pattern[i]);
@@ -743,8 +749,9 @@ void destory_Debug_reads(Debug_reads* x)
void init_all_ul_t(all_ul_t *x, All_reads *hR) { void init_all_ul_t(all_ul_t *x, All_reads *hR) {
memset(x, 0, sizeof(*x)); memset(x, 0, sizeof(*x));
x->hR = hR; x->mm = 0x7fffffff; x->hR = hR; x->mm = 0x7fffffff;
init_aux_table();
} }
void destory_all_ul_t(all_ul_t *x, All_reads *hR) { void destory_all_ul_t(all_ul_t *x) {
uint64_t i; uint64_t i;
for (i = 0; i < x->n; i++) { for (i = 0; i < x->n; i++) {
free(x->a[i].n_n); free(x->a[i].N_site.a); free(x->a[i].n_n); free(x->a[i].N_site.a);
@@ -753,7 +760,7 @@ void destory_all_ul_t(all_ul_t *x, All_reads *hR) {
free(x->a); free(x->a);
} }
void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset)
{ {
uint64_t i = 0; uint64_t i = 0;
uint64_t dest_i = 0; uint64_t dest_i = 0;
@@ -765,28 +772,28 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn)
c = seq_nt6_table[(uint8_t)src[i]]; c = seq_nt6_table[(uint8_t)src[i]];
if (c >= 4) { if (c >= 4) {
c = 0; kv_push(uint32_t, *nn, i); c = 0; kv_push(uint32_t, *nn, i+nn_offset);
} }
i++; i++;
tmp = tmp | (c<<6); tmp = tmp | (c<<6);
c = seq_nt6_table[(uint8_t)src[i]]; c = seq_nt6_table[(uint8_t)src[i]];
if (c >= 4) { if (c >= 4) {
c = 0; kv_push(uint32_t, *nn, i); c = 0; kv_push(uint32_t, *nn, i+nn_offset);
} }
i++; i++;
tmp = tmp | (c<<4); tmp = tmp | (c<<4);
c = seq_nt6_table[(uint8_t)src[i]]; c = seq_nt6_table[(uint8_t)src[i]];
if (c >= 4) { if (c >= 4) {
c = 0; kv_push(uint32_t, *nn, i); c = 0; kv_push(uint32_t, *nn, i+nn_offset);
} }
i++; i++;
tmp = tmp | (c<<2); tmp = tmp | (c<<2);
c = seq_nt6_table[(uint8_t)src[i]]; c = seq_nt6_table[(uint8_t)src[i]];
if (c >= 4) { if (c >= 4) {
c = 0; kv_push(uint32_t, *nn, i); c = 0; kv_push(uint32_t, *nn, i+nn_offset);
} }
i++; i++;
tmp = tmp | c; tmp = tmp | c;
@@ -803,7 +810,7 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn)
while (i < src_l) { while (i < src_l) {
c = seq_nt6_table[(uint8_t)src[i]]; c = seq_nt6_table[(uint8_t)src[i]];
if (c >= 4) { if (c >= 4) {
c = 0; kv_push(uint32_t, *nn, i); c = 0; kv_push(uint32_t, *nn, i+nn_offset);
} }
i++; i++;
tmp = tmp | (c << shift); tmp = tmp | (c << shift);
@@ -816,14 +823,33 @@ void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn)
} }
#define B4L(x) (((x)>>2)+(((x)&3)?1:0)) #define B4L(x) (((x)>>2)+(((x)&3)?1:0))
void push_subblock_original_bases(char* str, all_ul_t *x, ul_vec_t *p, uint32_t s, uint32_t e, uint32_t subLen)//for debug
{
uc_block_t *b = NULL;
uint32_t qs = s, qe = e;
while (qs < e) {
qe = qs + subLen; if(qe > e) qe = e;
kv_pushp(uc_block_t, p->bb, &b);
b->hid = x->mm; b->rev = 0;
b->qs = qs; b->qe = qe;
b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs);
kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te;
ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs);
qs = qe;
}
}
void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) { void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) {
int64_t i, mine, maxs, ovlp, end; int64_t i, mine, maxs, ovlp, end;
ul_vec_t *p = NULL; ul_vec_t *p = NULL;
ul_ov_t *z = NULL, *zp = NULL; ul_ov_t *z = NULL, *zp = NULL;
uc_block_t *b = NULL; uc_block_t *b = NULL;
if(rid) fprintf(stderr, "rid:%lu\n", *rid);
if(rid == NULL) { if(rid == NULL) {
kv_pushp(ul_vec_t, *x, &p); kv_pushp(ul_vec_t, *x, &p);
memset(p, 0, sizeof(*p)); memset(p, 0, sizeof(*p));
fprintf(stderr, "x->n:%u\n", x->n);
} }
else { else {
p = &(x->a[*rid]); p = &(x->a[*rid]);
@@ -854,7 +880,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
b->qs = end; b->qe = b->qs - ovlp; b->qs = end; b->qe = b->qs - ovlp;
b->ts = p->r_base.n; b->te = b->ts + B4L(-ovlp); b->ts = p->r_base.n; b->te = b->ts + B4L(-ovlp);
kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te;
ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs);
} }
///push ovlp bases ///push ovlp bases
@@ -876,72 +902,348 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
b->qs = end; b->qe = str_l; b->qs = end; b->qe = str_l;
b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs); b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs);
kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te;
ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site), b->qs);
// push_subblock_original_bases(str, x, p, end, str_l, 321);//for debug
} }
// char *sst = NULL; CALLOC(sst, str_l);//for debug
// retrieve_ul_t(NULL, sst, x, rid?*rid:x->n-1, 0, 0, -1);
// if(memcmp(sst, str, str_l)) {
// fprintf(stderr, "ap-Wrong read, id: %ld, [%d, %ld)\n", (int64_t)(rid?*rid:x->n-1), 0, str_l);
// for (i = 0; i < str_l; i++) {
// if(sst[i] != str[i]) fprintf(stderr, "[%ld] input:%c, decompress:%c\n", i, str[i], sst[i]);
// }
// }
// free(sst);
} }
} }
void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) { void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l) {
ul_vec_t *p = &(ref->a[ID]); ul_vec_t *p = &(ref->a[ID]);
if(l < 0) l = p->rlen;
uc_block_t *b = NULL; uc_block_t *b = NULL;
char *a = NULL; char *r = NULL;
uint8_t *src = NULL; uint8_t *src = NULL;
int64_t k, i, a_n, l_chr, idx; int64_t k, i, a_n, e = s + l, ssp, sep, sl, rts, rte;
r->length = p->rlen; r->RID = ID; int64_t offset, begLen, tailLen, src_i, des_i;
if(i_r) {
i_r->length = l; i_r->RID = ID;
if(i_r->length > i_r->size) {
i_r->size = i_r->length;
i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
}
r = i_r->seq;
}
if(i_s) r = i_s;
if(r->length + 4 > r->size) {
r->size = r->length + 4;
r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size));
}
a = NULL;
if(strand == 0) { if(strand == 0) {
for (k = 0; k < (int64_t)p->bb.n; k++) { for (k = 0, des_i = 0; k < (int64_t)p->bb.n; k++) {
b = &(p->bb.a[k]); a = r->seq + b->qs; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; b = &(p->bb.a[k]);
if(b->qe <= s) continue;
if(b->qs >= e) break;
src = p->r_base.a + b->ts;
ssp = MAX(s, b->qs) - b->qs;
sep = MIN(e, b->qe) - b->qs;
sl = sep - ssp;
if(b->hid == ref->mm){///original bases if(b->hid == ref->mm){///original bases
i = 0; offset = ssp&3;
while (i < a_n) { begLen = 4-offset;
memcpy(a+i, bit_t_seq_table[src[i>>2]], 4); if(begLen > sl) begLen = sl;
i += 4; tailLen = (sl-begLen)&3;
a_n = (sl - begLen - tailLen)>>2;
src_i = ssp; i = 0;
if(begLen > 0) {
memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]]+offset, begLen);
des_i += begLen; src_i += begLen;
}
for (i = 0; i < a_n; i++) {
memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], 4);
des_i += 4; src_i += 4;
}
if(tailLen > 0) {
memcpy(r+des_i, bit_t_seq_table[src[src_i>>2]], tailLen);
des_i += tailLen; src_i += tailLen;
} }
} else {///ovlps } else {///ovlps
recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev, ref->hR, b->hid); if(b->rev == 0) {
recover_UC_Read_sub_region(r+des_i, b->ts + ssp, sep - ssp, b->rev, ref->hR, b->hid);
} else{
rts = b->ts + (b->qe - sep); rte = rts + sep - ssp;
recover_UC_Read_sub_region(r+des_i, Get_READ_LENGTH((*ref->hR), b->hid) - rte,
sep - ssp, b->rev, ref->hR, b->hid);
}
des_i += sep - ssp;
}
}
for (k = 0; k < (int64_t)p->N_site.n; k++) {
if(p->N_site.a[k] >= s && p->N_site.a[k] < e){
r[p->N_site.a[k]-s] = 'N';
}
else if(p->N_site.a[k] >= e) {
break;
} }
} }
for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[p->N_site.a[k]] = 'N';
} else { } else {
for (k = 0; k < (int64_t)p->bb.n; k++) { sep = p->rlen - s;
b = &(p->bb.a[k]); a = r->seq + r->length - b->qe; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; ssp = p->rlen - e;
s = ssp; e = sep;
///[s, e)
for (k = ((int64_t)p->bb.n)-1, des_i = 0; k >= 0; k--) {
b = &(p->bb.a[k]);
if(b->qe <= s) break;
if(b->qs >= e) continue;
src = p->r_base.a + b->ts;
ssp = MAX(s, b->qs) - b->qs;
sep = MIN(e, b->qe) - b->qs;
sl = sep - ssp;
///[ssp, sep)
if(b->hid == ref->mm){///original bases if(b->hid == ref->mm){///original bases
l_chr = (b->qe - b->qs)&((uint32_t)3); begLen = sep&3;
i = ((b->qe - b->qs)>>2)+(l_chr!= 0)-1; offset = 4 - begLen;
idx = 0; if(begLen > sl) begLen = sl;
if(l_chr > 0) { tailLen = (sl-begLen)&3;
memcpy(a + idx, bit_t_seq_table_rc[src[i]]+4-l_chr, l_chr); a_n = (sl - begLen - tailLen)>>2;
idx = l_chr; i--; src_i = sep-1; i = 0;
if(begLen > 0) {
memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]]+offset, begLen);
des_i += begLen; src_i -= begLen;
} }
while (i >= 0) {
memcpy(a + idx, bit_t_seq_table_rc[src[i]], 4); for (i = 0; i < a_n; i++) {
i--; idx += 4; memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], 4);
des_i += 4; src_i -= 4;
}
if(tailLen > 0) {
memcpy(r+des_i, bit_t_seq_table_rc[src[src_i>>2]], tailLen);
des_i += tailLen; src_i -= tailLen;
} }
} else {///ovlps } else {///ovlps
recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev==strand?0:1, ref->hR, b->hid); if(b->rev == 0) {///b->rev != strand
rts = b->ts + ssp; rte = rts + sep - ssp;
recover_UC_Read_sub_region(r+des_i, Get_READ_LENGTH((*ref->hR), b->hid) - rte, sep - ssp, 1, ref->hR, b->hid);
} else {///b->rev == strand
rts = b->ts + (b->qe - sep); rte = rts + sep - ssp;
recover_UC_Read_sub_region(r+des_i, rts, sep - ssp, 0, ref->hR, b->hid);
}
des_i += sep - ssp;
} }
} }
for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[r->length - p->N_site.a[k] - 1] = 'N';
sep = p->rlen - s;
ssp = p->rlen - e;
s = ssp; e = sep;
for (k = 0; k < (int64_t)p->N_site.n; k++) {
sl = p->rlen - p->N_site.a[k] - 1;
if(sl >= s && sl < e) r[sl-s] = 'N';
else if(sl < s) {
break;
}
}
} }
} }
void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand)
void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l)
{ {
if(u->m == 0 || u->n == 0) return;
if(l < 0) l = u->len;
char *r = NULL, *a = NULL;
int64_t e = s + l, ssp, sep, rs, re, des_i;
uint64_t k, rId, ori, r_l;
if(i_r) { if(i_r) {
i_r->length = l; i_r->RID = 0;
if(i_r->length > i_r->size) {
i_r->size = i_r->length;
i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
}
r = i_r->seq;
}
if(i_s) r = i_s;
if(strand == 1) {
sep = u->len - s;
ssp = u->len - e;
s = ssp; e = sep;
}
for (k = l = des_i = 0; k < u->n; k++) {
rId = u->a[k]>>33;
ori = u->a[k]>>32&1;
r_l = (uint32_t)u->a[k];
if(r_l == 0) continue;
ssp = l; sep = l + r_l;
l += r_l;
if(sep <= s) continue;
if(ssp >= e) break;
rs = MAX(ssp, s); re = MIN(sep, e);
a = r + des_i; des_i += re - rs;
recover_UC_Read_sub_region(a, rs-ssp, re-rs, ori, &R_INF, rId);
}
if(strand == 1) {
char t;
re = (e - s);
l = re>>1;
for (k = 0; k < (uint64_t)l; k++) {
des_i = re - k - 1;
t = r[des_i];
r[des_i] = RC_CHAR(r[k]);
r[k] = RC_CHAR(t);
}
if(re&1) r[l] = RC_CHAR(r[l]);
}
}
void produce_u_seq(char* r, ma_utg_t *u, UC_Read *buf)
{
if(u->m == 0 || u->n == 0) return;
uint32_t j, k, l = 0;
uint32_t rId, ori, start, eLen, readLen;
char *readS = NULL;
memset(r, 'N', u->len);
for (j = 0; j < u->n; ++j) {
rId = u->a[j]>>33;
///uId = i;
ori = u->a[j]>>32&1;
start = l;
eLen = (uint32_t)u->a[j];
l += eLen;
if(eLen == 0) continue;
recover_UC_Read(buf, &R_INF, rId);
readS = buf->seq;
readLen = Get_READ_LENGTH(R_INF, rId);
if (!ori) // forward strand
{
for (k = 0; k < eLen; k++)
{
r[start + k] = readS[k];
}
}
else
{
for (k = 0; k < eLen; k++)
{
uint8_t c = (uint8_t)readS[readLen - 1 - k];
r[start + k] = c >= 128? 'N' : RC_CHAR(c);
}
}
}
}
void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step)
{
uint64_t i, step, s, e, occ;
UC_Read f, r;
init_UC_Read(&f); init_UC_Read(&r);
kvec_t(char) ss; kv_init(ss);
if(ref) {
for (i = 0, occ = 0; i < ref->n; i++) {
retrieve_ul_t(&f, NULL, ref, i, 0, 0, -1);
retrieve_ul_t(&r, NULL, ref, i, 1, 0, -1);
kv_resize(char, ss, ref->a[i].rlen);
ss.n = ref->a[i].rlen;
memcpy(ss.a, r.seq, ss.n);
reverse_complement(ss.a, ss.n);
if(memcmp(ss.a, f.seq, ss.n)) {
fprintf(stderr, "1-Wrong whole reverse-read, id: %lu\n", i);
// for (s = 0; s < ss.n; s++) {
// if(ss.a[s] != f.seq[s]) {
// fprintf(stderr,"s:%lu, ss.a[s]:%c, f.seq[s]:%c, r.seq[rs]:%c\n", s, ss.a[s], f.seq[s], r.seq[ref->a[i].rlen - s - 1]);
// }
// }
}
step = ss.n/n_step;
if(step <= 0) step = 1;
for (s = 0; s < ss.n; s += step) {
e = MIN(s+step, ss.n);
retrieve_ul_t(NULL, ss.a, ref, i, 0, s, e-s);
if(memcmp(ss.a, f.seq + s, e - s)) {
fprintf(stderr, "1-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e);
}
retrieve_ul_t(NULL, ss.a, ref, i, 1, s, e-s);
if(memcmp(ss.a, r.seq + s, e - s)) {
fprintf(stderr, "1-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e);
// uint64_t dk; char *tf = ss.a; char *rf = r.seq + s;
// for (dk = 0; dk < e - s; dk++) {
// if(tf[dk] != rf[dk]) {
// fprintf(stderr,"dk:%lu, tf[dk]:%c, rf[dk]:%c\n",
// dk, tf[dk], rf[dk]);
// }
// }
}
occ++;
}
}
fprintf(stderr, "[M::%s::# checking: %lu] ==> all_ul_t\n", __func__, occ);
} }
if (r->length + 4 > r->size) if(R_INF) {
{ for (i = 0, occ = 0; i < R_INF->total_reads; i++) {
r->size = r->length + 4; recover_UC_Read(&f, R_INF, i);
r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); recover_UC_Read_RC(&r, (All_reads*)R_INF, i);
kv_resize(char, ss, Get_READ_LENGTH((*R_INF), i));
ss.n = Get_READ_LENGTH((*R_INF), i);
memcpy(ss.a, r.seq, ss.n);
reverse_complement(ss.a, ss.n);
if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "2-Wrong whole reverse-read, id: %lu\n", i);
step = ss.n/n_step;
if(step <= 0) step = 1;
for (s = 0; s < ss.n; s += step) {
e = MIN(s+step, ss.n);
recover_UC_Read_sub_region(ss.a, s, e-s, 0, (All_reads*)R_INF, i);
if(memcmp(ss.a, f.seq + s, e - s)) fprintf(stderr, "2-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e);
recover_UC_Read_sub_region(ss.a, s, e-s, 1, (All_reads*)R_INF, i);
if(memcmp(ss.a, r.seq + s, e - s)) fprintf(stderr, "2-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e);
occ++;
}
}
fprintf(stderr, "[M::%s::# checking: %lu] ==> All_reads\n", __func__, occ);
} }
if(u) {
for (i = 0, occ = 0; i < u->n; i++) {
kv_resize(char, ss, u->a[i].len); ss.n = u->a[i].len;
retrieve_u_seq(&f, NULL, &(u->a[i]), 0, 0, -1);
produce_u_seq(ss.a, &(u->a[i]), &r);
if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "4-Wrong whole reverse-read, id: %lu\n", i);
retrieve_u_seq(&r, NULL, &(u->a[i]), 1, 0, -1);
memcpy(ss.a, r.seq, ss.n);
reverse_complement(ss.a, ss.n);
if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "3-Wrong whole reverse-read, id: %lu\n", i);
step = ss.n/n_step;
if(step <= 0) step = 1;
for (s = 0; s < ss.n; s += step) {
e = MIN(s+step, ss.n);
retrieve_u_seq(NULL, ss.a, &(u->a[i]), 0, s, e-s);
if(memcmp(ss.a, f.seq + s, e - s)) fprintf(stderr, "3-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e);
retrieve_u_seq(NULL, ss.a, &(u->a[i]), 1, s, e-s);
if(memcmp(ss.a, r.seq + s, e - s)) fprintf(stderr, "3-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e);
occ++;
}
}
fprintf(stderr, "[M::%s::# checking: %lu] ==> ma_utg_v\n", __func__, occ);
}
destory_UC_Read(&f); destory_UC_Read(&r);
kv_destroy(ss);
} }
Process_Read.h
+6 -4
View File
@@ -29,7 +29,7 @@ extern uint8_t seq_nt6_table[256];
extern char bit_t_seq_table[256][4]; extern char bit_t_seq_table[256][4];
extern char bit_t_seq_table_rc[256][4]; extern char bit_t_seq_table_rc[256][4];
extern char s_H[5]; extern char s_H[5];
extern char rc_Table[5]; extern char rc_Table[6];
#define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]]
@@ -198,7 +198,7 @@ void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_sit
void init_UC_Read(UC_Read* r); void init_UC_Read(UC_Read* r);
void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID); void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID);
void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID);
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); void recover_UC_Read_sub_region(char* r, int64_t start_pos, int64_t length, uint8_t strand, All_reads* R_INF, int64_t ID);
void destory_UC_Read(UC_Read* r); void destory_UC_Read(UC_Read* r);
void reverse_complement(char* pattern, uint64_t length); void reverse_complement(char* pattern, uint64_t length);
void write_All_reads(All_reads* r, char* read_file_name); void write_All_reads(All_reads* r, char* read_file_name);
@@ -210,8 +210,10 @@ void destory_Debug_reads(Debug_reads* x);
void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID);
void init_all_ul_t(all_ul_t *x, All_reads *hR); void init_all_ul_t(all_ul_t *x, All_reads *hR);
void destory_all_ul_t(all_ul_t *x, All_reads *hR); void destory_all_ul_t(all_ul_t *x);
void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on); void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on);
void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand); void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l);
void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l);
void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step);
#endif #endif
htab.cpp
+1 -1
View File
@@ -733,7 +733,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
}\ }\
MALLOC(s->seq[s->n_seq], u->len);\ MALLOC(s->seq[s->n_seq], u->len);\
if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\ if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\
else memcpy(s->seq[s->n_seq], u->s, u->len);\ else retrieve_u_seq(NULL, s->seq[s->n_seq], u, 0, 0, -1);\
s->len[s->n_seq++] = u->len;\ s->len[s->n_seq++] = u->len;\
++p->n_seq;\ ++p->n_seq;\
s->sum_len += u->len;\ s->sum_len += u->len;\
inter.cpp
+29 -28
View File
@@ -288,10 +288,10 @@ typedef struct { // data structure for each step in kt_pipeline()
int n, m, sum_len; int n, m, sum_len;
uint64_t *len, id; uint64_t *len, id;
char **seq; char **seq;
ha_mzl_v *mzs; ha_mzl_v *mzs;///useless
st_mt_t *sps; st_mt_t *sps;///useless
mg_gchains_t **gcs; mg_gchains_t **gcs;///useless
mg_tbuf_t **buf; mg_tbuf_t **buf;///useless
ha_ovec_buf_t **hab; ha_ovec_buf_t **hab;
} utepdat_t; } utepdat_t;
@@ -2320,6 +2320,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
} }
**/ **/
for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]); for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]);
free(s->hab);
// free(s->buf); free(s->mzs); free(s->sps); // free(s->buf); free(s->mzs); free(s->sps);
return s; return s;
} }
@@ -2337,13 +2338,14 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
**/ **/
rid = s->id + i; rid = s->id + i;
append_ul_t(&UL_INF, &rid, NULL, 0, s->seq[i], s->len[i], NULL, 0); append_ul_t(&UL_INF, &rid, NULL, 0, s->seq[i], s->len[i], NULL, 0);
// fprintf(stderr, "%.*s\n", (int)s->len[i], s->seq[i]);
free(s->seq[i]); p->total_base += s->len[i]; free(s->seq[i]); p->total_base += s->len[i];
} }
///debug ///debug
/** /**
free(s->gcs); free(s->gcs);
**/ **/
free(s); free(s->len); free(s->seq); free(s);
} }
return 0; return 0;
} }
@@ -2380,7 +2382,7 @@ int print_ul_rs(all_ul_t *U_INF)
init_UC_Read(&ur); init_UC_Read(&ur);
for (i = 0; i < U_INF->n; i++) { for (i = 0; i < U_INF->n; i++) {
p = &(U_INF->a[i]); p = &(U_INF->a[i]);
retrieve_ul_t(&ur, U_INF, i, 0); retrieve_ul_t(&ur, NULL, U_INF, i, 0, 0, -1);
fprintf(stderr, ">%s\n", p->n_n); fprintf(stderr, ">%s\n", p->n_n);
fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq); fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq);
} }
@@ -2870,19 +2872,19 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n)
uidx_destory(); uidx_destory();
} }
int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug) int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug)
{ {
uldat_t sl; memset(&sl, 0, sizeof(sl)); uldat_t sl; memset(&sl, 0, sizeof(sl));
sl.ha_flt_tab = ha_flt_tab; sl.ha_flt_tab = ha_flt_tab;
sl.ha_idx = ha_idx; sl.ha_idx = ha_idx;
sl.opt = opt; sl.opt = opt;
sl.chunk_size = 500000000; sl.chunk_size = 100000000;
sl.n_thread = asm_opt.thread_num; sl.n_thread = asm_opt.thread_num;
sl.ug = ug; sl.ug = ug;
sl.rg = rg;
sl.uopt = uopt; sl.uopt = uopt;
scall_ul_pipeline(&sl, fn); scall_ul_pipeline(&sl, fn);
print_ul_rs(&UL_INF); // print_ul_rs(&UL_INF);
debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100);
// if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) { // if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) {
// scall_ul_pipeline(&sl, fn); // scall_ul_pipeline(&sl, fn);
// write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name); // write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name);
@@ -2891,13 +2893,13 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enz
return 1; return 1;
} }
ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, uint64_t n_read, int64_t readLen) ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz)
{ {
uint64_t i, k, qn, tn; uint64_t i, k, qn, tn, n_read = R_INF.total_reads;
int32_t r; asg_arc_t t, *p = NULL; int32_t r; asg_arc_t t, *p = NULL;
asg_t *rg = asg_init(); asg_t *rg = asg_init();
rg->m_seq = n_read; MALLOC(rg->seq, rg->m_seq); rg->m_seq = rg->n_seq = n_read; MALLOC(rg->seq, rg->m_seq);
for (i = 0; i < n_read; ++i) rg->seq[i].len = Get_READ_LENGTH(R_INF, i); for (i = 0; i < n_read; ++i) rg->seq[i].len = Get_READ_LENGTH(R_INF, i);
for (i = 0; i < n_read; i++) { for (i = 0; i < n_read; i++) {
@@ -2906,6 +2908,8 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui
if(!src[i].buffer[k].el) continue; if(!src[i].buffer[k].el) continue;
qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]);
if(rg->seq[qn].del || rg->seq[tn].del) continue; if(rg->seq[qn].del || rg->seq[tn].del) continue;
if((Get_qe(src[i].buffer[k]) - Get_qs(src[i].buffer[k])) < min_ovlp) continue;
if((Get_te(src[i].buffer[k]) - Get_ts(src[i].buffer[k])) < min_ovlp) continue;
r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if (r == MA_HT_QCONT) { if (r == MA_HT_QCONT) {
rg->seq[qn].del = 1; rg->seq[qn].del = 1;
@@ -2922,6 +2926,8 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui
if(!src[i].buffer[k].el) continue; if(!src[i].buffer[k].el) continue;
qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]);
if(rg->seq[qn].del || rg->seq[tn].del) continue; if(rg->seq[qn].del || rg->seq[tn].del) continue;
if((Get_qe(src[i].buffer[k]) - Get_qs(src[i].buffer[k])) < min_ovlp) continue;
if((Get_te(src[i].buffer[k]) - Get_ts(src[i].buffer[k])) < min_ovlp) continue;
r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if (r >= 0) { if (r >= 0) {
p = asg_arc_pushp(rg); p = asg_arc_pushp(rg);
@@ -2931,34 +2937,29 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, ui
} }
asg_cleanup(rg); asg_symm(rg); asg_cleanup(rg); asg_symm(rg);
asg_arc_del_trans(rg, gap_fuzz);
ma_ug_t *ug = NULL; ma_ug_t *ug = NULL;
ug = ma_ug_gen(rg); ug = ma_ug_gen(rg);
asg_destroy(rg);
for (i = k = 0; i < ug->u.n; i++) k += ug->u.a[i].len;
fprintf(stderr, "[M::%s::] # unitigs: %lu, # bases: %lu\n", __func__, (uint64_t)ug->u.n, k);
return ug;
asg_destroy(rg); ma_ug_destroy(ug);
} }
void ul_load(const ug_opt_t *uopt) void ul_load(const ug_opt_t *uopt)
{ {
fprintf(stderr, "[M::%s::] ==> UL\n", __func__); fprintf(stderr, "[M::%s::] ==> UL\n", __func__);
mg_idxopt_t opt; mg_idxopt_t opt;
ma_ug_t *ug = NULL; ma_ug_t *ug = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz);
int cutoff = asm_opt.hom_cov * asm_opt.high_factor; int cutoff = asm_opt.hom_cov * asm_opt.high_factor;
init_aux_table(); init_aux_table();
init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff); init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff);
/**
int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0);
if(exist == 0) uidx_l_build(ug, &opt, cutoff); if(exist == 0) uidx_l_build(ug, &opt, cutoff);
if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name);
**/
ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, ug);
ma_ug_destroy(ug); destory_all_ul_t(&UL_INF);
ul_v_call(&opt, uopt, /**rg**/NULL, asm_opt.ar, ha_flt_tab, ha_idx, /**ug**/NULL);
} }