sketch.c

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#define __STDC_LIMIT_MACROS
#include "kvec.h"
#include "umpriv.h"
#include "ksort.h"

unsigned char seq_nt4_table[256] = {
	0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  3, 3, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  3, 3, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

typedef struct { // a simplified version of kdq
	int front, count;
	int a[32];
} tiny_queue_t;

static inline void tq_push(tiny_queue_t *q, int x)
{
	q->a[((q->count++) + q->front) & 0x1f] = x;
}

static inline int tq_shift(tiny_queue_t *q)
{
	int x;
	if (q->count == 0) return -1;
	x = q->a[q->front++];
	q->front &= 0x1f;
	--q->count;
	return x;
}

static inline int mzcmp(const mm128_t *a, const mm128_t *b) // TODO: we only need <=
{
	int32_t ya = a->y>>32, yb = b->y>>32;
	return ya < yb? -1 : ya > yb? 1 : ((a->x > b->x) - (a->x < b->x));
}

#define MAX_MAX_HIGH_OCC 128

#define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0)
KSORT_INIT(mz, mm128_t, mz_lt)

void mm_select_mz(mm128_v *p, int n0, int len, int dist)
{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak
	int32_t i, last0, n = (int32_t)p->n, m = 0;
	mm128_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation

	if (n - n0 == 0 || n - n0 == 1) return;
	for (i = n0; i < n; ++i)
		if (p->a[i].y>>32 > 2) ++m;
	if (m == 0) return; // no high-frequency k-mers; do nothing
	for (i = n0, last0 = n0 - 1; i <= n; ++i) {
		if (i == n || p->a[i].y>>32 <= 2) {
			if (i - last0 > 1) {
				int32_t ps = last0 < n0? 0 : (uint32_t)p->a[last0].y>>1;
				int32_t pe = i == n? len : (uint32_t)p->a[i].y>>1;
				int32_t j, k, st = last0 + 1, en = i;
				int32_t max_high_occ = (int32_t)((double)(pe - ps) / dist + .499);
				if (max_high_occ > MAX_MAX_HIGH_OCC)
					max_high_occ = MAX_MAX_HIGH_OCC;
				for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
					b[k] = p->a[j], b[k].y = (b[k].y&0xFFFFFFFF00000000ULL) | j; // lower 32 bits of b[].y keep the index in p->a[]
				ks_heapmake_mz(k, b); // initialize the binomial heap
				for (; j < en; ++j) { // if there are more, choose top max_high_occ
					if (mz_lt(p->a[j], b[0])) { // then update the heap
						b[0] = p->a[j], b[0].y = (b[0].y&0xFFFFFFFF00000000ULL) | j;
						ks_heapdown_mz(0, k, b);
					}
				}
				//ks_heapsort_mz(k, b); // sorting is not needed for now
				for (j = 0; j < k; ++j)
					p->a[(uint32_t)b[j].y].y |= 1ULL<<63;
				for (j = st; j < en; ++j) p->a[j].y ^= 1ULL<<63;
			}
			last0 = i;
		}
	}
	for (i = n = n0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
		if ((p->a[i].y&1ULL<<63) == 0)
			p->a[n++] = p->a[i];
	p->n = n;
}

/**
 * Find symmetric (w,k)-minimizers on a DNA sequence
 *
 * @param km     thread-local memory pool; using NULL falls back to malloc()
 * @param str    DNA sequence
 * @param len    length of $str
 * @param w      find a minimizer for every $w consecutive k-mers
 * @param k      k-mer size
 * @param rid    reference ID; will be copied to the output $p array
 * @param is_hpc homopolymer-compressed or not
 * @param p      minimizers
 *               p->a[i].x = kMer<<8 | kmerSpan
 *               p->a[i].y = rid<<32 | lastPos<<1 | strand
 *               where lastPos is the position of the last base of the i-th minimizer,
 *               and strand indicates whether the minimizer comes from the top or the bottom strand.
 *               Callers may want to set "p->n = 0"; otherwise results are appended to p
 */
void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p, const void *di, int skip_bf, int adap_occ, int adap_dist)
{
	uint64_t shift1 = 2 * (k - 1), mask = (1ULL<<2*k) - 1, kmer[2] = {0,0};
	int i, j, l, buf_pos, min_pos, kmer_span = 0, n0 = p->n;
	mm128_t buf[256], min = { UINT64_MAX, UINT64_MAX };
	tiny_queue_t tq;

	assert(len > 0 && (w > 0 && w < 256) && (k > 0 && k <= 28)); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice
	memset(buf, 0xff, w * 16);
	memset(&tq, 0, sizeof(tiny_queue_t));
	kv_resize(mm128_t, km, *p, p->n + len/w);

	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
		int c = seq_nt4_table[(uint8_t)str[i]];
		mm128_t info = { UINT64_MAX, UINT64_MAX };
		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) {
				int32_t c = 0;
				uint64_t hash;
				hash = um_hash64(kmer[z], mask);
				if (di) {
					c = um_didx_get(di, hash, skip_bf);
					if (c > 2 && c <= adap_occ) c = 2;
				}
				if (c >= 0) {
					info.x = hash << 8 | kmer_span;
					info.y = (uint64_t)c << 32 | (uint32_t)i<<1 | z;
				}
			}
		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
			for (j = buf_pos + 1; j < w; ++j)
				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
			for (j = 0; j < buf_pos; ++j)
				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
		}
		if (info.x <= min.x) { // a new minimum; then write the old min
			if (l >= w + k && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min);
			min = info, min_pos = buf_pos;
		} else if (buf_pos == min_pos) { // old min has moved outside the window
			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min);
			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
				if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
			for (j = 0; j <= buf_pos; ++j)
				if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j;
			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
				for (j = 0; j <= buf_pos; ++j)
					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
			}
		}
		if (++buf_pos == w) buf_pos = 0;
	}
	if (min.x != UINT64_MAX)
		kv_push(mm128_t, km, *p, min);
	if (p->n - n0 >= 2) {
		int n_del = 0;
		for (i = n0 + 1, j = n0; i <= (int)p->n; ++i) {
			if (i == (int)p->n || p->a[i].x != p->a[j].x) {
				if (i - j >= 2) {
					int st, en;
					if (i == (int)p->n || (j > n0 && p->a[j-1].x < p->a[i].x))
						st = j + 1, en = i;
					else st = j, en = i - 1;
					for (l = st; l < en; ++l)
						p->a[l].x = (uint64_t)-1;
					n_del += i - j - 1;
				}
				j = i;
			}
		}
		if (n_del > 0) {
			for (i = n0, j = n0; i < (int)p->n; ++i)
				if (p->a[i].x != (uint64_t)-1)
					p->a[j++] = p->a[i];
			p->n = j;
		}
	}
//	for (i = n0; i < (int)p->n; ++i) fprintf(stderr, "X\t%d\t%d\n", (int32_t)p->a[i].y, (int32_t)(p->a[i].y>>32));
	if (adap_dist > w && di) mm_select_mz(p, n0, len, adap_dist);
//	for (i = n0; i < (int)p->n; ++i) fprintf(stderr, "Y\t%d\t%d\n", (int32_t)p->a[i].y, (int32_t)(p->a[i].y>>32));
	for (i = n0; i < (int)p->n; ++i)
		p->a[i].y = p->a[i].y << 32 >> 32 | (uint64_t)rid << 32;
}