fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-28 10:38:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1149: implemented random open syncmer

Browse Source 
On the mm2-update dataset, -j8 leads to sparser k-mer selection at higher
accuracy. The speed becomes a little slower. There seems a benefit but not a
big one.
This commit is contained in:
Heng Li
2022-10-21 19:07:28 -04:00
parent cefd0d9f6c
commit c2f07ff2ac
9 changed files with 84 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files
sketch.c
+55
View File
@@ -141,3 +141,58 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
if (min.x != UINT64_MAX)
kv_push(mm128_t, km, *p, min);
}
void mm_sketch_syncmer(void *km, const char *str, int len, int smer, int k, uint32_t rid, int is_hpc, mm128_v *p)
{
uint64_t shift1 = 2 * (k - 1), mask = (1ULL<<2*k) - 1, smask = (1ULL<<2*smer) - 1, kmer[2] = {0,0};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
tiny_queue_t tq;
assert(len > 0 && (smer > 0 && smer <= k) && (k > 0 && k <= 28)); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice
memset(&tq, 0, sizeof(tiny_queue_t));
kv_resize(mm128_t, km, *p, p->n + len/(k - smer));
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
kmer[0] = (kmer[0] << 2 | c) & mask; // forward k-mer
kmer[1] = (kmer[1] >> 2) | (3ULL^c) << shift1; // reverse k-mer
if (kmer[0] == kmer[1]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[0] < kmer[1]? 0 : 1; // strand
++l;
if (l >= k && kmer_span < 256) {
uint64_t x, min = UINT64_MAX;
x = hash64(kmer[z], mask);
for (j = 0; j < k - smer; ++j) {
uint64_t y = x >> (j + j) & smask;
min = min < y? min : y;
}
if ((x & smask) == min) {
mm128_t t;
t.x = x << 8 | kmer_span;
t.y = (uint64_t)rid<<32 | (uint32_t)i<<1 | z;
kv_push(mm128_t, km, *p, t);
}
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
}
}
void mm_sketch2(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, int is_syncmer, mm128_v *p)
{
if (is_syncmer) mm_sketch_syncmer(km, str, len, w, k, rid, is_hpc, p);
else mm_sketch(km, str, len, w, k, rid, is_hpc, p);
}