8000 minibwa/l2bit.c at dist · nh13/minibwa · GitHub
[go: up one dir, main page]
Skip to content

l2bit.c

Latest commit

 

History

History
810 lines (737 loc) · 25.2 KB

l2bit.c

File metadata and controls

810 lines (737 loc) · 25.2 KB
Raw
#define _POSIX_C_SOURCE 200809L
#include <zlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "kommon.h"
#include "l2bit.h"
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)
static int64_t l2b_pos2cid(const l2b_t *l2b, int64_t s, int64_t len, int64_t *cst)
{
int64_t lo = 0, hi = l2b->n_ctg, mid;
while (lo < hi) {
const l2b_ctg_t *ctg;
mid = (lo + hi) / 2;
ctg = &l2b->ctg[mid];
if (ctg->off <= s && s < ctg->off + ctg->len) {
*cst = s - ctg->off;
return s + len <= ctg->off + ctg->len? mid : -1;
} else if (s < ctg->off) hi = mid;
else lo = mid + 1;
}
return -1; // s is in no contig
}
int64_t l2b_intv2cid(const l2b_t *l2b, uint64_t st, uint64_t en, int64_t *cst, int *rev)
{
int64_t s;
assert(st < en);
if (en > l2b->tot_len * 2) return -3;
if (st < l2b->tot_len && l2b->tot_len < en) return -2;
*rev = (st >= l2b->tot_len);
s = st < l2b->tot_len? st : l2b->tot_len * 2 - en;
return l2b_pos2cid(l2b, s, en - st, cst);
}
int64_t l2b_intv2cid_meth(const l2b_t *l2b, uint64_t st, uint64_t en, l2b_meth_t *mt, int64_t *cst, int *rev)
{
int64_t s, len = en - st;
int32_t copy;
uint64_t tot_len = l2b->tot_len;
assert(st < en);
if (en > tot_len * 4) return -3;
copy = st / tot_len; // 0: c2t_f, 1: g2a_f, 2: g2a_r, 3: c2t_r
*mt = (copy == 0 || copy == 3)? L2B_METH_C2T : L2B_METH_G2A;
*rev = (copy >= 2);
s = st - tot_len * copy;
if (copy >= 2) s = tot_len - len - s; // flip for reverse copies
return s < 0? -2 : l2b_pos2cid(l2b, s, len, cst);
}
int64_t l2b_getseq(const l2b_t *l2b, int64_t tid, int64_t st, int64_t en, uint8_t *seq)
{
const l2b_ctg_t *ctg;
int64_t i, aid;
int32_t n_ambi;
if (tid < 0 || tid >= l2b->n_ctg) return -1;
ctg = &l2b->ctg[tid];
if (st < 0) st = 0;
if (en > ctg->len) en = ctg->len;
st += ctg->off;
en += ctg->off;
for (i = st; i < en; ++i)
seq[i - st] = l2b_get0(l2b, i);
// retrieve ambiguous bases
aid = l2b_getambi(l2b, tid, st - ctg->off, en - ctg->off, &n_ambi);
if (aid >= 0) {
for (i = 0; i < n_ambi; ++i) {
const l2b_intv_t *iv = &l2b->ambi[aid + i];
int64_t s = iv->st, e = iv->en;
if (s < st) s = st;
if (e > en) e = en;
if (s < e) memset(&seq[s - st], 4, e - s);
}
}
return en - st;
}
void l2b_meth_convert(l2b_meth_t mt, int64_t len, uint8_t *seq)
{
int64_t i;
if (mt == L2B_METH_C2T) {
for (i = 0; i < len; ++i)
if (seq[i] == 1) seq[i] = 3;
} else if (mt == L2B_METH_G2A) {
for (i = 0; i < len; ++i)
if (seq[i] == 2) seq[i] = 0;
}
}
int64_t l2b_getambi(const l2b_t *l2b, int64_t tid, int64_t st, int64_t en, int32_t *n_ambi)
{
int64_t g_beg, g_end, lo, hi, mid, i_st, i_en;
*n_ambi = 0;
if (tid < 0 || tid >= l2b->n_ctg) return -1;
if (st < 0) st = 0;
if (en > l2b->ctg[tid].len) en = l2b->ctg[tid].len;
if (st >= en) return -1;
g_beg = l2b->ctg[tid].off + st;
g_end = l2b->ctg[tid].off + en;
lo = 0, hi = l2b->n_ambi;
while (lo < hi) {
mid = (lo + hi) / 2;
if (l2b->ambi[mid].en > g_beg) hi = mid;
else lo = mid + 1;
}
i_st = lo;
lo = i_st, hi = l2b->n_ambi;
while (lo < hi) {
mid = (lo + hi) / 2;
if (l2b->ambi[mid].st >= g_end) hi = mid;
else lo = mid + 1;
}
i_en = lo;
*n_ambi = i_en - i_st;
if (*n_ambi == 0) return -1;
return i_st;
}
static void l2b_format_seq(uint64_t len, char *seq, uint64_t *rng)
{
uint64_t i;
for (i = 0; i < len; ++i) {
int b = (uint8_t)seq[i];
uint8_t c = kom_nt4_table[b];
if (c > 4) c = 4;
if (c == 4) c |= kom_splitmix64(rng) & 3;
if (b < 'A' || b > 'Z') c |= 1<<3;
seq[i] = c;
}
}
static void l2b_add_seq(l2b_t *l2b, uint64_t len, const char *seq, const char *name, const char *comm, uint64_t *rng)
{
uint64_t i, ambi_len, mask_len, off, m_pac_old;
l2b_ctg_t *ctg;
off = l2b->tot_len;
kom_grow(l2b_ctg_t, l2b->ctg, l2b->n_ctg, l2b->m_ctg);
ctg = &l2b->ctg[l2b->n_ctg++];
ctg->name = kom_strdup(name);
ctg->comm = comm? kom_strdup(comm) : 0;
ctg->len = len;
ctg->off = l2b->tot_len;
/* kom_grow() uses realloc(), which does NOT zero new slots, so the ALT fields
* must be initialized explicitly here. Without this, l2b_destroy() frees an
* uninitialized ctg->lift on the plain `minibwa index` path (l2b_set_alt is
* never called there), which crashes on any multi-contig reference. */
ctg->is_alt = 0;
ctg->n_lift = 0;
ctg->lift = 0;
l2b->tot_len += len;
m_pac_old = l2b->m_pac;
l2b->n_pac = (l2b->tot_len + 31) / 32;
kom_grow(uint64_t, l2b->pac, l2b->n_pac, l2b->m_pac);
if (m_pac_old < l2b->m_pac) // zero out newly allocated part
memset(&l2b->pac[m_pac_old], 0, (l2b->m_pac - m_pac_old) * 8);
for (i = 0, ambi_len = mask_len = 0; i < len; ++i) {
uint64_t c = (uint8_t)seq[i], x = off + i;
if (c & 1<<3) { // soft-masked base
++mask_len;
} else if (mask_len > 0) {
kom_grow(l2b_intv_t, l2b->mask, l2b->n_mask, l2b->m_mask);
l2b->mask[l2b->n_mask].st = x - mask_len;
l2b->mask[l2b->n_mask].en = x;
l2b->n_mask++;
mask_len = 0;
}
if (c & 1<<2) { // ambiguous base
++ambi_len;
} else if (ambi_len > 0) {
kom_grow(l2b_intv_t, l2b->ambi, l2b->n_ambi, l2b->m_ambi);
l2b->ambi[l2b->n_ambi].st = x - ambi_len;
l2b->ambi[l2b->n_ambi].en = x;
l2b->n_ambi++;
ambi_len = 0;
}
l2b->pac[x>>5] |= (c&3) << (x&0x1f)*2;
}
}
static void l2b_collate_str(l2b_t *l2b)
{
uint64_t i, tot_name = 0, tot_comm = 0;
char *p_name, *p_comm;
if (l2b->cat_name || l2b->cat_comm) return;
for (i = 0; i < l2b->n_ctg; ++i) {
l2b_ctg_t *ctg = &l2b->ctg[i];
tot_name += strlen(ctg->name) + 1;
tot_comm += ctg->comm? strlen(ctg->comm) + 1 : 1;
}
p_name = l2b->cat_name = kom_calloc(char, tot_name);
p_comm = l2b->cat_comm = kom_calloc(char, tot_comm);
for (i = 0; i < l2b->n_ctg; ++i) {
l2b_ctg_t *ctg = &l2b->ctg[i];
uint64_t len;
len = strlen(ctg->name);
memcpy(p_name, ctg->name, len + 1);
free(ctg->name);
ctg->name = p_name;
p_name += len + 1;
if (ctg->comm) {
len = strlen(ctg->comm);
memcpy(p_comm, ctg->comm, len + 1);
free(ctg->comm);
ctg->comm = p_comm;
} else len = 0;
p_comm += len + 1;
}
}
l2b_t *l2b_import(const char *fn, uint64_t seed)
{
gzFile fp;
kseq_t *ks;
l2b_t *l2b;
uint64_t rng = seed;
fp = fn == 0 || strcmp(fn, "-") == 0? gzdopen(0, "r") : gzopen(fn, "r");
if (fp == 0) return 0;
ks = kseq_init(fp);
l2b = kom_calloc(l2b_t, 1);
while (kseq_read(ks) >= 0) {
l2b_format_seq(ks->seq.l, ks->seq.s, &rng);
l2b_add_seq(l2b, ks->seq.l, ks->seq.s, ks->name.s, ks->comment.l? ks->comment.s : 0, &rng);
}
kseq_destroy(ks);
gzclose(fp);
l2b_collate_str(l2b);
return l2b;
}
#ifdef MB_HAVE_B2
/* bntseq.h is bwa-mem3's reference-sequence API; found via the
* -I$(BWAMEM3_DIR)/src path this translation unit's compile rule adds under
* `make B2=1` (Makefile). Guarded by MB_HAVE_B2 so a native (no-B2) build
* never needs that path and l2bit.o's command line is unchanged. */
#include "bntseq.h"
/* Phase-2: reconstruct an l2b_t from bwa-mem3's bns+pac instead of reading
* <prefix>.l2b. `bns_` is the opaque bntseq_t* from mb_b2_load_bns (b2idx.h);
* `pac`/`l_pac` are its packed 2-bit reference buffer and base count.
*
* Converts bwa's high-bit-first packed reference (_get_pac) to minibwa's
* low-bit-first uint64 word layout (l2b_get0/l2b_pac's convention: base i
* lives at bit offset (i&31)*2 within pac[i>>5]), and bwa's N-run holes
* (bns->ambs[], offset+len) into l2b's ambi[] (st/en), which is how
* l2b_getseq marks N bases (the base value packed at those pac positions is
* irrelevant -- l2b_getseq overwrites them with 4 via ambi -- matching the
* benign N-substitution difference already characterized between the two
* indexes' pac content).
*
* bns has no soft-mask concept, so l->mask/n_mask stay empty (kom_calloc
* zeroes the whole struct); nothing downstream of seeding+chaining+extension
* consults l2b's mask table. Deep-copies everything, so the caller may free
* bns/pac immediately after this returns. */
l2b_t *l2b_from_bns(const void *bns_, const uint8_t *pac, int64_t l_pac)
{
const bntseq_t *bns = (const bntseq_t*)bns_;
l2b_t *l = kom_calloc(l2b_t, 1);
int64_t i;
int t;
l->tot_len = l_pac;
l->n_ctg = l->m_ctg = bns->n_seqs;
l->ctg = kom_calloc(l2b_ctg_t, l->n_ctg);
for (t = 0; t < bns->n_seqs; ++t) {
const bntann1_t *ann = &bns->anns[t];
l->ctg[t].name = kom_strdup(ann->name);
l->ctg[t].comm = (ann->anno && ann->anno[0])? kom_strdup(ann->anno) : 0;
l->ctg[t].len = ann->len;
l->ctg[t].off = ann->offset;
}
l->n_pac = l->m_pac = (l_pac + 31) >> 5; // 32 bases per uint64
l->pac = kom_calloc(uint64_t, l->n_pac);
for (i = 0; i < l_pac; ++i) { // bit-order conversion: high-bit-first -> low-bit-first
int c = _get_pac(pac, i);
l->pac[i>>5] |= (uint64_t)c << ((i & 31) << 1);
}
l->n_ambi = l->m_ambi = bns->n_holes;
l->ambi = kom_calloc(l2b_intv_t, l->n_ambi? l->n_ambi : 1);
for (t = 0; t < bns->n_holes; ++t) {
l->ambi[t].st = bns->ambs[t].offset;
l->ambi[t].en = bns->ambs[t].offset + bns->ambs[t].len;
}
l->n_mask = l->m_mask = 0; // bns has no soft-mask
l2b_collate_str(l); // collate ctg[].name/comm into cat_name/cat_comm (matches l2b_import)
return l;
}
#endif
void l2b_destroy(l2b_t *l2b)
{
uint64_t i;
// lift[] is built by l2b_set_alt() with realloc() after the index is loaded,
// so it is heap-allocated in the mmap case too (ctg[] itself always is).
for (i = 0; i < l2b->n_ctg; ++i)
free(l2b->ctg[i].lift);
if (l2b->mmap) { // ambi/mask/pac/cat_name/cat_comm point into the mapped file
free(l2b->ctg); // ctg[] is always heap-allocated
kom_munmap(l2b->mmap, l2b->mmap_len);
} else {
free(l2b->cat_name); free(l2b->cat_comm);
free(l2b->pac); free(l2b->ambi); free(l2b->mask); free(l2b->ctg);
}
free(l2b);
}
/****************************
* ALT liftover index *
****************************/
/* Order lift blocks by ALT-forward start (then end) so l2b_lift()'s binary
* search over [alt_st, alt_en) is valid. Blocks are appended in file order,
* which is non-monotonic for reverse records (emitted descending after the
* RC->forward remap) and for ALT contigs carrying multiple .alt records. */
static int l2b_lift_cmp(const void *a, const void *b)
{
const l2b_lift_t *x = (const l2b_lift_t*)a, *y = (const l2b_lift_t*)b;
if (x->alt_st != y->alt_st) return x->alt_st < y->alt_st? -1 : 1;
if (x->alt_en != y->alt_en) return x->alt_en < y->alt_en? -1 : 1;
/* qsort is not stable, so equal keys would reorder arbitrarily between libc
* implementations and make the lift -- and therefore the SAM -- differ across
* platforms. Two records CAN tie on the ALT span and differ in target, which
* multi-record ALT contigs make reachable. Break every tie. */
if (x->pri_tid != y->pri_tid) return x->pri_tid < y->pri_tid? -1 : 1;
if (x->pri_st != y->pri_st) return x->pri_st < y->pri_st? -1 : 1;
if (x->rev != y->rev) return x->rev < y->rev? -1 : 1;
return 0;
}
int l2b_set_alt(l2b_t *l2b, const char *fn)
{
FILE *fp;
char *line = 0;
size_t line_cap = 0;
ssize_t line_len;
int n_alt = 0;
uint64_t i;
/* Reset: free any existing lift blocks and clear flags. */
for (i = 0; i < l2b->n_ctg; ++i) {
free(l2b->ctg[i].lift);
l2b->ctg[i].lift = 0;
l2b->ctg[i].n_lift = 0;
l2b->ctg[i].is_alt = 0;
}
l2b->n_alt_ctg = 0;
fp = fopen(fn, "r");
if (fp == 0) return -1;
while ((line_len = getline(&line, &line_cap, fp)) > 0) {
char *p, *q;
char *fields[12];
int nf;
int64_t alt_tid, pri_tid;
long pos1;
uint64_t pri_pos, alt_cursor, pri_cursor;
uint32_t flag;
uint8_t rev;
uint32_t m_lift, n_lift0;
int bad;
l2b_ctg_t *ctg;
/* Skip SAM header lines. */
if (line[0] == '@') continue;
/* Trim trailing newline. */
if (line_len > 0 && line[line_len-1] == '\n') line[--line_len] = '\0';
if (line_len > 0 && line[line_len-1] == '\r') line[--line_len] = '\0';
/* Split into at least 12 tab-delimited fields. */
for (nf = 0, p = line; nf < 12; ++nf) {
fields[nf] = p;
q = strchr(p, '\t');
if (q) { *q = '\0'; p = q + 1; }
else { ++nf; break; }
}
if (nf < 6) continue; /* need QNAME, FLAG, RNAME, POS, MAPQ, CIGAR */
/* Locate ALT contig by QNAME (col 0). */
alt_tid = -1;
for (i = 0; i < l2b->n_ctg; ++i)
if (strcmp(l2b->ctg[i].name, fields[0]) == 0) { alt_tid = i; break; }
if (alt_tid < 0) {
if (kom_verbose >= 2)
fprintf(stderr, "[W::%s] QNAME '%s' not in index, skipping\n", __func__, fields[0]);
continue;
}
/* Unmapped ALT/decoy entries (FLAG 0x4 or RNAME '*') carry no primary
* span; a standard hs38DH .alt ships ~2400 such decoy lines. They were
* never loadable, so skip them silently — warning on each is just noise
* and would flood stderr on every real run. */
if ((atol(fields[1]) & 0x4) || strcmp(fields[2], "*") == 0)
continue;
/* Locate primary contig by RNAME (col 2). */
pri_tid = -1;
for (i = 0; i < l2b->n_ctg; ++i)
if (strcmp(l2b->ctg[i].name, fields[2]) == 0) { pri_tid = i; break; }
if (pri_tid < 0) {
if (kom_verbose >= 2)
fprintf(stderr, "[W::%s] RNAME '%s' not in index, skipping\n", __func__, fields[2]);
continue;
}
flag = (uint32_t)atol(fields[1]);
rev = (flag & 0x10) ? 1 : 0;
/* POS is 1-based; a valid mapped record has POS >= 1. A malformed POS=0
* (or negative) would underflow to a garbage primary coordinate, so skip. */
pos1 = atol(fields[3]);
if (pos1 < 1) {
if (kom_verbose >= 2)
fprintf(stderr, "[W::%s] QNAME '%s' has invalid POS '%s', skipping\n", __func__, fields[0], fields[3]);
continue;
}
pri_pos = (uint64_t)(pos1 - 1);
ctg = &l2b->ctg[alt_tid];
/* n_alt_ctg counts contigs (a contig may carry several records via
* supplementary lines), so only the first record for a contig advances it. */
if (!ctg->is_alt) ++l2b->n_alt_ctg;
ctg->is_alt = 1;
++n_alt;
/* Walk CIGAR to build lift blocks. alt_cursor tracks the position along
* the CIGAR's query axis: for a forward record this IS the ALT-forward
* coordinate; for a reverse record it is the reverse-complement coordinate,
* so each M block is remapped to ALT-forward below. A final sort over all
* of the contig's blocks restores the alt_st ordering l2b_lift() needs. */
alt_cursor = 0;
pri_cursor = pri_pos;
m_lift = ctg->n_lift;
n_lift0 = ctg->n_lift; /* rollback point if the CIGAR is malformed */
bad = 0;
p = fields[5];
while (*p) {
uint64_t len = 0;
int op;
if (*p < '0' || *p > '9') { bad = 1; break; } /* expected an op length */
while (*p >= '0' && *p <= '9') len = len * 10 + (*p++ - '0');
if (*p == '\0') { bad = 1; break; } /* length with no operator (ran off the field) */
op = *p++;
if (op == 'M' || op == '=' || op == 'X') {
/* Emit one lift block. */
l2b_lift_t blk;
blk.alt_st = alt_cursor;
blk.alt_en = alt_cursor + len;
blk.pri_tid = pri_tid;
blk.pri_st = pri_cursor;
blk.pri_en = pri_cursor + len;
blk.rev = rev;
/* Reverse record: the CIGAR walks the reverse-complement of the ALT
* contig, so [alt_st, alt_en) is in RC coordinates. Remap to
* ALT-forward ([len_alt-alt_en, len_alt-alt_st)) while leaving
* pri_st/pri_en unchanged; l2b_lift()'s reverse branch then maps the
* forward query position back with its decreasing correspondence. */
if (rev) {
/* The remap subtracts from ctg->len, so a record whose query span
* runs past the contig it names would wrap to ~UINT64_MAX and sort
* to the end of lift[], where the binary search reads it as a valid
* block -- silently wrong lifts, no crash. Reject the record. */
if (blk.alt_en > ctg->len) { bad = 1; break; }
uint64_t fst = ctg->len - blk.alt_en;
uint64_t fen = ctg->len - blk.alt_st;
blk.alt_st = fst;
blk.alt_en = fen;
}
kom_grow(l2b_lift_t, ctg->lift, ctg->n_lift, m_lift);
ctg->lift[ctg->n_lift++] = blk;
alt_cursor += len;
pri_cursor += len;
} else if (op == 'I' || op == 'S' || op == 'H') {
/* Query/ALT-consuming ops with no primary footprint. A hard clip
* removes bases from SEQ but they still occupy the ALT contig, so
* (like S and I) it advances the ALT axis -- essential for a
* supplementary .alt record whose leading H is its offset into the
* ALT contig. */
alt_cursor += len;
} else if (op == 'D' || op == 'N') {
/* Deletion from ALT: advance only pri_cursor. */
pri_cursor += len;
} else if (op == 'P') {
/* Padding: consumes neither coordinate. */
} else {
bad = 1; break; /* unknown CIGAR operator */
}
}
if (bad) {
/* Roll back this record's partial blocks. The buffer keeps whatever
* capacity the partial append grew it to; shrink so the next record's
* `m_lift = ctg->n_lift` still describes the real allocation, as it does
* on every non-rollback path. */
ctg->n_lift = n_lift0;
if (ctg->n_lift > 0 && ctg->n_lift < m_lift)
ctg->lift = kom_realloc(l2b_lift_t, ctg->lift, ctg->n_lift);
if (kom_verbose >= 2)
fprintf(stderr, "[W::%s] QNAME '%s' has an unparsable CIGAR '%s', skipping record\n", __func__, fields[0], fields[5]);
continue;
}
/* Shrink the lift array to exact size. */
if (ctg->n_lift > 0 && ctg->n_lift < m_lift)
ctg->lift = kom_realloc(l2b_lift_t, ctg->lift, ctg->n_lift);
}
free(line);
fclose(fp);
/* Blocks were appended in file order. Reverse records and multi-record ALT
* contigs make lift[] non-monotonic in alt_st, but l2b_lift()'s binary search
* requires it sorted. Sort each ALT contig's blocks now that every record is
* loaded. */
for (i = 0; i < l2b->n_ctg; ++i) {
l2b_ctg_t *c = &l2b->ctg[i];
if (c->n_lift > 1)
qsort(c->lift, c->n_lift, sizeof(l2b_lift_t), l2b_lift_cmp);
}
return n_alt;
}
int l2b_lift(const l2b_t *l2b, int64_t alt_tid, uint64_t alt_pos,
int64_t *pri_tid, uint64_t *pri_pos, uint8_t *rev)
{
const l2b_ctg_t *ctg;
int64_t lo, hi, mid;
if (alt_tid < 0 || alt_tid >= (int64_t)l2b->n_ctg) return 0;
ctg = &l2b->ctg[alt_tid];
if (!ctg->is_alt || ctg->n_lift == 0) return 0;
/* Binary search for the block whose [alt_st, alt_en) contains alt_pos. */
lo = 0; hi = ctg->n_lift;
while (lo < hi) {
mid = (lo + hi) / 2;
if (ctg->lift[mid].alt_en <= alt_pos) lo = mid + 1;
else hi = mid;
}
if (lo >= (int64_t)ctg->n_lift) return 0;
{
const l2b_lift_t *blk = &ctg->lift[lo];
if (alt_pos < blk->alt_st) return 0; /* in a hole */
*pri_tid = blk->pri_tid;
*rev = blk->rev;
if (!blk->rev) {
/* Forward: primary coord increases with alt coord. */
*pri_pos = blk->pri_st + (alt_pos - blk->alt_st);
} else {
/* Reverse: primary coord decreases as alt coord increases. */
*pri_pos = blk->pri_en - 1 - (alt_pos - blk->alt_st);
}
}
return 1;
}
int l2b_save(const char *fn, const l2b_t *l2b)
{
FILE *fp;
uint64_t i, len_name = 0, len_comm = 0;
uint32_t dummy = 0;
fp = fn == 0 || strcmp(fn, "-") == 0? stdout : fopen(fn, "wb");
if (fp == 0) return -1;
for (i = 0; i < l2b->n_ctg; ++i) {
const l2b_ctg_t *ctg = &l2b->ctg[i];
len_name += strlen(ctg->name) + 1;
len_comm += ctg->comm? strlen(ctg->comm) + 1 : 1;
}
fwrite(L2B_MAGIC, 1, 4, fp);
fwrite(&dummy, 4, 1, fp);
fwrite(&l2b->n_ctg, 8, 1, fp);
fwrite(&l2b->tot_len, 8, 1, fp);
fwrite(&l2b->n_ambi, 8, 1, fp);
fwrite(&l2b->n_mask, 8, 1, fp);
fwrite(&len_name, 8, 1, fp);
fwrite(&len_comm, 8, 1, fp);
fwrite(&l2b->n_pac, 8, 1, fp);
for (i = 0; i < l2b->n_ctg; ++i)
fwrite(&l2b->ctg[i].len, 8, 1, fp);
fwrite(l2b->ambi, 16, l2b->n_ambi, fp);
fwrite(l2b->mask, 16, l2b->n_mask, fp);
fwrite(l2b->pac, 8, l2b->n_pac, fp);
fwrite(l2b->cat_name, 1, len_name, fp); // put strings at the end to make sure uint64_t are all aligned
fwrite(l2b->cat_comm, 1, len_comm, fp);
fclose(fp);
return 0;
}
l2b_t *l2b_load(const char *fn)
{
FILE *fp;
char magic[4], *p_name, *p_comm;
uint32_t dummy;
uint64_t off, i, len_name, len_comm;
l2b_t *l2b;
fp = fn == 0 || strcmp(fn, "-") == 0? stdin : fopen(fn, "rb");
if (fp == 0) return 0;
fread(magic, 1, 4, fp);
if (strncmp(magic, L2B_MAGIC, 4) != 0) {
if (fp != stdin) fclose(fp);
return 0;
}
l2b = kom_calloc(l2b_t, 1);
fread(&dummy, 4, 1, fp);
fread(&l2b->n_ctg, 8, 1, fp);
fread(&l2b->tot_len, 8, 1, fp);
fread(&l2b->n_ambi, 8, 1, fp);
fread(&l2b->n_mask, 8, 1, fp);
fread(&len_name, 8, 1, fp);
fread(&len_comm, 8, 1, fp);
fread(&l2b->n_pac, 8, 1, fp);
l2b->ctg = kom_calloc(l2b_ctg_t, l2b->n_ctg);
for (i = 0, off = 0; i < l2b->n_ctg; ++i) { // read contig lengths
fread(&l2b->ctg[i].len, 8, 1, fp);
l2b->ctg[i].off = off;
off += l2b->ctg[i].len;
}
if (off != l2b->tot_len) goto load_failure;
l2b->ambi = kom_malloc(l2b_intv_t, l2b->n_ambi);
l2b->mask = kom_malloc(l2b_intv_t, l2b->n_mask);
l2b->pac = kom_malloc(uint64_t, l2b->n_pac);
l2b->cat_name = kom_malloc(char, len_name);
l2b->cat_comm = kom_malloc(char, len_comm);
fread(l2b->ambi, 16, l2b->n_ambi, fp);
fread(l2b->mask, 16, l2b->n_mask, fp);
fread(l2b->pac, 8, l2b->n_pac, fp);
fread(l2b->cat_name, 1, len_name, fp);
fread(l2b->cat_comm, 1, len_comm, fp);
p_name = l2b->cat_name, p_comm = l2b->cat_comm;
for (i = 0; i < l2b->n_ctg; ++i) { // synchronize contig names and comments
l2b_ctg_t *ctg = &l2b->ctg[i];
ctg->name = p_name;
p_name += strlen(p_name) + 1;
ctg->comm = *p_comm? p_comm : 0;
p_comm += *p_comm? strlen(p_comm) + 1 : 1;
}
if (p_name - l2b->cat_name != len_name || p_comm - l2b->cat_comm != len_comm) goto load_failure;
if (fp != stdin) fclose(fp);
return l2b;
load_failure:
if (fp != stdin) fclose(fp);
l2b_destroy(l2b);
return 0;
}
l2b_t *l2b_load_mmap(const char *fn, int preload)
{
uint8_t *base;
const uint64_t *hdr, *lens;
uint64_t off, i, len_name, len_comm, foff;
size_t map_len;
char *p_name, *p_comm;
l2b_t *l2b;
base = (uint8_t*)kom_mmap_file(fn, &map_len, preload);
if (base == 0) return 0;
if (map_len < 64 || strncmp((const char*)base, L2B_MAGIC, 4) != 0) { kom_munmap(base, map_len); return 0; }
l2b = kom_calloc(l2b_t, 1);
l2b->mmap = base;
l2b->mmap_len = map_len;
hdr = (const uint64_t*)base; // hdr[0] is magic+dummy; fields start at hdr[1]
l2b->n_ctg = hdr[1];
l2b->tot_len = hdr[2];
l2b->n_ambi = hdr[3];
l2b->n_mask = hdr[4];
len_name = hdr[5];
len_comm = hdr[6];
l2b->n_pac = hdr[7];
l2b->ctg = kom_calloc(l2b_ctg_t, l2b->n_ctg); // ctg[] is not stored in the file
lens = &hdr[8]; // contig lengths follow the 64-byte header
for (i = 0, off = 0; i < l2b->n_ctg; ++i) {
l2b->ctg[i].len = lens[i];
l2b->ctg[i].off = off;
off += lens[i];
}
if (off != l2b->tot_len) goto mmap_failure;
foff = 64 + l2b->n_ctg * 8; // point large arrays into the mapped file
l2b->ambi = (l2b_intv_t*)(base + foff); foff += l2b->n_ambi * 16;
l2b->mask = (l2b_intv_t*)(base + foff); foff += l2b->n_mask * 16;
l2b->pac = (uint64_t*) (base + foff); foff += l2b->n_pac * 8;
l2b->cat_name = (char*)(base + foff); foff += len_name;
l2b->cat_comm = (char*)(base + foff); foff += len_comm;
if (foff != map_len) goto mmap_failure;
p_name = l2b->cat_name, p_comm = l2b->cat_comm;
for (i = 0; i < l2b->n_ctg; ++i) { // synchronize contig names and comments
l2b_ctg_t *ctg = &l2b->ctg[i];
ctg->name = p_name;
p_name += strlen(p_name) + 1;
ctg->comm = *p_comm? p_comm : 0;
p_comm += *p_comm? strlen(p_comm) + 1 : 1;
}
if (p_name - l2b->cat_name != len_name || p_comm - l2b->cat_comm != len_comm) goto mmap_failure;
return l2b;
mmap_failure:
l2b_destroy(l2b);
return 0;
}
/******************************
* Save .paf files for bwtgen *
******************************/
int l2b_save_pac(const char *fn, const l2b_t *l2b, int both_strand)
{
FILE *fp;
uint64_t n_pac, x;
int64_t i;
uint8_t *pac, ct;
fp = fn == 0 || strcmp(fn, "-") == 0? stdout : fopen(fn, "wb");
if (fp == 0) return -1;
// fill pac[]
n_pac = ((both_strand? l2b->tot_len * 2 : l2b->tot_len) + 3) / 4;
pac = kom_calloc(uint8_t, n_pac);
for (i = 0, x = 0; i < l2b->tot_len; ++i, ++x)
pac[x>>2] |= l2b_get0(l2b, i) << (~x&3) * 2;
if (both_strand)
for (i = l2b->tot_len - 1; i >= 0; --i, ++x)
pac[x>>2] |= (3 - l2b_get0(l2b, i)) << (~x&3) * 2;
// write pac
fwrite(pac, 1, (x>>2) + ((x&3) == 0? 0 : 1), fp);
// the following codes make the pac file size always (x/4+1+1)
if (x % 4 == 0) {
ct = 0;
fwrite(&ct, 1, 1, fp);
}
ct = x % 4;
fwrite(&ct, 1, 1, fp);
fclose(fp);
free(pac);
return 0;
}
static inline uint8_t l2b_c2t(uint8_t b) { return b == 1? 3 : b; } // C(1) -> T(3)
static inline uint8_t l2b_g2a(uint8_t b) { return b == 2? 0 : b; } // G(2) -> A(0)
int l2b_save_pac_meth(const char *fn, const l2b_t *l2b, int both_strand)
{
FILE *fp;
uint64_t n_pac, x, len;
int64_t i;
uint8_t *pac, ct, b;
fp = fn == 0 || strcmp(fn, "-") == 0? stdout : fopen(fn, "wb");
if (fp == 0) return -1;
len = l2b->tot_len * 2; // c2t + g2a
if (both_strand) len *= 2; // forward + reverse
n_pac = (len + 3) / 4;
pac = kom_calloc(uint8_t, n_pac);
// c2t forward
for (i = 0, x = 0; i < l2b->tot_len; ++i, ++x) {
b = l2b_c2t(l2b_get0(l2b, i));
pac[x>>2] |= b << ((~x&3) * 2);
}
// g2a forward
for (i = 0; i < l2b->tot_len; ++i, ++x) {
b = l2b_g2a(l2b_get0(l2b, i));
pac[x>>2] |= b << ((~x&3) * 2);
}
if (both_strand) {
// g2a reverse (reverse complement of g2a converted sequence)
for (i = l2b->tot_len - 1; i >= 0; --i, ++x) {
b = 3 - l2b_g2a(l2b_get0(l2b, i));
pac[x>>2] |= b << ((~x&3) * 2);
}
// c2t reverse (reverse complement of c2t converted sequence)
for (i = l2b->tot_len - 1; i >= 0; --i, ++x) {
b = 3 - l2b_c2t(l2b_get0(l2b, i));
pac[x>>2] |= b << ((~x&3) * 2);
}
}
fwrite(pac, 1, (x>>2) + ((x&3) == 0? 0 : 1), fp);
if (x % 4 == 0) {
ct = 0;
fwrite(&ct, 1, 1, fp);
}
ct = x % 4;
fwrite(&ct, 1, 1, fp);
fclose(fp);
free(pac);
return 0;
}
0