c_cpp 从fermi2指数获得k-mer计数

Posted 2021-05-14

/* To compile this program, you need rld0.{h,c} from ropebwt2 or fermi2:
 *
 *   gcc -g -O2 -Wall -o occ1 rld0.c occ1.c
 */

#include <string.h>
#include <stdint.h>
#include <stdio.h>
#include "rld0.h"

static unsigned char seq_nt6_table[128] = { // lookup table
    0, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 1, 5, 2,  5, 5, 5, 3,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  4, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 1, 5, 2,  5, 5, 5, 3,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  4, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5
};

int main(int argc, char *argv[])
{
	int i, j;
	rld_t *e;
	if (argc < 3) {
		fprintf(stderr, "Usage: occ1 <index.fmd> <string1> [string2 [...]]\n");
		return 1;
	}
	e = rld_restore_mmap(argv[1]);
	for (i = 2; i < argc; ++i) {
		char *s = argv[i];
		uint64_t k = 0, l = e->mcnt[0]; // e->mcnt[0]: the total number of symbols
		for (j = strlen(s) - 1; j >= 0; --j) {
			int c = s[j] >= 0 && s[j] <= 127? seq_nt6_table[(int)s[j]] : 5; // nt6 base encoding
			k = e->cnt[c] + rld_rank11(e, k, c); // using rld_rank21() will be faster
			l = e->cnt[c] + rld_rank11(e, l, c);
			if (k == l) break; // no occurrences
		}
		printf("%lld\t%s\n", (long long)(l - k), s);
	}
	rld_destroy(e);
	return 0;
}