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

r255: faster way to extract reads

Browse Source
This commit is contained in:
Heng Li
2020-06-03 21:38:59 -04:00
parent eddbc61952
commit 97f6efc4da
2 changed files with 35 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CommandLines.h
View File

@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.7-dirty-r254"
#define HA_VERSION "0.7-dirty-r255"
#define VERBOSE 0

60
extract.cpp
View File

@@ -5,7 +5,8 @@
#include "kseq.h"
typedef const char *cstr_t;
KHASHL_CMAP_INIT(KH_LOCAL, strmap_t, ss, cstr_t, int, kh_hash_str, kh_eq_str)
KHASHL_CSET_INIT(KH_LOCAL, strset_t, ss, cstr_t, kh_hash_str, kh_eq_str)
KHASHL_MAP_INIT(KH_LOCAL, hm64_t, h64, uint64_t, int, kh_hash_uint64, kh_eq_generic)
KSTREAM_INIT(gzFile, gzread, 65536)
#define GFA_MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
@@ -76,25 +77,36 @@ char **gv_read_list(const char *o, int *n_)
void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
{
strmap_t *h;
hm64_t *h;
khint_t k;
int i, absent, m, l, max_len = 0;
int i, absent, m, l;
uint64_t j;
char *s = 0;
const ma_hit_t_alloc *ov[2] = { rs->paf, rs->reverse_paf };
FILE *fp = stdout;
for (j = 0; j < rs->total_reads; ++j) {
if (n > 0) {
int max_len = 0;
char *s = 0;
strset_t *ss;
ss = ss_init();
for (i = 0; i < n; ++i)
ss_put(ss, list[i], &absent);
for (j = 0; j < rs->total_reads; ++j)
if (max_len < (int)Get_NAME_LENGTH(*rs, j))
max_len = Get_NAME_LENGTH(*rs, j);
}
GFA_MALLOC(s, max_len + 1);
h = ss_init();
for (i = 0; i < n; ++i) {
k = ss_put(h, gfa_strdup(list[i]), &absent);
h = h64_init();
for (j = 0; j < rs->total_reads; ++j) {
strncpy(s, Get_NAME(*rs, j), Get_NAME_LENGTH(*rs, j));
s[Get_NAME_LENGTH(*rs, j)] = 0;
if (ss_get(ss, s) != kh_end(ss)) {
k = h64_put(h, j, &absent);
kh_val(h, k) = -1;
}
}
free(s);
ss_destroy(ss);
} else return;
for (m = 0; m < n_rounds; ++m) {
for (j = 0; j < rs->total_reads; ++j) {
@@ -104,24 +116,18 @@ void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
uint64_t q = Get_qn(o->buffer[i]);
uint64_t t = Get_tn(o->buffer[i]);
int q_hit = 0, t_hit = 0;
strncpy(s, Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q)); s[Get_NAME_LENGTH(*rs, q)] = 0;
k = ss_get(h, s);
k = h64_get(h, q);
q_hit = (k < kh_end(h) && kh_val(h, k) < m);
strncpy(s, Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t)); s[Get_NAME_LENGTH(*rs, t)] = 0;
k = ss_get(h, s);
k = h64_get(h, t);
t_hit = (k < kh_end(h) && kh_val(h, k) < m);
if ((!q_hit && !t_hit) || (q_hit && t_hit)) continue;
if (!q_hit) {
char *tmp = gfa_strndup(Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q));
k = ss_put(h, tmp, &absent);
k = h64_put(h, q, &absent);
if (absent) kh_val(h, k) = m;
else free(tmp);
}
if (!t_hit) {
char *tmp = gfa_strndup(Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t));
k = ss_put(h, tmp, &absent);
k = h64_put(h, t, &absent);
if (absent) kh_val(h, k) = m;
else free(tmp);
}
}
}
@@ -135,12 +141,8 @@ void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
uint64_t q = Get_qn(o->buffer[i]);
uint64_t t = Get_tn(o->buffer[i]);
int q_hit = 0, t_hit = 0;
strncpy(s, Get_NAME(*rs, q), Get_NAME_LENGTH(*rs, q)); s[Get_NAME_LENGTH(*rs, q)] = 0;
k = ss_get(h, s);
q_hit = (k < kh_end(h));
strncpy(s, Get_NAME(*rs, t), Get_NAME_LENGTH(*rs, t)); s[Get_NAME_LENGTH(*rs, t)] = 0;
k = ss_get(h, s);
t_hit = (k < kh_end(h));
q_hit = (h64_get(h, q) < kh_end(h));
t_hit = (h64_get(h, t) < kh_end(h));
if (!q_hit && !t_hit) continue;
fwrite(Get_NAME(*rs, q), 1, Get_NAME_LENGTH(*rs, q), fp);
fwrite("\t", 1, 1, fp);
@@ -157,11 +159,7 @@ void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
}
}
for (k = 0; k != kh_end(h); ++k)
if (kh_exist(h, k))
free((char*)kh_key(h, k));
ss_destroy(h);
free(s);
h64_destroy(h);
}
void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o)