zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Use wrapper functions to catch system errors 

Use the wrapper functions in utils.c plus a few extra bits of error
checking code to catch system errors and exit non-zero when they occur.
This commit is contained in:
Rob Davies 2012-12-16 10:34:57 +00:00
parent 752ce69b78
commit b081ac9b8b
32 changed files with 411 additions and 355 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
Makefile
View File

@ -31,19 +31,38 @@ bwa:libbwa.a $(AOBJS) main.o
libbwa.a:$(LOBJS) libbwa.a:$(LOBJS)
$(AR) -csru $@ $(LOBJS) $(AR) -csru $@ $(LOBJS)
bwa.o:bwa.h
QSufSort.o:QSufSort.h
bwt.o:bwt.h
bwtio.o:bwt.h
bwtaln.o:bwt.h bwtaln.h kseq.h
bntseq.o:bntseq.h
bwtgap.o:bwtgap.h bwtaln.h bwt.h
bwtsw2_core.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h
bwtsw2_aux.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h
bwtsw2_main.o:bwtsw2.h
clean: clean:
rm -f gmon.out *.o a.out $(PROG) *~ *.a rm -f gmon.out *.o a.out $(PROG) *~ *.a
QSufSort.o: QSufSort.h
bamlite.o: bamlite.h utils.h
bntseq.o: bntseq.h kseq.h main.h utils.h
bwa.o: bntseq.h bwa.h bwt.h bwtaln.h bwtgap.h stdaln.h utils.h
bwape.o: bntseq.h bwase.h bwt.h bwtaln.h khash.h ksort.h kvec.h stdaln.h
bwape.o: utils.h
bwase.o: bntseq.h bwase.h bwt.h bwtaln.h kstring.h stdaln.h utils.h
bwaseqio.o: bamlite.h bwt.h bwtaln.h kseq.h stdaln.h utils.h
bwt.o: bwt.h kvec.h utils.h
bwt_gen.o: QSufSort.h utils.h
bwt_lite.o: bwt_lite.h utils.h
bwtaln.o: bwt.h bwtaln.h bwtgap.h stdaln.h utils.h
bwtgap.o: bwt.h bwtaln.h bwtgap.h stdaln.h utils.h
bwtindex.o: bntseq.h bwt.h main.h utils.h
bwtio.o: bwt.h utils.h
bwtmisc.o: bntseq.h bwt.h main.h utils.h
bwtsw2_aux.o: bntseq.h bwt.h bwt_lite.h bwtsw2.h kseq.h ksort.h kstring.h
bwtsw2_aux.o: stdaln.h utils.h
bwtsw2_chain.o: bntseq.h bwt.h bwt_lite.h bwtsw2.h ksort.h utils.h
bwtsw2_core.o: bntseq.h bwt.h bwt_lite.h bwtsw2.h khash.h ksort.h kvec.h
bwtsw2_core.o: utils.h
bwtsw2_main.o: bntseq.h bwt.h bwt_lite.h bwtsw2.h utils.h
bwtsw2_pair.o: bntseq.h bwt.h bwt_lite.h bwtsw2.h kstring.h ksw.h utils.h
cs2nt.o: bwt.h bwtaln.h stdaln.h utils.h
fastmap.o: bntseq.h bwt.h kseq.h kvec.h utils.h
is.o: utils.h
kstring.o: kstring.h utils.h
ksw.o: ksw.h utils.h
main.o: main.h utils.h
simple_dp.o: kseq.h stdaln.h utils.h
stdaln.o: stdaln.h utils.h
utils.o: utils.h

38
bamlite.c
View File

@ -2,6 +2,7 @@
#include <ctype.h> #include <ctype.h>
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include "utils.h"
#include "bamlite.h" #include "bamlite.h"
/********************* /*********************
@ -53,7 +54,7 @@ int bam_is_be;
bam_header_t *bam_header_init() bam_header_t *bam_header_init()
{ {
bam_is_be = bam_is_big_endian(); bam_is_be = bam_is_big_endian();
return (bam_header_t*)calloc(1, sizeof(bam_header_t)); return (bam_header_t*)xcalloc(1, sizeof(bam_header_t));
} }
void bam_header_destroy(bam_header_t *header) void bam_header_destroy(bam_header_t *header)
@ -62,11 +63,11 @@ void bam_header_destroy(bam_header_t *header)
if (header == 0) return; if (header == 0) return;
if (header->target_name) { if (header->target_name) {
for (i = 0; i < header->n_targets; ++i) for (i = 0; i < header->n_targets; ++i)
free(header->target_name[i]); if (header->target_name[i]) free(header->target_name[i]);
if (header->target_len) free(header->target_len);
free(header->target_name); free(header->target_name);
free(header->target_len);
} }
free(header->text); if (header->text) free(header->text);
free(header); free(header);
} }
@ -80,28 +81,33 @@ bam_header_t *bam_header_read(bamFile fp)
magic_len = bam_read(fp, buf, 4); magic_len = bam_read(fp, buf, 4);
if (magic_len != 4 || strncmp(buf, "BAM\001", 4) != 0) { if (magic_len != 4 || strncmp(buf, "BAM\001", 4) != 0) {
fprintf(stderr, "[bam_header_read] invalid BAM binary header (this is not a BAM file).\n"); fprintf(stderr, "[bam_header_read] invalid BAM binary header (this is not a BAM file).\n");
return 0; return NULL;
} }
header = bam_header_init(); header = bam_header_init();
// read plain text and the number of reference sequences // read plain text and the number of reference sequences
bam_read(fp, &header->l_text, 4); if (bam_read(fp, &header->l_text, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&header->l_text); if (bam_is_be) bam_swap_endian_4p(&header->l_text);
header->text = (char*)calloc(header->l_text + 1, 1); header->text = (char*)xcalloc(header->l_text + 1, 1);
bam_read(fp, header->text, header->l_text); if (bam_read(fp, header->text, header->l_text) != header->l_text) goto fail;
bam_read(fp, &header->n_targets, 4); if (bam_read(fp, &header->n_targets, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&header->n_targets); if (bam_is_be) bam_swap_endian_4p(&header->n_targets);
// read reference sequence names and lengths // read reference sequence names and lengths
header->target_name = (char**)calloc(header->n_targets, sizeof(char*)); header->target_name = (char**)xcalloc(header->n_targets, sizeof(char*));
header->target_len = (uint32_t*)calloc(header->n_targets, 4); header->target_len = (uint32_t*)xcalloc(header->n_targets, 4);
for (i = 0; i != header->n_targets; ++i) { for (i = 0; i != header->n_targets; ++i) {
bam_read(fp, &name_len, 4); if (bam_read(fp, &name_len, 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&name_len); if (bam_is_be) bam_swap_endian_4p(&name_len);
header->target_name[i] = (char*)calloc(name_len, 1); header->target_name[i] = (char*)xcalloc(name_len, 1);
bam_read(fp, header->target_name[i], name_len); if (bam_read(fp, header->target_name[i], name_len) != name_len) {
bam_read(fp, &header->target_len[i], 4); goto fail;
}
if (bam_read(fp, &header->target_len[i], 4) != 4) goto fail;
if (bam_is_be) bam_swap_endian_4p(&header->target_len[i]); if (bam_is_be) bam_swap_endian_4p(&header->target_len[i]);
} }
return header; return header;
fail:
bam_header_destroy(header);
return NULL;
} }
static void swap_endian_data(const bam1_core_t *c, int data_len, uint8_t *data) static void swap_endian_data(const bam1_core_t *c, int data_len, uint8_t *data)
@ -146,7 +152,7 @@ int bam_read1(bamFile fp, bam1_t *b)
if (b->m_data < b->data_len) { if (b->m_data < b->data_len) {
b->m_data = b->data_len; b->m_data = b->data_len;
kroundup32(b->m_data); kroundup32(b->m_data);
b->data = (uint8_t*)realloc(b->data, b->m_data); b->data = (uint8_t*)xrealloc(b->data, b->m_data);
} }
if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4; if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4;
b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2; b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;

7
bamlite.h
View File

@ -3,12 +3,13 @@
#include <stdint.h> #include <stdint.h>
#include <zlib.h> #include <zlib.h>
#include "utils.h"
typedef gzFile bamFile; typedef gzFile bamFile;
#define bam_open(fn, mode) gzopen(fn, mode) #define bam_open(fn, mode) xzopen(fn, mode)
#define bam_dopen(fd, mode) gzdopen(fd, mode) #define bam_dopen(fd, mode) gzdopen(fd, mode)
#define bam_close(fp) gzclose(fp) #define bam_close(fp) gzclose(fp)
#define bam_read(fp, buf, size) gzread(fp, buf, size) #define bam_read(fp, buf, size) err_gzread(fp, buf, size)
typedef struct { typedef struct {
int32_t n_targets; int32_t n_targets;
@ -71,7 +72,7 @@ typedef struct {
#define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf) #define bam1_seqi(s, i) ((s)[(i)/2] >> 4*(1-(i)%2) & 0xf)
#define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2) #define bam1_aux(b) ((b)->data + (b)->core.n_cigar*4 + (b)->core.l_qname + (b)->core.l_qseq + ((b)->core.l_qseq + 1)/2)
#define bam_init1() ((bam1_t*)calloc(1, sizeof(bam1_t))) #define bam_init1() ((bam1_t*)xcalloc(1, sizeof(bam1_t)))
#define bam_destroy1(b) do { \ #define bam_destroy1(b) do { \
if (b) { free((b)->data); free(b); } \ if (b) { free((b)->data); free(b); } \
} while (0) } while (0)

101
bntseq.c
View File

@ -30,12 +30,13 @@
#include <string.h> #include <string.h>
#include <zlib.h> #include <zlib.h>
#include <unistd.h> #include <unistd.h>
#include <errno.h>
#include "bntseq.h" #include "bntseq.h"
#include "main.h" #include "main.h"
#include "utils.h" #include "utils.h"
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) KSEQ_INIT(gzFile, err_gzread)
unsigned char nst_nt4_table[256] = { unsigned char nst_nt4_table[256] = {
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,
@ -64,25 +65,25 @@ void bns_dump(const bntseq_t *bns, const char *prefix)
{ // dump .ann { // dump .ann
strcpy(str, prefix); strcat(str, ".ann"); strcpy(str, prefix); strcat(str, ".ann");
fp = xopen(str, "w"); fp = xopen(str, "w");
fprintf(fp, "%lld %d %u\n", (long long)bns->l_pac, bns->n_seqs, bns->seed); err_fprintf(fp, "%lld %d %u\n", (long long)bns->l_pac, bns->n_seqs, bns->seed);
for (i = 0; i != bns->n_seqs; ++i) { for (i = 0; i != bns->n_seqs; ++i) {
bntann1_t *p = bns->anns + i; bntann1_t *p = bns->anns + i;
fprintf(fp, "%d %s", p->gi, p->name); err_fprintf(fp, "%d %s", p->gi, p->name);
if (p->anno[0]) fprintf(fp, " %s\n", p->anno); if (p->anno[0]) err_fprintf(fp, " %s\n", p->anno);
else fprintf(fp, "\n"); else err_fprintf(fp, "\n");
fprintf(fp, "%lld %d %d\n", (long long)p->offset, p->len, p->n_ambs); err_fprintf(fp, "%lld %d %d\n", (long long)p->offset, p->len, p->n_ambs);
} }
fclose(fp); err_fclose(fp);
} }
{ // dump .amb { // dump .amb
strcpy(str, prefix); strcat(str, ".amb"); strcpy(str, prefix); strcat(str, ".amb");
fp = xopen(str, "w"); fp = xopen(str, "w");
fprintf(fp, "%lld %d %u\n", (long long)bns->l_pac, bns->n_seqs, bns->n_holes); err_fprintf(fp, "%lld %d %u\n", (long long)bns->l_pac, bns->n_seqs, bns->n_holes);
for (i = 0; i != bns->n_holes; ++i) { for (i = 0; i != bns->n_holes; ++i) {
bntamb1_t *p = bns->ambs + i; bntamb1_t *p = bns->ambs + i;
fprintf(fp, "%lld %d %c\n", (long long)p->offset, p->len, p->amb); err_fprintf(fp, "%lld %d %c\n", (long long)p->offset, p->len, p->amb);
} }
fclose(fp); err_fclose(fp);
} }
} }
@ -90,53 +91,71 @@ bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, c
{ {
char str[1024]; char str[1024];
FILE *fp; FILE *fp;
const char *fname;
bntseq_t *bns; bntseq_t *bns;
long long xx; long long xx;
int i; int i;
bns = (bntseq_t*)calloc(1, sizeof(bntseq_t)); int scanres;
bns = (bntseq_t*)xcalloc(1, sizeof(bntseq_t));
{ // read .ann { // read .ann
fp = xopen(ann_filename, "r"); fp = xopen(fname = ann_filename, "r");
fscanf(fp, "%lld%d%u", &xx, &bns->n_seqs, &bns->seed); scanres = fscanf(fp, "%lld%d%u", &xx, &bns->n_seqs, &bns->seed);
if (scanres != 3) goto badread;
bns->l_pac = xx; bns->l_pac = xx;
bns->anns = (bntann1_t*)calloc(bns->n_seqs, sizeof(bntann1_t)); bns->anns = (bntann1_t*)xcalloc(bns->n_seqs, sizeof(bntann1_t));
for (i = 0; i < bns->n_seqs; ++i) { for (i = 0; i < bns->n_seqs; ++i) {
bntann1_t *p = bns->anns + i; bntann1_t *p = bns->anns + i;
char *q = str; char *q = str;
int c; int c;
// read gi and sequence name // read gi and sequence name
fscanf(fp, "%u%s", &p->gi, str); scanres = fscanf(fp, "%u%s", &p->gi, str);
p->name = strdup(str); if (scanres != 2) goto badread;
p->name = xstrdup(str);
// read fasta comments // read fasta comments
while ((c = fgetc(fp)) != '\n' && c != EOF) *q++ = c; while (str - q < sizeof(str) - 1 && (c = fgetc(fp)) != '\n' && c != EOF) *q++ = c;
while (c != '\n' && c != EOF) c = fgetc(fp);
if (c == EOF) {
scanres = EOF;
goto badread;
}
*q = 0; *q = 0;
if (q - str > 1) p->anno = strdup(str + 1); // skip leading space if (q - str > 1) p->anno = xstrdup(str + 1); // skip leading space
else p->anno = strdup(""); else p->anno = xstrdup("");
// read the rest // read the rest
fscanf(fp, "%lld%d%d", &xx, &p->len, &p->n_ambs); scanres = fscanf(fp, "%lld%d%d", &xx, &p->len, &p->n_ambs);
if (scanres != 3) goto badread;
p->offset = xx; p->offset = xx;
} }
fclose(fp); err_fclose(fp);
} }
{ // read .amb { // read .amb
int64_t l_pac; int64_t l_pac;
int32_t n_seqs; int32_t n_seqs;
fp = xopen(amb_filename, "r"); fp = xopen(fname = amb_filename, "r");
fscanf(fp, "%lld%d%d", &xx, &n_seqs, &bns->n_holes); scanres = fscanf(fp, "%lld%d%d", &xx, &n_seqs, &bns->n_holes);
if (scanres != 3) goto badread;
l_pac = xx; l_pac = xx;
xassert(l_pac == bns->l_pac && n_seqs == bns->n_seqs, "inconsistent .ann and .amb files."); xassert(l_pac == bns->l_pac && n_seqs == bns->n_seqs, "inconsistent .ann and .amb files.");
bns->ambs = (bntamb1_t*)calloc(bns->n_holes, sizeof(bntamb1_t)); bns->ambs = (bntamb1_t*)xcalloc(bns->n_holes, sizeof(bntamb1_t));
for (i = 0; i < bns->n_holes; ++i) { for (i = 0; i < bns->n_holes; ++i) {
bntamb1_t *p = bns->ambs + i; bntamb1_t *p = bns->ambs + i;
fscanf(fp, "%lld%d%s", &xx, &p->len, str); scanres = fscanf(fp, "%lld%d%s", &xx, &p->len, str);
if (scanres != 3) goto badread;
p->offset = xx; p->offset = xx;
p->amb = str[0]; p->amb = str[0];
} }
fclose(fp); err_fclose(fp);
} }
{ // open .pac { // open .pac
bns->fp_pac = xopen(pac_filename, "rb"); bns->fp_pac = xopen(pac_filename, "rb");
} }
return bns; return bns;
badread:
if (EOF == scanres) {
err_fatal(__func__, "Error reading %s : %s\n", fname, ferror(fp) ? strerror(errno) : "Unexpected end of file");
}
err_fatal(__func__, "Parse error reading %s\n", fname);
} }
bntseq_t *bns_restore(const char *prefix) bntseq_t *bns_restore(const char *prefix)
@ -153,7 +172,7 @@ void bns_destroy(bntseq_t *bns)
if (bns == 0) return; if (bns == 0) return;
else { else {
int i; int i;
if (bns->fp_pac) fclose(bns->fp_pac); if (bns->fp_pac) err_fclose(bns->fp_pac);
free(bns->ambs); free(bns->ambs);
for (i = 0; i < bns->n_seqs; ++i) { for (i = 0; i < bns->n_seqs; ++i) {
free(bns->anns[i].name); free(bns->anns[i].name);
@ -173,11 +192,11 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
int i, lasts; int i, lasts;
if (bns->n_seqs == *m_seqs) { if (bns->n_seqs == *m_seqs) {
*m_seqs <<= 1; *m_seqs <<= 1;
bns->anns = (bntann1_t*)realloc(bns->anns, *m_seqs * sizeof(bntann1_t)); bns->anns = (bntann1_t*)xrealloc(bns->anns, *m_seqs * sizeof(bntann1_t));
} }
p = bns->anns + bns->n_seqs; p = bns->anns + bns->n_seqs;
p->name = strdup((char*)seq->name.s); p->name = xstrdup((char*)seq->name.s);
p->anno = seq->comment.s? strdup((char*)seq->comment.s) : strdup("(null)"); p->anno = seq->comment.s? xstrdup((char*)seq->comment.s) : xstrdup("(null)");
p->gi = 0; p->len = seq->seq.l; p->gi = 0; p->len = seq->seq.l;
p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len; p->offset = (bns->n_seqs == 0)? 0 : (p-1)->offset + (p-1)->len;
p->n_ambs = 0; p->n_ambs = 0;
@ -189,7 +208,7 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
} else { } else {
if (bns->n_holes == *m_holes) { if (bns->n_holes == *m_holes) {
(*m_holes) <<= 1; (*m_holes) <<= 1;
bns->ambs = (bntamb1_t*)realloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t)); bns->ambs = (bntamb1_t*)xrealloc(bns->ambs, (*m_holes) * sizeof(bntamb1_t));
} }
*q = bns->ambs + bns->n_holes; *q = bns->ambs + bns->n_holes;
(*q)->len = 1; (*q)->len = 1;
@ -204,7 +223,7 @@ static uint8_t *add1(const kseq_t *seq, bntseq_t *bns, uint8_t *pac, int64_t *m_
if (c >= 4) c = lrand48()&3; if (c >= 4) c = lrand48()&3;
if (bns->l_pac == *m_pac) { // double the pac size if (bns->l_pac == *m_pac) { // double the pac size
*m_pac <<= 1; *m_pac <<= 1;
pac = realloc(pac, *m_pac/4); pac = xrealloc(pac, *m_pac/4);
memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4); memset(pac + bns->l_pac/4, 0, (*m_pac - bns->l_pac)/4);
} }
_set_pac(pac, bns->l_pac, c); _set_pac(pac, bns->l_pac, c);
@ -229,13 +248,13 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
// initialization // initialization
seq = kseq_init(fp_fa); seq = kseq_init(fp_fa);
bns = (bntseq_t*)calloc(1, sizeof(bntseq_t)); bns = (bntseq_t*)xcalloc(1, sizeof(bntseq_t));
bns->seed = 11; // fixed seed for random generator bns->seed = 11; // fixed seed for random generator
srand48(bns->seed); srand48(bns->seed);
m_seqs = m_holes = 8; m_pac = 0x10000; m_seqs = m_holes = 8; m_pac = 0x10000;
bns->anns = (bntann1_t*)calloc(m_seqs, sizeof(bntann1_t)); bns->anns = (bntann1_t*)xcalloc(m_seqs, sizeof(bntann1_t));
bns->ambs = (bntamb1_t*)calloc(m_holes, sizeof(bntamb1_t)); bns->ambs = (bntamb1_t*)xcalloc(m_holes, sizeof(bntamb1_t));
pac = calloc(m_pac/4, 1); pac = xcalloc(m_pac/4, 1);
q = bns->ambs; q = bns->ambs;
strcpy(name, prefix); strcat(name, ".pac"); strcpy(name, prefix); strcat(name, ".pac");
fp = xopen(name, "wb"); fp = xopen(name, "wb");
@ -243,7 +262,7 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q); while (kseq_read(seq) >= 0) pac = add1(seq, bns, pac, &m_pac, &m_seqs, &m_holes, &q);
if (!for_only) { // add the reverse complemented sequence if (!for_only) { // add the reverse complemented sequence
m_pac = (bns->l_pac * 2 + 3) / 4 * 4; m_pac = (bns->l_pac * 2 + 3) / 4 * 4;
pac = realloc(pac, m_pac/4); pac = xrealloc(pac, m_pac/4);
memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4); memset(pac + (bns->l_pac+3)/4, 0, (m_pac - (bns->l_pac+3)/4*4) / 4);
for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac) for (l = bns->l_pac - 1; l >= 0; --l, ++bns->l_pac)
_set_pac(pac, bns->l_pac, 3-_get_pac(pac, l)); _set_pac(pac, bns->l_pac, 3-_get_pac(pac, l));
@ -251,16 +270,16 @@ int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only)
ret = bns->l_pac; ret = bns->l_pac;
{ // finalize .pac file { // finalize .pac file
ubyte_t ct; ubyte_t ct;
fwrite(pac, 1, (bns->l_pac>>2) + ((bns->l_pac&3) == 0? 0 : 1), fp); err_fwrite(pac, 1, (bns->l_pac>>2) + ((bns->l_pac&3) == 0? 0 : 1), fp);
// the following codes make the pac file size always (l_pac/4+1+1) // the following codes make the pac file size always (l_pac/4+1+1)
if (bns->l_pac % 4 == 0) { if (bns->l_pac % 4 == 0) {
ct = 0; ct = 0;
fwrite(&ct, 1, 1, fp); err_fwrite(&ct, 1, 1, fp);
} }
ct = bns->l_pac % 4; ct = bns->l_pac % 4;
fwrite(&ct, 1, 1, fp); err_fwrite(&ct, 1, 1, fp);
// close .pac file // close .pac file
fclose(fp); err_fclose(fp);
} }
bns_dump(bns, prefix); bns_dump(bns, prefix);
bns_destroy(bns); bns_destroy(bns);

28
bwa.c
View File

@ -1,7 +1,9 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include <math.h> #include <math.h>
#include "utils.h"
#include "bwa.h" #include "bwa.h"
#include "bwt.h" #include "bwt.h"
#include "bwtgap.h" #include "bwtgap.h"
@ -38,8 +40,8 @@ bwa_idx_t *bwa_idx_load(const char *prefix)
int l; int l;
char *str; char *str;
l = strlen(prefix); l = strlen(prefix);
p = calloc(1, sizeof(bwa_idx_t)); p = xcalloc(1, sizeof(bwa_idx_t));
str = malloc(l + 10); str = xmalloc(l + 10);
strcpy(str, prefix); strcpy(str, prefix);
p->bns = bns_restore(str); p->bns = bns_restore(str);
strcpy(str + l, ".bwt"); strcpy(str + l, ".bwt");
@ -48,9 +50,9 @@ bwa_idx_t *bwa_idx_load(const char *prefix)
strcpy(str + l, ".sa"); strcpy(str + l, ".sa");
bwt_restore_sa(str, p->bwt); bwt_restore_sa(str, p->bwt);
free(str); free(str);
p->pac = calloc(p->bns->l_pac/4+1, 1); p->pac = xcalloc(p->bns->l_pac/4+1, 1);
fread(p->pac, 1, p->bns->l_pac/4+1, p->bns->fp_pac); err_fread_noeof(p->pac, 1, p->bns->l_pac/4+1, p->bns->fp_pac);
fclose(p->bns->fp_pac); err_fclose(p->bns->fp_pac);
p->bns->fp_pac = 0; p->bns->fp_pac = 0;
return p; return p;
} }
@ -69,7 +71,7 @@ bwa_buf_t *bwa_buf_init(const bwa_opt_t *opt, int max_score)
extern int bwa_cal_maxdiff(int l, double err, double thres); extern int bwa_cal_maxdiff(int l, double err, double thres);
int i; int i;
bwa_buf_t *p; bwa_buf_t *p;
p = malloc(sizeof(bwa_buf_t)); p = xmalloc(sizeof(bwa_buf_t));
p->stack = gap_init_stack2(max_score); p->stack = gap_init_stack2(max_score);
p->opt = gap_init_opt(); p->opt = gap_init_opt();
p->opt->s_gapo = opt->s_gapo; p->opt->s_gapo = opt->s_gapo;
@ -80,10 +82,10 @@ bwa_buf_t *bwa_buf_init(const bwa_opt_t *opt, int max_score)
p->opt->seed_len = opt->seed_len; p->opt->seed_len = opt->seed_len;
p->opt->max_seed_diff = opt->max_seed_diff; p->opt->max_seed_diff = opt->max_seed_diff;
p->opt->fnr = opt->fnr; p->opt->fnr = opt->fnr;
p->diff_tab = calloc(BWA_MAX_QUERY_LEN, sizeof(int)); p->diff_tab = xcalloc(BWA_MAX_QUERY_LEN, sizeof(int));
for (i = 1; i < BWA_MAX_QUERY_LEN; ++i) for (i = 1; i < BWA_MAX_QUERY_LEN; ++i)
p->diff_tab[i] = bwa_cal_maxdiff(i, BWA_AVG_ERR, opt->fnr); p->diff_tab[i] = bwa_cal_maxdiff(i, BWA_AVG_ERR, opt->fnr);
p->logn = calloc(256, sizeof(int)); p->logn = xcalloc(256, sizeof(int));
for (i = 1; i != 256; ++i) for (i = 1; i != 256; ++i)
p->logn[i] = (int)(4.343 * log(i) + 0.499); p->logn[i] = (int)(4.343 * log(i) + 0.499);
return p; return p;
@ -111,7 +113,7 @@ bwa_sai_t bwa_sai(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq)
if (buf_len > buf->max_buf) { if (buf_len > buf->max_buf) {
buf->max_buf = buf_len; buf->max_buf = buf_len;
kroundup32(buf->max_buf); kroundup32(buf->max_buf);
buf->buf = realloc(buf->buf, buf->max_buf); buf->buf = xrealloc(buf->buf, buf->max_buf);
} }
memset(buf->buf, 0, buf_len); memset(buf->buf, 0, buf_len);
seed_w = (bwt_width_t*)buf->buf; seed_w = (bwt_width_t*)buf->buf;
@ -166,7 +168,7 @@ void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t
if (seq_len<<1 > buf->max_buf) { if (seq_len<<1 > buf->max_buf) {
buf->max_buf = seq_len<<1; buf->max_buf = seq_len<<1;
kroundup32(buf->max_buf); kroundup32(buf->max_buf);
buf->buf = realloc(buf->buf, buf->max_buf); buf->buf = xrealloc(buf->buf, buf->max_buf);
} }
s[0] = buf->buf; s[0] = buf->buf;
s[1] = s[0] + seq_len; s[1] = s[0] + seq_len;
@ -180,7 +182,7 @@ void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t
bwa_cigar_t *cigar16; bwa_cigar_t *cigar16;
cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, seq_len, s[strand], &pac_pos, strand? n_gaps : -n_gaps, &n_cigar, 1); cigar16 = bwa_refine_gapped_core(idx->bns->l_pac, idx->pac, seq_len, s[strand], &pac_pos, strand? n_gaps : -n_gaps, &n_cigar, 1);
aln->n_cigar = n_cigar; aln->n_cigar = n_cigar;
aln->cigar = malloc(n_cigar * 4); aln->cigar = xmalloc(n_cigar * 4);
for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) { for (i = 0, pos3 = pac_pos; i < n_cigar; ++i) {
int op = cigar16[i]>>14; int op = cigar16[i]>>14;
int len = cigar16[i]&0x3fff; int len = cigar16[i]&0x3fff;
@ -191,7 +193,7 @@ void bwa_sa2aln(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, uint64_t
free(cigar16); free(cigar16);
} else { // ungapped } else { // ungapped
aln->n_cigar = 1; aln->n_cigar = 1;
aln->cigar = malloc(4); aln->cigar = xmalloc(4);
aln->cigar[0] = seq_len<<4 | 0; aln->cigar[0] = seq_len<<4 | 0;
pos3 = pac_pos + seq_len; pos3 = pac_pos + seq_len;
} }
@ -214,7 +216,7 @@ bwa_one_t *bwa_se(const bwa_idx_t *idx, bwa_buf_t *buf, const char *seq, int gen
int best, cnt, i, seq_len; int best, cnt, i, seq_len;
seq_len = strlen(seq); seq_len = strlen(seq);
one = calloc(1, sizeof(bwa_one_t)); one = xcalloc(1, sizeof(bwa_one_t));
one->sai = bwa_sai(idx, buf, seq); one->sai = bwa_sai(idx, buf, seq);
if (one->sai.n == 0) return one; if (one->sai.n == 0) return one;
// count number of hits; randomly select one alignment // count number of hits; randomly select one alignment

40
bwape.c
View File

@ -58,7 +58,7 @@ void bwa_print_sam_PG();
pe_opt_t *bwa_init_pe_opt() pe_opt_t *bwa_init_pe_opt()
{ {
pe_opt_t *po; pe_opt_t *po;
po = (pe_opt_t*)calloc(1, sizeof(pe_opt_t)); po = (pe_opt_t*)xcalloc(1, sizeof(pe_opt_t));
po->max_isize = 500; po->max_isize = 500;
po->force_isize = 0; po->force_isize = 0;
po->max_occ = 100000; po->max_occ = 100000;
@ -104,7 +104,7 @@ static int infer_isize(int n_seqs, bwa_seq_t *seqs[2], isize_info_t *ii, double
ii->avg = ii->std = -1.0; ii->avg = ii->std = -1.0;
ii->low = ii->high = ii->high_bayesian = 0; ii->low = ii->high = ii->high_bayesian = 0;
isizes = (uint64_t*)calloc(n_seqs, 8); isizes = (uint64_t*)xcalloc(n_seqs, 8);
for (i = 0, tot = 0; i != n_seqs; ++i) { for (i = 0, tot = 0; i != n_seqs; ++i) {
bwa_seq_t *p[2]; bwa_seq_t *p[2];
p[0] = seqs[0] + i; p[1] = seqs[1] + i; p[0] = seqs[0] + i; p[1] = seqs[1] + i;
@ -292,9 +292,9 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
pe_data_t *d; pe_data_t *d;
aln_buf_t *buf[2]; aln_buf_t *buf[2];
d = (pe_data_t*)calloc(1, sizeof(pe_data_t)); d = (pe_data_t*)xcalloc(1, sizeof(pe_data_t));
buf[0] = (aln_buf_t*)calloc(n_seqs, sizeof(aln_buf_t)); buf[0] = (aln_buf_t*)xcalloc(n_seqs, sizeof(aln_buf_t));
buf[1] = (aln_buf_t*)calloc(n_seqs, sizeof(aln_buf_t)); buf[1] = (aln_buf_t*)xcalloc(n_seqs, sizeof(aln_buf_t));
if (_bwt == 0) { // load forward SA if (_bwt == 0) { // load forward SA
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
@ -309,11 +309,11 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
p[j] = seqs[j] + i; p[j] = seqs[j] + i;
p[j]->n_multi = 0; p[j]->n_multi = 0;
p[j]->extra_flag |= SAM_FPD | (j == 0? SAM_FR1 : SAM_FR2); p[j]->extra_flag |= SAM_FPD | (j == 0? SAM_FR1 : SAM_FR2);
fread(&n_aln, 4, 1, fp_sa[j]); err_fread_noeof(&n_aln, 4, 1, fp_sa[j]);
if (n_aln > kv_max(d->aln[j])) if (n_aln > kv_max(d->aln[j]))
kv_resize(bwt_aln1_t, d->aln[j], n_aln); kv_resize(bwt_aln1_t, d->aln[j], n_aln);
d->aln[j].n = n_aln; d->aln[j].n = n_aln;
fread(d->aln[j].a, sizeof(bwt_aln1_t), n_aln, fp_sa[j]); err_fread_noeof(d->aln[j].a, sizeof(bwt_aln1_t), n_aln, fp_sa[j]);
kv_copy(bwt_aln1_t, buf[j][i].aln, d->aln[j]); // backup d->aln[j] kv_copy(bwt_aln1_t, buf[j][i].aln, d->aln[j]); // backup d->aln[j]
// generate SE alignment and mapping quality // generate SE alignment and mapping quality
bwa_aln2seq(n_aln, d->aln[j].a, p[j]); bwa_aln2seq(n_aln, d->aln[j].a, p[j]);
@ -367,7 +367,7 @@ int bwa_cal_pac_pos_pe(const bntseq_t *bns, const char *prefix, bwt_t *const _bw
if (ret) { // not in the hash table; ret must equal 1 as we never remove elements if (ret) { // not in the hash table; ret must equal 1 as we never remove elements
poslist_t *z = &kh_val(g_hash, iter); poslist_t *z = &kh_val(g_hash, iter);
z->n = r->l - r->k + 1; z->n = r->l - r->k + 1;
z->a = (bwtint_t*)malloc(sizeof(bwtint_t) * z->n); z->a = (bwtint_t*)xmalloc(sizeof(bwtint_t) * z->n);
for (l = r->k; l <= r->l; ++l) { for (l = r->k; l <= r->l; ++l) {
int strand; int strand;
z->a[l - r->k] = bwa_sa2pos(bns, bwt, l, p[j]->len, &strand)<<1; z->a[l - r->k] = bwa_sa2pos(bns, bwt, l, p[j]->len, &strand)<<1;
@ -448,10 +448,10 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u
if ((float)x/len >= 0.25 || len - x < SW_MIN_MATCH_LEN) return 0; if ((float)x/len >= 0.25 || len - x < SW_MIN_MATCH_LEN) return 0;
// get reference subsequence // get reference subsequence
ref_seq = (ubyte_t*)calloc(reglen, 1); ref_seq = (ubyte_t*)xcalloc(reglen, 1);
for (k = *beg, l = 0; l < reglen && k < l_pac; ++k) for (k = *beg, l = 0; l < reglen && k < l_pac; ++k)
ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3;
path = (path_t*)calloc(l+len, sizeof(path_t)); path = (path_t*)xcalloc(l+len, sizeof(path_t));
// do alignment // do alignment
ret = aln_local_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len, 1, &subo); ret = aln_local_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len, 1, &subo);
@ -480,7 +480,7 @@ bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const u
*beg += (p->i? p->i : 1) - 1; *beg += (p->i? p->i : 1) - 1;
start = (p->j? p->j : 1) - 1; start = (p->j? p->j : 1) - 1;
end = path->j; end = path->j;
cigar = (bwa_cigar_t*)realloc(cigar, sizeof(bwa_cigar_t) * (*n_cigar + 2)); cigar = (bwa_cigar_t*)xrealloc(cigar, sizeof(bwa_cigar_t) * (*n_cigar + 2));
if (start) { if (start) {
memmove(cigar + 1, cigar, sizeof(bwa_cigar_t) * (*n_cigar)); memmove(cigar + 1, cigar, sizeof(bwa_cigar_t) * (*n_cigar));
cigar[0] = __cigar_create(3, start); cigar[0] = __cigar_create(3, start);
@ -525,9 +525,9 @@ ubyte_t *bwa_paired_sw(const bntseq_t *bns, const ubyte_t *_pacseq, int n_seqs,
// load reference sequence // load reference sequence
if (_pacseq == 0) { if (_pacseq == 0) {
pacseq = (ubyte_t*)calloc(bns->l_pac/4+1, 1); pacseq = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1);
rewind(bns->fp_pac); err_rewind(bns->fp_pac);
fread(pacseq, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac);
} else pacseq = (ubyte_t*)_pacseq; } else pacseq = (ubyte_t*)_pacseq;
if (!popt->is_sw || ii->avg < 0.0) return pacseq; if (!popt->is_sw || ii->avg < 0.0) return pacseq;
@ -683,10 +683,10 @@ void bwa_sai2sam_pe_core(const char *prefix, char *const fn_sa[2], char *const f
g_hash = kh_init(b128); g_hash = kh_init(b128);
last_ii.avg = -1.0; last_ii.avg = -1.0;
fread(&opt, sizeof(gap_opt_t), 1, fp_sa[0]); err_fread_noeof(&opt, sizeof(gap_opt_t), 1, fp_sa[0]);
ks[0] = bwa_open_reads(opt.mode, fn_fa[0]); ks[0] = bwa_open_reads(opt.mode, fn_fa[0]);
opt0 = opt; opt0 = opt;
fread(&opt, sizeof(gap_opt_t), 1, fp_sa[1]); // overwritten! err_fread_noeof(&opt, sizeof(gap_opt_t), 1, fp_sa[1]); // overwritten!
ks[1] = bwa_open_reads(opt.mode, fn_fa[1]); ks[1] = bwa_open_reads(opt.mode, fn_fa[1]);
if (!(opt.mode & BWA_MODE_COMPREAD)) { if (!(opt.mode & BWA_MODE_COMPREAD)) {
popt->type = BWA_PET_SOLID; popt->type = BWA_PET_SOLID;
@ -695,9 +695,9 @@ void bwa_sai2sam_pe_core(const char *prefix, char *const fn_sa[2], char *const f
if (popt->is_preload) { if (popt->is_preload) {
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt); strcpy(str, prefix); strcat(str, ".sa"); bwt_restore_sa(str, bwt);
pac = (ubyte_t*)calloc(bns->l_pac/4+1, 1); pac = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1);
rewind(bns->fp_pac); err_rewind(bns->fp_pac);
fread(pac, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
} }
} }
@ -752,7 +752,7 @@ void bwa_sai2sam_pe_core(const char *prefix, char *const fn_sa[2], char *const f
if (ntbns) bns_destroy(ntbns); if (ntbns) bns_destroy(ntbns);
for (i = 0; i < 2; ++i) { for (i = 0; i < 2; ++i) {
bwa_seq_close(ks[i]); bwa_seq_close(ks[i]);
fclose(fp_sa[i]); err_fclose(fp_sa[i]);
} }
for (iter = kh_begin(g_hash); iter != kh_end(g_hash); ++iter) for (iter = kh_begin(g_hash); iter != kh_end(g_hash); ++iter)
if (kh_exist(g_hash, iter)) free(kh_val(g_hash, iter).a); if (kh_exist(g_hash, iter)) free(kh_val(g_hash, iter).a);

48
bwase.c
View File

@ -59,7 +59,7 @@ void bwa_aln2seq_core(int n_aln, const bwt_aln1_t *aln, bwa_seq_t *s, int set_ma
* simply output all hits, but the following samples "rest" * simply output all hits, but the following samples "rest"
* number of random hits. */ * number of random hits. */
rest = n_occ > n_multi + 1? n_multi + 1 : n_occ; // find one additional for ->sa rest = n_occ > n_multi + 1? n_multi + 1 : n_occ; // find one additional for ->sa
s->multi = calloc(rest, sizeof(bwt_multi1_t)); s->multi = xcalloc(rest, sizeof(bwt_multi1_t));
for (k = 0; k < n_aln; ++k) { for (k = 0; k < n_aln; ++k) {
const bwt_aln1_t *q = aln + k; const bwt_aln1_t *q = aln + k;
if (q->l - q->k + 1 <= rest) { if (q->l - q->k + 1 <= rest) {
@ -172,16 +172,16 @@ bwa_cigar_t *bwa_refine_gapped_core(bwtint_t l_pac, const ubyte_t *pacseq, int l
ref_len = len + abs(ext); ref_len = len + abs(ext);
if (ext > 0) { if (ext > 0) {
ref_seq = (ubyte_t*)calloc(ref_len, 1); ref_seq = (ubyte_t*)xcalloc(ref_len, 1);
for (k = __pos; k < __pos + ref_len && k < l_pac; ++k) for (k = __pos; k < __pos + ref_len && k < l_pac; ++k)
ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3;
} else { } else {
int64_t x = __pos + (is_end_correct? len : ref_len); int64_t x = __pos + (is_end_correct? len : ref_len);
ref_seq = (ubyte_t*)calloc(ref_len, 1); ref_seq = (ubyte_t*)xcalloc(ref_len, 1);
for (l = 0, k = x - ref_len > 0? x - ref_len : 0; k < x && k < l_pac; ++k) for (l = 0, k = x - ref_len > 0? x - ref_len : 0; k < x && k < l_pac; ++k)
ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3; ref_seq[l++] = pacseq[k>>2] >> ((~k&3)<<1) & 3;
} }
path = (path_t*)calloc(l+len, sizeof(path_t)); path = (path_t*)xcalloc(l+len, sizeof(path_t));
aln_global_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len); aln_global_core(ref_seq, l, (ubyte_t*)seq, len, &ap, path, &path_len);
cigar = bwa_aln_path2cigar(path, path_len, n_cigar); cigar = bwa_aln_path2cigar(path, path_len, n_cigar);
@ -257,7 +257,7 @@ char *bwa_cal_md1(int n_cigar, bwa_cigar_t *cigar, int len, bwtint_t pos, ubyte_
} }
ksprintf(str, "%d", u); ksprintf(str, "%d", u);
*_nm = nm; *_nm = nm;
return strdup(str->s); return xstrdup(str->s);
} }
void bwa_correct_trimmed(bwa_seq_t *s) void bwa_correct_trimmed(bwa_seq_t *s)
@ -269,11 +269,11 @@ void bwa_correct_trimmed(bwa_seq_t *s)
} else { } else {
if (s->cigar == 0) { if (s->cigar == 0) {
s->n_cigar = 2; s->n_cigar = 2;
s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t)); s->cigar = xcalloc(s->n_cigar, sizeof(bwa_cigar_t));
s->cigar[0] = __cigar_create(0, s->len); s->cigar[0] = __cigar_create(0, s->len);
} else { } else {
++s->n_cigar; ++s->n_cigar;
s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); s->cigar = xrealloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t));
} }
s->cigar[s->n_cigar-1] = __cigar_create(3, (s->full_len - s->len)); s->cigar[s->n_cigar-1] = __cigar_create(3, (s->full_len - s->len));
} }
@ -283,11 +283,11 @@ void bwa_correct_trimmed(bwa_seq_t *s)
} else { } else {
if (s->cigar == 0) { if (s->cigar == 0) {
s->n_cigar = 2; s->n_cigar = 2;
s->cigar = calloc(s->n_cigar, sizeof(bwa_cigar_t)); s->cigar = xcalloc(s->n_cigar, sizeof(bwa_cigar_t));
s->cigar[1] = __cigar_create(0, s->len); s->cigar[1] = __cigar_create(0, s->len);
} else { } else {
++s->n_cigar; ++s->n_cigar;
s->cigar = realloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t)); s->cigar = xrealloc(s->cigar, s->n_cigar * sizeof(bwa_cigar_t));
memmove(s->cigar + 1, s->cigar, (s->n_cigar-1) * sizeof(bwa_cigar_t)); memmove(s->cigar + 1, s->cigar, (s->n_cigar-1) * sizeof(bwa_cigar_t));
} }
s->cigar[0] = __cigar_create(3, (s->full_len - s->len)); s->cigar[0] = __cigar_create(3, (s->full_len - s->len));
@ -303,15 +303,15 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t
kstring_t *str; kstring_t *str;
if (ntbns) { // in color space if (ntbns) { // in color space
ntpac = (ubyte_t*)calloc(ntbns->l_pac/4+1, 1); ntpac = (ubyte_t*)xcalloc(ntbns->l_pac/4+1, 1);
rewind(ntbns->fp_pac); err_rewind(ntbns->fp_pac);
fread(ntpac, 1, ntbns->l_pac/4 + 1, ntbns->fp_pac); err_fread_noeof(ntpac, 1, ntbns->l_pac/4 + 1, ntbns->fp_pac);
} }
if (!_pacseq) { if (!_pacseq) {
pacseq = (ubyte_t*)calloc(bns->l_pac/4+1, 1); pacseq = (ubyte_t*)xcalloc(bns->l_pac/4+1, 1);
rewind(bns->fp_pac); err_rewind(bns->fp_pac);
fread(pacseq, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pacseq, 1, bns->l_pac/4+1, bns->fp_pac);
} else pacseq = _pacseq; } else pacseq = _pacseq;
for (i = 0; i != n_seqs; ++i) { for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *s = seqs + i; bwa_seq_t *s = seqs + i;
@ -351,7 +351,7 @@ void bwa_refine_gapped(const bntseq_t *bns, int n_seqs, bwa_seq_t *seqs, ubyte_t
} }
#endif #endif
// generate MD tag // generate MD tag
str = (kstring_t*)calloc(1, sizeof(kstring_t)); str = (kstring_t*)xcalloc(1, sizeof(kstring_t));
for (i = 0; i != n_seqs; ++i) { for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *s = seqs + i; bwa_seq_t *s = seqs + i;
if (s->type != BWA_TYPE_NO_MATCH) { if (s->type != BWA_TYPE_NO_MATCH) {
@ -523,7 +523,7 @@ bntseq_t *bwa_open_nt(const char *prefix)
{ {
bntseq_t *ntbns; bntseq_t *ntbns;
char *str; char *str;
str = (char*)calloc(strlen(prefix) + 10, 1); str = (char*)xcalloc(strlen(prefix) + 10, 1);
strcat(strcpy(str, prefix), ".nt"); strcat(strcpy(str, prefix), ".nt");
ntbns = bns_restore(str); ntbns = bns_restore(str);
free(str); free(str);
@ -566,14 +566,14 @@ int bwa_set_rg(const char *s)
if (strstr(s, "@RG") != s) return -1; if (strstr(s, "@RG") != s) return -1;
if (bwa_rg_line) free(bwa_rg_line); if (bwa_rg_line) free(bwa_rg_line);
if (bwa_rg_id) free(bwa_rg_id); if (bwa_rg_id) free(bwa_rg_id);
bwa_rg_line = strdup(s); bwa_rg_line = xstrdup(s);
bwa_rg_id = 0; bwa_rg_id = 0;
bwa_escape(bwa_rg_line); bwa_escape(bwa_rg_line);
p = strstr(bwa_rg_line, "\tID:"); p = strstr(bwa_rg_line, "\tID:");
if (p == 0) return -1; if (p == 0) return -1;
p += 4; p += 4;
for (q = p; *q && *q != '\t' && *q != '\n'; ++q); for (q = p; *q && *q != '\t' && *q != '\n'; ++q);
bwa_rg_id = calloc(q - p + 1, 1); bwa_rg_id = xcalloc(q - p + 1, 1);
for (q = p, r = bwa_rg_id; *q && *q != '\t' && *q != '\n'; ++q) for (q = p, r = bwa_rg_id; *q && *q != '\t' && *q != '\n'; ++q)
*r++ = *q; *r++ = *q;
return 0; return 0;
@ -598,7 +598,7 @@ void bwa_sai2sam_se_core(const char *prefix, const char *fn_sa, const char *fn_f
fp_sa = xopen(fn_sa, "r"); fp_sa = xopen(fn_sa, "r");
m_aln = 0; m_aln = 0;
fread(&opt, sizeof(gap_opt_t), 1, fp_sa); err_fread_noeof(&opt, sizeof(gap_opt_t), 1, fp_sa);
if (!(opt.mode & BWA_MODE_COMPREAD)) // in color space; initialize ntpac if (!(opt.mode & BWA_MODE_COMPREAD)) // in color space; initialize ntpac
ntbns = bwa_open_nt(prefix); ntbns = bwa_open_nt(prefix);
bwa_print_sam_SQ(bns); bwa_print_sam_SQ(bns);
@ -614,12 +614,12 @@ void bwa_sai2sam_se_core(const char *prefix, const char *fn_sa, const char *fn_f
for (i = 0; i < n_seqs; ++i) { for (i = 0; i < n_seqs; ++i) {
bwa_seq_t *p = seqs + i; bwa_seq_t *p = seqs + i;
int n_aln; int n_aln;
fread(&n_aln, 4, 1, fp_sa); err_fread_noeof(&n_aln, 4, 1, fp_sa);
if (n_aln > m_aln) { if (n_aln > m_aln) {
m_aln = n_aln; m_aln = n_aln;
aln = (bwt_aln1_t*)realloc(aln, sizeof(bwt_aln1_t) * m_aln); aln = (bwt_aln1_t*)xrealloc(aln, sizeof(bwt_aln1_t) * m_aln);
} }
fread(aln, sizeof(bwt_aln1_t), n_aln, fp_sa); err_fread_noeof(aln, sizeof(bwt_aln1_t), n_aln, fp_sa);
bwa_aln2seq_core(n_aln, aln, p, 1, n_occ); bwa_aln2seq_core(n_aln, aln, p, 1, n_occ);
} }
@ -644,7 +644,7 @@ void bwa_sai2sam_se_core(const char *prefix, const char *fn_sa, const char *fn_f
bwa_seq_close(ks); bwa_seq_close(ks);
if (ntbns) bns_destroy(ntbns); if (ntbns) bns_destroy(ntbns);
bns_destroy(bns); bns_destroy(bns);
fclose(fp_sa); err_fclose(fp_sa);
free(aln); free(aln);
} }

26
bwaseqio.c
View File

@ -5,7 +5,7 @@
#include "bamlite.h" #include "bamlite.h"
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) KSEQ_INIT(gzFile, err_gzread)
extern unsigned char nst_nt4_table[256]; extern unsigned char nst_nt4_table[256];
static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
@ -22,7 +22,7 @@ bwa_seqio_t *bwa_bam_open(const char *fn, int which)
{ {
bwa_seqio_t *bs; bwa_seqio_t *bs;
bam_header_t *h; bam_header_t *h;
bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t)); bs = (bwa_seqio_t*)xcalloc(1, sizeof(bwa_seqio_t));
bs->is_bam = 1; bs->is_bam = 1;
bs->which = which; bs->which = which;
bs->fp = bam_open(fn, "r"); bs->fp = bam_open(fn, "r");
@ -35,7 +35,7 @@ bwa_seqio_t *bwa_seq_open(const char *fn)
{ {
gzFile fp; gzFile fp;
bwa_seqio_t *bs; bwa_seqio_t *bs;
bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t)); bs = (bwa_seqio_t*)xcalloc(1, sizeof(bwa_seqio_t));
fp = xzopen(fn, "r"); fp = xzopen(fn, "r");
bs->ks = kseq_init(fp); bs->ks = kseq_init(fp);
return bs; return bs;
@ -93,7 +93,7 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
b = bam_init1(); b = bam_init1();
n_seqs = 0; n_seqs = 0;
seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t)); seqs = (bwa_seq_t*)xcalloc(n_needed, sizeof(bwa_seq_t));
while (bam_read1(bs->fp, b) >= 0) { while (bam_read1(bs->fp, b) >= 0) {
uint8_t *s, *q; uint8_t *s, *q;
int go = 0; int go = 0;
@ -108,8 +108,8 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
p->full_len = p->clip_len = p->len = l; p->full_len = p->clip_len = p->len = l;
n_tot += p->full_len; n_tot += p->full_len;
s = bam1_seq(b); q = bam1_qual(b); s = bam1_seq(b); q = bam1_qual(b);
p->seq = (ubyte_t*)calloc(p->len + 1, 1); p->seq = (ubyte_t*)xcalloc(p->len + 1, 1);
p->qual = (ubyte_t*)calloc(p->len + 1, 1); p->qual = (ubyte_t*)xcalloc(p->len + 1, 1);
for (i = 0; i != p->full_len; ++i) { for (i = 0; i != p->full_len; ++i) {
p->seq[i] = bam_nt16_nt4_table[(int)bam1_seqi(s, i)]; p->seq[i] = bam_nt16_nt4_table[(int)bam1_seqi(s, i)];
p->qual[i] = q[i] + 33 < 126? q[i] + 33 : 126; p->qual[i] = q[i] + 33 < 126? q[i] + 33 : 126;
@ -119,11 +119,11 @@ static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_com
seq_reverse(p->len, p->qual, 0); seq_reverse(p->len, p->qual, 0);
} }
if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p);
p->rseq = (ubyte_t*)calloc(p->full_len, 1); p->rseq = (ubyte_t*)xcalloc(p->full_len, 1);
memcpy(p->rseq, p->seq, p->len); memcpy(p->rseq, p->seq, p->len);
seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped()
seq_reverse(p->len, p->rseq, is_comp); seq_reverse(p->len, p->rseq, is_comp);
p->name = strdup((const char*)bam1_qname(b)); p->name = xstrdup((const char*)bam1_qname(b));
if (n_seqs == n_needed) break; if (n_seqs == n_needed) break;
} }
*n = n_seqs; *n = n_seqs;
@ -153,7 +153,7 @@ bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int tri
} }
if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input
n_seqs = 0; n_seqs = 0;
seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t)); seqs = (bwa_seq_t*)xcalloc(n_needed, sizeof(bwa_seq_t));
while ((l = kseq_read(seq)) >= 0) { while ((l = kseq_read(seq)) >= 0) {
if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) { if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) {
// skip reads that are marked to be filtered by Casava // skip reads that are marked to be filtered by Casava
@ -184,18 +184,18 @@ bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int tri
p->qual = 0; p->qual = 0;
p->full_len = p->clip_len = p->len = l; p->full_len = p->clip_len = p->len = l;
n_tot += p->full_len; n_tot += p->full_len;
p->seq = (ubyte_t*)calloc(p->len, 1); p->seq = (ubyte_t*)xcalloc(p->full_len, 1);
for (i = 0; i != p->full_len; ++i) for (i = 0; i != p->full_len; ++i)
p->seq[i] = nst_nt4_table[(int)seq->seq.s[i]]; p->seq[i] = nst_nt4_table[(int)seq->seq.s[i]];
if (seq->qual.l) { // copy quality if (seq->qual.l) { // copy quality
p->qual = (ubyte_t*)strdup((char*)seq->qual.s); p->qual = (ubyte_t*)xstrdup((char*)seq->qual.s);
if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p);
} }
p->rseq = (ubyte_t*)calloc(p->full_len, 1); p->rseq = (ubyte_t*)xcalloc(p->full_len, 1);
memcpy(p->rseq, p->seq, p->len); memcpy(p->rseq, p->seq, p->len);
seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped()
seq_reverse(p->len, p->rseq, is_comp); seq_reverse(p->len, p->rseq, is_comp);
p->name = strdup((const char*)seq->name.s); p->name = xstrdup((const char*)seq->name.s);
{ // trim /[12]$ { // trim /[12]$
int t = strlen(p->name); int t = strlen(p->name);
if (t > 2 && p->name[t-2] == '/' && (p->name[t-1] == '1' || p->name[t-1] == '2')) p->name[t-2] = '\0'; if (t > 2 && p->name[t-2] == '/' && (p->name[t-1] == '1' || p->name[t-1] == '2')) p->name[t-2] = '\0';

6
bwt.c
View File

@ -58,11 +58,7 @@ void bwt_cal_sa(bwt_t *bwt, int intv)
if (bwt->sa) free(bwt->sa); if (bwt->sa) free(bwt->sa);
bwt->sa_intv = intv; bwt->sa_intv = intv;
bwt->n_sa = (bwt->seq_len + intv) / intv; bwt->n_sa = (bwt->seq_len + intv) / intv;
bwt->sa = (bwtint_t*)calloc(bwt->n_sa, sizeof(bwtint_t)); bwt->sa = (bwtint_t*)xcalloc(bwt->n_sa, sizeof(bwtint_t));
if (bwt->sa == 0) {
fprintf(stderr, "[%s] Fail to allocate %.3fMB memory. Abort!\n", __func__, bwt->n_sa * sizeof(bwtint_t) / 1024.0/1024.0);
abort();
}
// calculate SA value // calculate SA value
isa = 0; sa = bwt->seq_len; isa = 0; sa = bwt->seq_len;
for (i = 0; i < bwt->seq_len; ++i) { for (i = 0; i < bwt->seq_len; ++i) {

43
bwt_gen.c
View File

@ -28,6 +28,7 @@
#include <assert.h> #include <assert.h>
#include <stdint.h> #include <stdint.h>
#include "QSufSort.h" #include "QSufSort.h"
#include "utils.h"
typedef uint64_t bgint_t; typedef uint64_t bgint_t;
typedef int64_t sbgint_t; typedef int64_t sbgint_t;
@ -319,25 +320,25 @@ BWT *BWTCreate(const bgint_t textLength, unsigned int *decodeTable)
{ {
BWT *bwt; BWT *bwt;
bwt = (BWT*)calloc(1, sizeof(BWT)); bwt = (BWT*)xcalloc(1, sizeof(BWT));
bwt->textLength = 0; bwt->textLength = 0;
bwt->cumulativeFreq = (bgint_t*)calloc((ALPHABET_SIZE + 1), sizeof(bgint_t)); bwt->cumulativeFreq = (bgint_t*)xcalloc((ALPHABET_SIZE + 1), sizeof(bgint_t));
initializeVAL_bg(bwt->cumulativeFreq, ALPHABET_SIZE + 1, 0); initializeVAL_bg(bwt->cumulativeFreq, ALPHABET_SIZE + 1, 0);
bwt->bwtSizeInWord = 0; bwt->bwtSizeInWord = 0;
// Generate decode tables // Generate decode tables
if (decodeTable == NULL) { if (decodeTable == NULL) {
bwt->decodeTable = (unsigned*)calloc(DNA_OCC_CNT_TABLE_SIZE_IN_WORD, sizeof(unsigned int)); bwt->decodeTable = (unsigned*)xcalloc(DNA_OCC_CNT_TABLE_SIZE_IN_WORD, sizeof(unsigned int));
GenerateDNAOccCountTable(bwt->decodeTable); GenerateDNAOccCountTable(bwt->decodeTable);
} else { } else {
bwt->decodeTable = decodeTable; bwt->decodeTable = decodeTable;
} }
bwt->occMajorSizeInWord = BWTOccValueMajorSizeInWord(textLength); bwt->occMajorSizeInWord = BWTOccValueMajorSizeInWord(textLength);
bwt->occValueMajor = (bgint_t*)calloc(bwt->occMajorSizeInWord, sizeof(bgint_t)); bwt->occValueMajor = (bgint_t*)xcalloc(bwt->occMajorSizeInWord, sizeof(bgint_t));
bwt->occSizeInWord = 0; bwt->occSizeInWord = 0;
bwt->occValue = NULL; bwt->occValue = NULL;
@ -353,16 +354,16 @@ BWTInc *BWTIncCreate(const bgint_t textLength, unsigned int initialMaxBuildSize,
if (textLength < incMaxBuildSize) incMaxBuildSize = textLength; if (textLength < incMaxBuildSize) incMaxBuildSize = textLength;
if (textLength < initialMaxBuildSize) initialMaxBuildSize = textLength; if (textLength < initialMaxBuildSize) initialMaxBuildSize = textLength;
bwtInc = (BWTInc*)calloc(1, sizeof(BWTInc)); bwtInc = (BWTInc*)xcalloc(1, sizeof(BWTInc));
bwtInc->numberOfIterationDone = 0; bwtInc->numberOfIterationDone = 0;
bwtInc->bwt = BWTCreate(textLength, NULL); bwtInc->bwt = BWTCreate(textLength, NULL);
bwtInc->initialMaxBuildSize = initialMaxBuildSize; bwtInc->initialMaxBuildSize = initialMaxBuildSize;
bwtInc->incMaxBuildSize = incMaxBuildSize; bwtInc->incMaxBuildSize = incMaxBuildSize;
bwtInc->cumulativeCountInCurrentBuild = (bgint_t*)calloc((ALPHABET_SIZE + 1), sizeof(bgint_t)); bwtInc->cumulativeCountInCurrentBuild = (bgint_t*)xcalloc((ALPHABET_SIZE + 1), sizeof(bgint_t));
initializeVAL_bg(bwtInc->cumulativeCountInCurrentBuild, ALPHABET_SIZE + 1, 0); initializeVAL_bg(bwtInc->cumulativeCountInCurrentBuild, ALPHABET_SIZE + 1, 0);
// Build frequently accessed data // Build frequently accessed data
bwtInc->packedShift = (unsigned*)calloc(CHAR_PER_WORD, sizeof(unsigned int)); bwtInc->packedShift = (unsigned*)xcalloc(CHAR_PER_WORD, sizeof(unsigned int));
for (i=0; i<CHAR_PER_WORD; i++) for (i=0; i<CHAR_PER_WORD; i++)
bwtInc->packedShift[i] = BITS_IN_WORD - (i+1) * BIT_PER_CHAR; bwtInc->packedShift[i] = BITS_IN_WORD - (i+1) * BIT_PER_CHAR;
@ -372,7 +373,7 @@ BWTInc *BWTIncCreate(const bgint_t textLength, unsigned int initialMaxBuildSize,
+ incMaxBuildSize/5 * 3 * (sizeof(bgint_t) / 4); // space for the 3 temporary arrays in each iteration + incMaxBuildSize/5 * 3 * (sizeof(bgint_t) / 4); // space for the 3 temporary arrays in each iteration
if (bwtInc->availableWord < MIN_AVAILABLE_WORD) bwtInc->availableWord = MIN_AVAILABLE_WORD; // lh3: otherwise segfaul when availableWord is too small if (bwtInc->availableWord < MIN_AVAILABLE_WORD) bwtInc->availableWord = MIN_AVAILABLE_WORD; // lh3: otherwise segfaul when availableWord is too small
fprintf(stderr, "[%s] textLength=%ld, availableWord=%ld\n", __func__, (long)textLength, (long)bwtInc->availableWord); fprintf(stderr, "[%s] textLength=%ld, availableWord=%ld\n", __func__, (long)textLength, (long)bwtInc->availableWord);
bwtInc->workingMemory = (unsigned*)calloc(bwtInc->availableWord, BYTES_IN_WORD); bwtInc->workingMemory = (unsigned*)xcalloc(bwtInc->availableWord, BYTES_IN_WORD);
return bwtInc; return bwtInc;
} }
@ -1447,9 +1448,9 @@ BWTInc *BWTIncConstructFromPacked(const char *inputFileName, bgint_t initialMaxB
exit(1); exit(1);
} }
fseek(packedFile, -1, SEEK_END); err_fseek(packedFile, -1, SEEK_END);
packedFileLen = ftell(packedFile); packedFileLen = ftell(packedFile);
fread(&lastByteLength, sizeof(unsigned char), 1, packedFile); err_fread_noeof(&lastByteLength, sizeof(unsigned char), 1, packedFile);
totalTextLength = TextLengthFromBytePacked(packedFileLen, BIT_PER_CHAR, lastByteLength); totalTextLength = TextLengthFromBytePacked(packedFileLen, BIT_PER_CHAR, lastByteLength);
bwtInc = BWTIncCreate(totalTextLength, initialMaxBuildSize, incMaxBuildSize); bwtInc = BWTIncCreate(totalTextLength, initialMaxBuildSize, incMaxBuildSize);
@ -1463,10 +1464,10 @@ BWTInc *BWTIncConstructFromPacked(const char *inputFileName, bgint_t initialMaxB
} }
textSizeInByte = textToLoad / CHAR_PER_BYTE; // excluded the odd byte textSizeInByte = textToLoad / CHAR_PER_BYTE; // excluded the odd byte
fseek(packedFile, -2, SEEK_CUR); err_fseek(packedFile, -2, SEEK_CUR);
fseek(packedFile, -((long)textSizeInByte), SEEK_CUR); err_fseek(packedFile, -((long)textSizeInByte), SEEK_CUR);
fread(bwtInc->textBuffer, sizeof(unsigned char), textSizeInByte + 1, packedFile); err_fread_noeof(bwtInc->textBuffer, sizeof(unsigned char), textSizeInByte + 1, packedFile);
fseek(packedFile, -((long)textSizeInByte + 1), SEEK_CUR); err_fseek(packedFile, -((long)textSizeInByte + 1), SEEK_CUR);
ConvertBytePackedToWordPacked(bwtInc->textBuffer, bwtInc->packedText, ALPHABET_SIZE, textToLoad); ConvertBytePackedToWordPacked(bwtInc->textBuffer, bwtInc->packedText, ALPHABET_SIZE, textToLoad);
BWTIncConstruct(bwtInc, textToLoad); BWTIncConstruct(bwtInc, textToLoad);
@ -1479,9 +1480,9 @@ BWTInc *BWTIncConstructFromPacked(const char *inputFileName, bgint_t initialMaxB
textToLoad = totalTextLength - processedTextLength; textToLoad = totalTextLength - processedTextLength;
} }
textSizeInByte = textToLoad / CHAR_PER_BYTE; textSizeInByte = textToLoad / CHAR_PER_BYTE;
fseek(packedFile, -((long)textSizeInByte), SEEK_CUR); err_fseek(packedFile, -((long)textSizeInByte), SEEK_CUR);
fread(bwtInc->textBuffer, sizeof(unsigned char), textSizeInByte, packedFile); err_fread_noeof(bwtInc->textBuffer, sizeof(unsigned char), textSizeInByte, packedFile);
fseek(packedFile, -((long)textSizeInByte), SEEK_CUR); err_fseek(packedFile, -((long)textSizeInByte), SEEK_CUR);
ConvertBytePackedToWordPacked(bwtInc->textBuffer, bwtInc->packedText, ALPHABET_SIZE, textToLoad); ConvertBytePackedToWordPacked(bwtInc->textBuffer, bwtInc->packedText, ALPHABET_SIZE, textToLoad);
BWTIncConstruct(bwtInc, textToLoad); BWTIncConstruct(bwtInc, textToLoad);
processedTextLength += textToLoad; processedTextLength += textToLoad;
@ -1530,11 +1531,11 @@ void BWTSaveBwtCodeAndOcc(const BWT *bwt, const char *bwtFileName, const char *o
exit(1); exit(1);
} }
fwrite(&bwt->inverseSa0, sizeof(bgint_t), 1, bwtFile); err_fwrite(&bwt->inverseSa0, sizeof(bgint_t), 1, bwtFile);
fwrite(bwt->cumulativeFreq + 1, sizeof(bgint_t), ALPHABET_SIZE, bwtFile); err_fwrite(bwt->cumulativeFreq + 1, sizeof(bgint_t), ALPHABET_SIZE, bwtFile);
bwtLength = BWTFileSizeInWord(bwt->textLength); bwtLength = BWTFileSizeInWord(bwt->textLength);
fwrite(bwt->bwtCode, sizeof(unsigned int), bwtLength, bwtFile); err_fwrite(bwt->bwtCode, sizeof(unsigned int), bwtLength, bwtFile);
fclose(bwtFile); err_fclose(bwtFile);
} }
void bwt_bwtgen(const char *fn_pac, const char *fn_bwt) void bwt_bwtgen(const char *fn_pac, const char *fn_bwt)

11
bwt_lite.c
View File

@ -2,6 +2,7 @@
#include <string.h> #include <string.h>
#include <stdio.h> #include <stdio.h>
#include "bwt_lite.h" #include "bwt_lite.h"
#include "utils.h"
int is_sa(const uint8_t *T, uint32_t *SA, int n); int is_sa(const uint8_t *T, uint32_t *SA, int n);
int is_bwt(uint8_t *T, int n); int is_bwt(uint8_t *T, int n);
@ -10,21 +11,21 @@ bwtl_t *bwtl_seq2bwtl(int len, const uint8_t *seq)
{ {
bwtl_t *b; bwtl_t *b;
int i; int i;
b = (bwtl_t*)calloc(1, sizeof(bwtl_t)); b = (bwtl_t*)xcalloc(1, sizeof(bwtl_t));
b->seq_len = len; b->seq_len = len;
{ // calculate b->bwt { // calculate b->bwt
uint8_t *s; uint8_t *s;
b->sa = (uint32_t*)calloc(len + 1, 4); b->sa = (uint32_t*)xcalloc(len + 1, 4);
is_sa(seq, b->sa, len); is_sa(seq, b->sa, len);
s = (uint8_t*)calloc(len + 1, 1); s = (uint8_t*)xcalloc(len + 1, 1);
for (i = 0; i <= len; ++i) { for (i = 0; i <= len; ++i) {
if (b->sa[i] == 0) b->primary = i; if (b->sa[i] == 0) b->primary = i;
else s[i] = seq[b->sa[i] - 1]; else s[i] = seq[b->sa[i] - 1];
} }
for (i = b->primary; i < len; ++i) s[i] = s[i + 1]; for (i = b->primary; i < len; ++i) s[i] = s[i + 1];
b->bwt_size = (len + 15) / 16; b->bwt_size = (len + 15) / 16;
b->bwt = (uint32_t*)calloc(b->bwt_size, 4); b->bwt = (uint32_t*)xcalloc(b->bwt_size, 4);
for (i = 0; i < len; ++i) for (i = 0; i < len; ++i)
b->bwt[i>>4] |= s[i] << ((15 - (i&15)) << 1); b->bwt[i>>4] |= s[i] << ((15 - (i&15)) << 1);
free(s); free(s);
@ -32,7 +33,7 @@ bwtl_t *bwtl_seq2bwtl(int len, const uint8_t *seq)
{ // calculate b->occ { // calculate b->occ
uint32_t c[4]; uint32_t c[4];
b->n_occ = (len + 15) / 16 * 4; b->n_occ = (len + 15) / 16 * 4;
b->occ = (uint32_t*)calloc(b->n_occ, 4); b->occ = (uint32_t*)xcalloc(b->n_occ, 4);
memset(c, 0, 16); memset(c, 0, 16);
for (i = 0; i < len; ++i) { for (i = 0; i < len; ++i) {
if (i % 16 == 0) if (i % 16 == 0)

14
bwtaln.c
View File

@ -19,7 +19,7 @@
gap_opt_t *gap_init_opt() gap_opt_t *gap_init_opt()
{ {
gap_opt_t *o; gap_opt_t *o;
o = (gap_opt_t*)calloc(1, sizeof(gap_opt_t)); o = (gap_opt_t*)xcalloc(1, sizeof(gap_opt_t));
/* IMPORTANT: s_mm*10 should be about the average base error /* IMPORTANT: s_mm*10 should be about the average base error
rate. Voilating this requirement will break pairing! */ rate. Voilating this requirement will break pairing! */
o->s_mm = 3; o->s_gapo = 11; o->s_gape = 4; o->s_mm = 3; o->s_gapo = 11; o->s_gape = 4;
@ -89,7 +89,7 @@ void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt, int n_seqs, bwa_seq_t *seqs,
if (local_opt.max_diff < local_opt.max_gapo) local_opt.max_gapo = local_opt.max_diff; if (local_opt.max_diff < local_opt.max_gapo) local_opt.max_gapo = local_opt.max_diff;
stack = gap_init_stack(local_opt.max_diff, local_opt.max_gapo, local_opt.max_gape, &local_opt); stack = gap_init_stack(local_opt.max_diff, local_opt.max_gapo, local_opt.max_gape, &local_opt);
seed_w = (bwt_width_t*)calloc(opt->seed_len+1, sizeof(bwt_width_t)); seed_w = (bwt_width_t*)xcalloc(opt->seed_len+1, sizeof(bwt_width_t));
w = 0; w = 0;
for (i = 0; i != n_seqs; ++i) { for (i = 0; i != n_seqs; ++i) {
bwa_seq_t *p = seqs + i; bwa_seq_t *p = seqs + i;
@ -99,7 +99,7 @@ void bwa_cal_sa_reg_gap(int tid, bwt_t *const bwt, int n_seqs, bwa_seq_t *seqs,
p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0; p->sa = 0; p->type = BWA_TYPE_NO_MATCH; p->c1 = p->c2 = 0; p->n_aln = 0; p->aln = 0;
if (max_l < p->len) { if (max_l < p->len) {
max_l = p->len; max_l = p->len;
w = (bwt_width_t*)realloc(w, (max_l + 1) * sizeof(bwt_width_t)); w = (bwt_width_t*)xrealloc(w, (max_l + 1) * sizeof(bwt_width_t));
memset(w, 0, (max_l + 1) * sizeof(bwt_width_t)); memset(w, 0, (max_l + 1) * sizeof(bwt_width_t));
} }
bwt_cal_width(bwt, p->len, p->seq, w); bwt_cal_width(bwt, p->len, p->seq, w);
@ -162,7 +162,7 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
ks = bwa_open_reads(opt->mode, fn_fa); ks = bwa_open_reads(opt->mode, fn_fa);
{ // load BWT { // load BWT
char *str = (char*)calloc(strlen(prefix) + 10, 1); char *str = (char*)xcalloc(strlen(prefix) + 10, 1);
strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str); strcpy(str, prefix); strcat(str, ".bwt"); bwt = bwt_restore_bwt(str);
free(str); free(str);
} }
@ -185,8 +185,8 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
int j; int j;
pthread_attr_init(&attr); pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t)); data = (thread_aux_t*)xcalloc(opt->n_threads, sizeof(thread_aux_t));
tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t));
for (j = 0; j < opt->n_threads; ++j) { for (j = 0; j < opt->n_threads; ++j) {
data[j].tid = j; data[j].bwt = bwt; data[j].tid = j; data[j].bwt = bwt;
data[j].n_seqs = n_seqs; data[j].seqs = seqs; data[j].opt = opt; data[j].n_seqs = n_seqs; data[j].seqs = seqs; data[j].opt = opt;
@ -225,7 +225,7 @@ char *bwa_infer_prefix(const char *hint)
int l_hint; int l_hint;
FILE *fp; FILE *fp;
l_hint = strlen(hint); l_hint = strlen(hint);
prefix = malloc(l_hint + 3 + 4 + 1); prefix = xmalloc(l_hint + 3 + 4 + 1);
strcpy(prefix, hint); strcpy(prefix, hint);
strcpy(prefix + l_hint, ".64.bwt"); strcpy(prefix + l_hint, ".64.bwt");
if ((fp = fopen(prefix, "rb")) != 0) { if ((fp = fopen(prefix, "rb")) != 0) {

11
bwtgap.c
View File

@ -3,6 +3,7 @@
#include <string.h> #include <string.h>
#include "bwtgap.h" #include "bwtgap.h"
#include "bwtaln.h" #include "bwtaln.h"
#include "utils.h"
#define STATE_M 0 #define STATE_M 0
#define STATE_I 1 #define STATE_I 1
@ -13,9 +14,9 @@
gap_stack_t *gap_init_stack2(int max_score) gap_stack_t *gap_init_stack2(int max_score)
{ {
gap_stack_t *stack; gap_stack_t *stack;
stack = (gap_stack_t*)calloc(1, sizeof(gap_stack_t)); stack = (gap_stack_t*)xcalloc(1, sizeof(gap_stack_t));
stack->n_stacks = max_score; stack->n_stacks = max_score;
stack->stacks = (gap_stack1_t*)calloc(stack->n_stacks, sizeof(gap_stack1_t)); stack->stacks = (gap_stack1_t*)xcalloc(stack->n_stacks, sizeof(gap_stack1_t));
return stack; return stack;
} }
@ -51,7 +52,7 @@ static inline void gap_push(gap_stack_t *stack, int i, bwtint_t k, bwtint_t l, i
q = stack->stacks + score; q = stack->stacks + score;
if (q->n_entries == q->m_entries) { if (q->n_entries == q->m_entries) {
q->m_entries = q->m_entries? q->m_entries<<1 : 4; q->m_entries = q->m_entries? q->m_entries<<1 : 4;
q->stack = (gap_entry_t*)realloc(q->stack, sizeof(gap_entry_t) * q->m_entries); q->stack = (gap_entry_t*)xrealloc(q->stack, sizeof(gap_entry_t) * q->m_entries);
} }
p = q->stack + q->n_entries; p = q->stack + q->n_entries;
p->info = (u_int32_t)score<<21 | i; p->k = k; p->l = l; p->info = (u_int32_t)score<<21 | i; p->k = k; p->l = l;
@ -110,7 +111,7 @@ bwt_aln1_t *bwt_match_gap(bwt_t *const bwt, int len, const ubyte_t *seq, bwt_wid
bwt_aln1_t *aln; bwt_aln1_t *aln;
m_aln = 4; n_aln = 0; m_aln = 4; n_aln = 0;
aln = (bwt_aln1_t*)calloc(m_aln, sizeof(bwt_aln1_t)); aln = (bwt_aln1_t*)xcalloc(m_aln, sizeof(bwt_aln1_t));
// check whether there are too many N // check whether there are too many N
for (j = _j = 0; j < len; ++j) for (j = _j = 0; j < len; ++j)
@ -177,7 +178,7 @@ bwt_aln1_t *bwt_match_gap(bwt_t *const bwt, int len, const ubyte_t *seq, bwt_wid
gap_shadow(l - k + 1, len, bwt->seq_len, e.last_diff_pos, width); gap_shadow(l - k + 1, len, bwt->seq_len, e.last_diff_pos, width);
if (n_aln == m_aln) { if (n_aln == m_aln) {
m_aln <<= 1; m_aln <<= 1;
aln = (bwt_aln1_t*)realloc(aln, m_aln * sizeof(bwt_aln1_t)); aln = (bwt_aln1_t*)xrealloc(aln, m_aln * sizeof(bwt_aln1_t));
memset(aln + m_aln/2, 0, m_aln/2*sizeof(bwt_aln1_t)); memset(aln + m_aln/2, 0, m_aln/2*sizeof(bwt_aln1_t));
} }
p = aln + n_aln; p = aln + n_aln;

10
bwtindex.c
View File

@ -54,7 +54,7 @@ int bwa_index(int argc, char *argv[])
else if (strcmp(optarg, "is") == 0) algo_type = 3; else if (strcmp(optarg, "is") == 0) algo_type = 3;
else err_fatal(__func__, "unknown algorithm: '%s'.", optarg); else err_fatal(__func__, "unknown algorithm: '%s'.", optarg);
break; break;
case 'p': prefix = strdup(optarg); break; case 'p': prefix = xstrdup(optarg); break;
case 'c': is_color = 1; break; case 'c': is_color = 1; break;
case '6': is_64 = 1; break; case '6': is_64 = 1; break;
default: return 1; default: return 1;
@ -75,13 +75,13 @@ int bwa_index(int argc, char *argv[])
return 1; return 1;
} }
if (prefix == 0) { if (prefix == 0) {
prefix = malloc(strlen(argv[optind]) + 4); prefix = xmalloc(strlen(argv[optind]) + 4);
strcpy(prefix, argv[optind]); strcpy(prefix, argv[optind]);
if (is_64) strcat(prefix, ".64"); if (is_64) strcat(prefix, ".64");
} }
str = (char*)calloc(strlen(prefix) + 10, 1); str = (char*)xcalloc(strlen(prefix) + 10, 1);
str2 = (char*)calloc(strlen(prefix) + 10, 1); str2 = (char*)xcalloc(strlen(prefix) + 10, 1);
str3 = (char*)calloc(strlen(prefix) + 10, 1); str3 = (char*)xcalloc(strlen(prefix) + 10, 1);
if (is_color == 0) { // nucleotide indexing if (is_color == 0) { // nucleotide indexing
gzFile fp = xzopen(argv[optind], "r"); gzFile fp = xzopen(argv[optind], "r");

54
bwtio.c
View File

@ -6,24 +6,24 @@
void bwt_dump_bwt(const char *fn, const bwt_t *bwt) void bwt_dump_bwt(const char *fn, const bwt_t *bwt)
{ {
FILE *fp; FILE *fp = NULL;
fp = xopen(fn, "wb"); fp = xopen(fn, "wb");
fwrite(&bwt->primary, sizeof(bwtint_t), 1, fp); err_fwrite(&bwt->primary, sizeof(bwtint_t), 1, fp);
fwrite(bwt->L2+1, sizeof(bwtint_t), 4, fp); err_fwrite(bwt->L2+1, sizeof(bwtint_t), 4, fp);
fwrite(bwt->bwt, 4, bwt->bwt_size, fp); err_fwrite(bwt->bwt, 4, bwt->bwt_size, fp);
fclose(fp); err_fclose(fp);
} }
void bwt_dump_sa(const char *fn, const bwt_t *bwt) void bwt_dump_sa(const char *fn, const bwt_t *bwt)
{ {
FILE *fp; FILE *fp;
fp = xopen(fn, "wb"); fp = xopen(fn, "wb");
fwrite(&bwt->primary, sizeof(bwtint_t), 1, fp); err_fwrite(&bwt->primary, sizeof(bwtint_t), 1, fp);
fwrite(bwt->L2+1, sizeof(bwtint_t), 4, fp); err_fwrite(bwt->L2+1, sizeof(bwtint_t), 4, fp);
fwrite(&bwt->sa_intv, sizeof(bwtint_t), 1, fp); err_fwrite(&bwt->sa_intv, sizeof(bwtint_t), 1, fp);
fwrite(&bwt->seq_len, sizeof(bwtint_t), 1, fp); err_fwrite(&bwt->seq_len, sizeof(bwtint_t), 1, fp);
fwrite(bwt->sa + 1, sizeof(bwtint_t), bwt->n_sa - 1, fp); err_fwrite(bwt->sa + 1, sizeof(bwtint_t), bwt->n_sa - 1, fp);
fclose(fp); err_fclose(fp);
} }
void bwt_restore_sa(const char *fn, bwt_t *bwt) void bwt_restore_sa(const char *fn, bwt_t *bwt)
@ -33,19 +33,19 @@ void bwt_restore_sa(const char *fn, bwt_t *bwt)
bwtint_t primary; bwtint_t primary;
fp = xopen(fn, "rb"); fp = xopen(fn, "rb");
fread(&primary, sizeof(bwtint_t), 1, fp); err_fread_noeof(&primary, sizeof(bwtint_t), 1, fp);
xassert(primary == bwt->primary, "SA-BWT inconsistency: primary is not the same."); xassert(primary == bwt->primary, "SA-BWT inconsistency: primary is not the same.");
fread(skipped, sizeof(bwtint_t), 4, fp); // skip err_fread_noeof(skipped, sizeof(bwtint_t), 4, fp); // skip
fread(&bwt->sa_intv, sizeof(bwtint_t), 1, fp); err_fread_noeof(&bwt->sa_intv, sizeof(bwtint_t), 1, fp);
fread(&primary, sizeof(bwtint_t), 1, fp); err_fread_noeof(&primary, sizeof(bwtint_t), 1, fp);
xassert(primary == bwt->seq_len, "SA-BWT inconsistency: seq_len is not the same."); xassert(primary == bwt->seq_len, "SA-BWT inconsistency: seq_len is not the same.");
bwt->n_sa = (bwt->seq_len + bwt->sa_intv) / bwt->sa_intv; bwt->n_sa = (bwt->seq_len + bwt->sa_intv) / bwt->sa_intv;
bwt->sa = (bwtint_t*)calloc(bwt->n_sa, sizeof(bwtint_t)); bwt->sa = (bwtint_t*)xcalloc(bwt->n_sa, sizeof(bwtint_t));
bwt->sa[0] = -1; bwt->sa[0] = -1;
fread(bwt->sa + 1, sizeof(bwtint_t), bwt->n_sa - 1, fp); err_fread_noeof(bwt->sa + 1, sizeof(bwtint_t), bwt->n_sa - 1, fp);
fclose(fp); err_fclose(fp);
} }
bwt_t *bwt_restore_bwt(const char *fn) bwt_t *bwt_restore_bwt(const char *fn)
@ -53,17 +53,17 @@ bwt_t *bwt_restore_bwt(const char *fn)
bwt_t *bwt; bwt_t *bwt;
FILE *fp; FILE *fp;
bwt = (bwt_t*)calloc(1, sizeof(bwt_t)); bwt = (bwt_t*)xcalloc(1, sizeof(bwt_t));
fp = xopen(fn, "rb"); fp = xopen(fn, "rb");
fseek(fp, 0, SEEK_END); err_fseek(fp, 0, SEEK_END);
bwt->bwt_size = (ftell(fp) - sizeof(bwtint_t) * 5) >> 2; bwt->bwt_size = (err_ftell(fp) - sizeof(bwtint_t) * 5) >> 2;
bwt->bwt = (uint32_t*)calloc(bwt->bwt_size, 4); bwt->bwt = (uint32_t*)xcalloc(bwt->bwt_size, 4);
fseek(fp, 0, SEEK_SET); err_fseek(fp, 0, SEEK_SET);
fread(&bwt->primary, sizeof(bwtint_t), 1, fp); err_fread_noeof(&bwt->primary, sizeof(bwtint_t), 1, fp);
fread(bwt->L2+1, sizeof(bwtint_t), 4, fp); err_fread_noeof(bwt->L2+1, sizeof(bwtint_t), 4, fp);
fread(bwt->bwt, 4, bwt->bwt_size, fp); err_fread_noeof(bwt->bwt, 4, bwt->bwt_size, fp);
bwt->seq_len = bwt->L2[4]; bwt->seq_len = bwt->L2[4];
fclose(fp); err_fclose(fp);
bwt_gen_cnt_table(bwt); bwt_gen_cnt_table(bwt);
return bwt; return bwt;

38
bwtmisc.c
View File

@ -46,10 +46,10 @@ int64_t bwa_seq_len(const char *fn_pac)
int64_t pac_len; int64_t pac_len;
ubyte_t c; ubyte_t c;
fp = xopen(fn_pac, "rb"); fp = xopen(fn_pac, "rb");
fseek(fp, -1, SEEK_END); err_fseek(fp, -1, SEEK_END);
pac_len = ftell(fp); pac_len = err_ftell(fp);
fread(&c, 1, 1, fp); err_fread_noeof(&c, 1, 1, fp);
fclose(fp); err_fclose(fp);
return (pac_len - 1) * 4 + (int)c; return (pac_len - 1) * 4 + (int)c;
} }
@ -61,18 +61,18 @@ bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is)
FILE *fp; FILE *fp;
// initialization // initialization
bwt = (bwt_t*)calloc(1, sizeof(bwt_t)); bwt = (bwt_t*)xcalloc(1, sizeof(bwt_t));
bwt->seq_len = bwa_seq_len(fn_pac); bwt->seq_len = bwa_seq_len(fn_pac);
bwt->bwt_size = (bwt->seq_len + 15) >> 4; bwt->bwt_size = (bwt->seq_len + 15) >> 4;
fp = xopen(fn_pac, "rb"); fp = xopen(fn_pac, "rb");
// prepare sequence // prepare sequence
pac_size = (bwt->seq_len>>2) + ((bwt->seq_len&3) == 0? 0 : 1); pac_size = (bwt->seq_len>>2) + ((bwt->seq_len&3) == 0? 0 : 1);
buf2 = (ubyte_t*)calloc(pac_size, 1); buf2 = (ubyte_t*)xcalloc(pac_size, 1);
fread(buf2, 1, pac_size, fp); err_fread_noeof(buf2, 1, pac_size, fp);
fclose(fp); err_fclose(fp);
memset(bwt->L2, 0, 5 * 4); memset(bwt->L2, 0, 5 * 4);
buf = (ubyte_t*)calloc(bwt->seq_len + 1, 1); buf = (ubyte_t*)xcalloc(bwt->seq_len + 1, 1);
for (i = 0; i < bwt->seq_len; ++i) { for (i = 0; i < bwt->seq_len; ++i) {
buf[i] = buf2[i>>2] >> ((3 - (i&3)) << 1) & 3; buf[i] = buf2[i>>2] >> ((3 - (i&3)) << 1) & 3;
++bwt->L2[1+buf[i]]; ++bwt->L2[1+buf[i]];
@ -90,7 +90,7 @@ bwt_t *bwt_pac2bwt(const char *fn_pac, int use_is)
err_fatal_simple("libdivsufsort is not compiled in."); err_fatal_simple("libdivsufsort is not compiled in.");
#endif #endif
} }
bwt->bwt = (u_int32_t*)calloc(bwt->bwt_size, 4); bwt->bwt = (u_int32_t*)xcalloc(bwt->bwt_size, 4);
for (i = 0; i < bwt->seq_len; ++i) for (i = 0; i < bwt->seq_len; ++i)
bwt->bwt[i>>4] |= buf[i] << ((15 - (i&15)) << 1); bwt->bwt[i>>4] |= buf[i] << ((15 - (i&15)) << 1);
free(buf); free(buf);
@ -126,7 +126,7 @@ void bwt_bwtupdate_core(bwt_t *bwt)
n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1;
bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size
buf = (uint32_t*)calloc(bwt->bwt_size, 4); // will be the new bwt buf = (uint32_t*)xcalloc(bwt->bwt_size, 4); // will be the new bwt
c[0] = c[1] = c[2] = c[3] = 0; c[0] = c[1] = c[2] = c[3] = 0;
for (i = k = 0; i < bwt->seq_len; ++i) { for (i = k = 0; i < bwt->seq_len; ++i) {
if (i % OCC_INTERVAL == 0) { if (i % OCC_INTERVAL == 0) {
@ -167,10 +167,10 @@ uint8_t *bwa_pac2cspac_core(const bntseq_t *bns)
uint8_t *pac, *cspac; uint8_t *pac, *cspac;
bwtint_t i; bwtint_t i;
int c1, c2; int c1, c2;
pac = (uint8_t*)calloc(bns->l_pac/4 + 1, 1); pac = (uint8_t*)xcalloc(bns->l_pac/4 + 1, 1);
cspac = (uint8_t*)calloc(bns->l_pac/4 + 1, 1); cspac = (uint8_t*)xcalloc(bns->l_pac/4 + 1, 1);
fread(pac, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
rewind(bns->fp_pac); err_rewind(bns->fp_pac);
c1 = pac[0]>>6; cspac[0] = c1<<6; c1 = pac[0]>>6; cspac[0] = c1<<6;
for (i = 1; i < bns->l_pac; ++i) { for (i = 1; i < bns->l_pac; ++i) {
c2 = pac[i>>2] >> (~i&3)*2 & 3; c2 = pac[i>>2] >> (~i&3)*2 & 3;
@ -196,13 +196,13 @@ int bwa_pac2cspac(int argc, char *argv[])
cspac = bwa_pac2cspac_core(bns); cspac = bwa_pac2cspac_core(bns);
bns_dump(bns, argv[2]); bns_dump(bns, argv[2]);
// now write cspac // now write cspac
str = (char*)calloc(strlen(argv[2]) + 5, 1); str = (char*)xcalloc(strlen(argv[2]) + 5, 1);
strcat(strcpy(str, argv[2]), ".pac"); strcat(strcpy(str, argv[2]), ".pac");
fp = xopen(str, "wb"); fp = xopen(str, "wb");
fwrite(cspac, 1, bns->l_pac/4 + 1, fp); err_fwrite(cspac, 1, bns->l_pac/4 + 1, fp);
ct = bns->l_pac % 4; ct = bns->l_pac % 4;
fwrite(&ct, 1, 1, fp); err_fwrite(&ct, 1, 1, fp);
fclose(fp); err_fclose(fp);
bns_destroy(bns); bns_destroy(bns);
free(cspac); free(cspac);
return 0; return 0;

74
bwtsw2_aux.c
View File

@ -15,7 +15,7 @@
#include "kstring.h" #include "kstring.h"
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) KSEQ_INIT(gzFile, err_gzread)
#include "ksort.h" #include "ksort.h"
#define __left_lt(a, b) ((a).end > (b).end) #define __left_lt(a, b) ((a).end > (b).end)
@ -47,7 +47,7 @@ extern int bsw2_resolve_query_overlaps(bwtsw2_t *b, float mask_level);
bsw2opt_t *bsw2_init_opt() bsw2opt_t *bsw2_init_opt()
{ {
bsw2opt_t *o = (bsw2opt_t*)calloc(1, sizeof(bsw2opt_t)); bsw2opt_t *o = (bsw2opt_t*)xcalloc(1, sizeof(bsw2opt_t));
o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30; o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30;
o->bw = 50; o->bw = 50;
o->max_ins = 20000; o->max_ins = 20000;
@ -72,11 +72,11 @@ void bsw2_destroy(bwtsw2_t *b)
bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b) bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b)
{ {
bwtsw2_t *p; bwtsw2_t *p;
p = calloc(1, sizeof(bwtsw2_t)); p = xcalloc(1, sizeof(bwtsw2_t));
p->max = p->n = b->n; p->max = p->n = b->n;
if (b->n) { if (b->n) {
kroundup32(p->max); kroundup32(p->max);
p->hits = calloc(p->max, sizeof(bsw2hit_t)); p->hits = xcalloc(p->max, sizeof(bsw2hit_t));
memcpy(p->hits, b->hits, p->n * sizeof(bsw2hit_t)); memcpy(p->hits, b->hits, p->n * sizeof(bsw2hit_t));
} }
return p; return p;
@ -100,10 +100,10 @@ void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq
par.matrix = matrix; par.matrix = matrix;
__gen_ap(par, opt); __gen_ap(par, opt);
query = calloc(lq, 1); query = xcalloc(lq, 1);
// sort according to the descending order of query end // sort according to the descending order of query end
ks_introsort(hit, b->n, b->hits); ks_introsort(hit, b->n, b->hits);
target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1); target = xcalloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
// reverse _query // reverse _query
for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i]; for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i];
// core loop // core loop
@ -146,7 +146,7 @@ void bsw2_extend_rght(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *query, int lq,
par.matrix = matrix; par.matrix = matrix;
__gen_ap(par, opt); __gen_ap(par, opt);
target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1); target = xcalloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
for (i = 0; i < b->n; ++i) { for (i = 0; i < b->n; ++i) {
bsw2hit_t *p = b->hits + i; bsw2hit_t *p = b->hits + i;
int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq; int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
@ -178,8 +178,8 @@ static void gen_cigar(const bsw2opt_t *opt, int lq, uint8_t *seq[2], const uint8
par.matrix = matrix; par.matrix = matrix;
__gen_ap(par, opt); __gen_ap(par, opt);
i = ((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq; // maximum possible target length i = ((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq; // maximum possible target length
target = calloc(i, 1); target = xcalloc(i, 1);
path = calloc(i + lq, sizeof(path_t)); path = xcalloc(i + lq, sizeof(path_t));
// generate CIGAR // generate CIGAR
for (i = 0; i < b->n; ++i) { for (i = 0; i < b->n; ++i) {
bsw2hit_t *p = b->hits + i; bsw2hit_t *p = b->hits + i;
@ -206,7 +206,7 @@ static void gen_cigar(const bsw2opt_t *opt, int lq, uint8_t *seq[2], const uint8
} }
#endif #endif
if (beg != 0 || end < lq) { // write soft clipping if (beg != 0 || end < lq) { // write soft clipping
q->cigar = realloc(q->cigar, 4 * (q->n_cigar + 2)); q->cigar = xrealloc(q->cigar, 4 * (q->n_cigar + 2));
if (beg != 0) { if (beg != 0) {
memmove(q->cigar + 1, q->cigar, q->n_cigar * 4); memmove(q->cigar + 1, q->cigar, q->n_cigar * 4);
q->cigar[0] = beg<<4 | 4; q->cigar[0] = beg<<4 | 4;
@ -238,7 +238,7 @@ static void merge_hits(bwtsw2_t *b[2], int l, int is_reverse)
int i; int i;
if (b[0]->n + b[1]->n > b[0]->max) { if (b[0]->n + b[1]->n > b[0]->max) {
b[0]->max = b[0]->n + b[1]->n; b[0]->max = b[0]->n + b[1]->n;
b[0]->hits = realloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t)); b[0]->hits = xrealloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t));
} }
for (i = 0; i < b[1]->n; ++i) { for (i = 0; i < b[1]->n; ++i) {
bsw2hit_t *p = b[0]->hits + b[0]->n + i; bsw2hit_t *p = b[0]->hits + b[0]->n + i;
@ -266,9 +266,9 @@ static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8
_b = bsw2_core(bns, opt, query, target, pool); _b = bsw2_core(bns, opt, query, target, pool);
bwtl_destroy(query); bwtl_destroy(query);
for (k = 0; k < 2; ++k) { for (k = 0; k < 2; ++k) {
bb[k] = calloc(2, sizeof(void*)); bb[k] = xcalloc(2, sizeof(void*));
bb[k][0] = calloc(1, sizeof(bwtsw2_t)); bb[k][0] = xcalloc(1, sizeof(bwtsw2_t));
bb[k][1] = calloc(1, sizeof(bwtsw2_t)); bb[k][1] = xcalloc(1, sizeof(bwtsw2_t));
} }
for (k = 0; k < 2; ++k) { // separate _b into bb[2] based on the strand for (k = 0; k < 2; ++k) { // separate _b into bb[2] based on the strand
for (j = 0; j < _b[k]->n; ++j) { for (j = 0; j < _b[k]->n; ++j) {
@ -276,7 +276,7 @@ static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8
p = bb[_b[k]->hits[j].is_rev][k]; p = bb[_b[k]->hits[j].is_rev][k];
if (p->n == p->max) { if (p->n == p->max) {
p->max = p->max? p->max<<1 : 8; p->max = p->max? p->max<<1 : 8;
p->hits = realloc(p->hits, p->max * sizeof(bsw2hit_t)); p->hits = xrealloc(p->hits, p->max * sizeof(bsw2hit_t));
} }
q = &p->hits[p->n++]; q = &p->hits[p->n++];
*q = _b[k]->hits[j]; *q = _b[k]->hits[j];
@ -355,7 +355,7 @@ static int fix_cigar(const bntseq_t *bns, bsw2hit_t *p, int n_cigar, uint32_t *c
uint32_t *cn; uint32_t *cn;
bwtint_t kk = 0; bwtint_t kk = 0;
nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0; nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0;
cn = calloc(n_cigar + 3, 4); cn = xcalloc(n_cigar + 3, 4);
x = coor; y = 0; x = coor; y = 0;
for (i = j = 0; i < n_cigar; ++i) { for (i = j = 0; i < n_cigar; ++i) {
int op = cigar[i]&0xf, ln = cigar[i]>>4; int op = cigar[i]&0xf, ln = cigar[i]>>4;
@ -434,9 +434,9 @@ static void write_aux(const bsw2opt_t *opt, const bntseq_t *bns, int qlen, uint8
if (b->n<<1 < b->max) { if (b->n<<1 < b->max) {
b->max = b->n; b->max = b->n;
kroundup32(b->max); kroundup32(b->max);
b->hits = realloc(b->hits, b->max * sizeof(bsw2hit_t)); b->hits = xrealloc(b->hits, b->max * sizeof(bsw2hit_t));
} }
b->aux = calloc(b->n, sizeof(bsw2aux_t)); b->aux = xcalloc(b->n, sizeof(bsw2aux_t));
// generate CIGAR // generate CIGAR
gen_cigar(opt, qlen, seq, pac, b, name); gen_cigar(opt, qlen, seq, pac, b, name);
// fix CIGAR, generate mapQ, and write chromosomal position // fix CIGAR, generate mapQ, and write chromosomal position
@ -596,7 +596,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
bsw2opt_t opt; bsw2opt_t opt;
bsw2global_t *pool = bsw2_global_init(); bsw2global_t *pool = bsw2_global_init();
bwtsw2_t **buf; bwtsw2_t **buf;
buf = calloc(_seq->n, sizeof(void*)); buf = xcalloc(_seq->n, sizeof(void*));
for (x = 0; x < _seq->n; ++x) { for (x = 0; x < _seq->n; ++x) {
bsw2seq1_t *p = _seq->seq + x; bsw2seq1_t *p = _seq->seq + x;
uint8_t *seq[2], *rseq[2]; uint8_t *seq[2], *rseq[2];
@ -607,10 +607,10 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
if (pool->max_l < l) { // then enlarge working space for aln_extend_core() if (pool->max_l < l) { // then enlarge working space for aln_extend_core()
int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l; int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
pool->max_l = l; pool->max_l = l;
pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24); pool->aln_mem = xrealloc(pool->aln_mem, (tmp + 2) * 24);
} }
// set seq[2] and rseq[2] // set seq[2] and rseq[2]
seq[0] = calloc(l * 4, 1); seq[0] = xcalloc(l * 4, 1);
seq[1] = seq[0] + l; seq[1] = seq[0] + l;
rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l; rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
// convert sequences to 2-bit representation // convert sequences to 2-bit representation
@ -623,7 +623,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
rseq[1][i] = c; rseq[1][i] = c;
} }
if (l - k < opt.t) { // too few unambiguous bases if (l - k < opt.t) { // too few unambiguous bases
buf[x] = calloc(1, sizeof(bwtsw2_t)); buf[x] = xcalloc(1, sizeof(bwtsw2_t));
free(seq[0]); continue; free(seq[0]); continue;
} }
// alignment // alignment
@ -655,7 +655,7 @@ static void bsw2_aln_core(bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t
bsw2seq1_t *p = _seq->seq + x; bsw2seq1_t *p = _seq->seq + x;
uint8_t *seq[2]; uint8_t *seq[2];
int i; int i;
seq[0] = malloc(p->l * 2); seq[1] = seq[0] + p->l; seq[0] = xmalloc(p->l * 2); seq[1] = seq[0] + p->l;
for (i = 0; i < p->l; ++i) { for (i = 0; i < p->l; ++i) {
int c = nst_nt4_table[(int)p->seq[i]]; int c = nst_nt4_table[(int)p->seq[i]];
if (c >= 4) c = (int)(drand48() * 4); if (c >= 4) c = (int)(drand48() * 4);
@ -711,16 +711,16 @@ static void process_seqs(bsw2seq_t *_seq, const bsw2opt_t *opt, const bntseq_t *
int j; int j;
pthread_attr_init(&attr); pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE); pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t)); data = (thread_aux_t*)xcalloc(opt->n_threads, sizeof(thread_aux_t));
tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t)); tid = (pthread_t*)xcalloc(opt->n_threads, sizeof(pthread_t));
for (j = 0; j < opt->n_threads; ++j) { for (j = 0; j < opt->n_threads; ++j) {
thread_aux_t *p = data + j; thread_aux_t *p = data + j;
p->tid = j; p->_opt = opt; p->bns = bns; p->is_pe = is_pe; p->tid = j; p->_opt = opt; p->bns = bns; p->is_pe = is_pe;
p->pac = pac; p->target = target; p->pac = pac; p->target = target;
p->_seq = calloc(1, sizeof(bsw2seq_t)); p->_seq = xcalloc(1, sizeof(bsw2seq_t));
p->_seq->max = (_seq->n + opt->n_threads - 1) / opt->n_threads + 1; p->_seq->max = (_seq->n + opt->n_threads - 1) / opt->n_threads + 1;
p->_seq->n = 0; p->_seq->n = 0;
p->_seq->seq = calloc(p->_seq->max, sizeof(bsw2seq1_t)); p->_seq->seq = xcalloc(p->_seq->max, sizeof(bsw2seq1_t));
} }
for (i = 0; i < _seq->n; ++i) { // assign sequences to each thread for (i = 0; i < _seq->n; ++i) { // assign sequences to each thread
bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq; bsw2seq_t *p = data[(i>>is_pe)%opt->n_threads]._seq;
@ -760,10 +760,10 @@ static void kseq_to_bsw2seq(const kseq_t *ks, bsw2seq1_t *p)
{ {
p->tid = -1; p->tid = -1;
p->l = ks->seq.l; p->l = ks->seq.l;
p->name = strdup(ks->name.s); p->name = xstrdup(ks->name.s);
p->seq = strdup(ks->seq.s); p->seq = xstrdup(ks->seq.s);
p->qual = ks->qual.l? strdup(ks->qual.s) : 0; p->qual = ks->qual.l? xstrdup(ks->qual.s) : 0;
p->comment = ks->comment.l? strdup(ks->comment.s) : 0; p->comment = ks->comment.l? xstrdup(ks->comment.s) : 0;
p->sam = 0; p->sam = 0;
} }
@ -775,17 +775,13 @@ void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target, c
uint8_t *pac; uint8_t *pac;
bsw2seq_t *_seq; bsw2seq_t *_seq;
pac = calloc(bns->l_pac/4+1, 1); pac = xcalloc(bns->l_pac/4+1, 1);
if (pac == 0) {
fprintf(stderr, "[bsw2_aln] insufficient memory!\n");
return;
}
for (l = 0; l < bns->n_seqs; ++l) for (l = 0; l < bns->n_seqs; ++l)
printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len); printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len);
fread(pac, 1, bns->l_pac/4+1, bns->fp_pac); err_fread_noeof(pac, 1, bns->l_pac/4+1, bns->fp_pac);
fp = xzopen(fn, "r"); fp = xzopen(fn, "r");
ks = kseq_init(fp); ks = kseq_init(fp);
_seq = calloc(1, sizeof(bsw2seq_t)); _seq = xcalloc(1, sizeof(bsw2seq_t));
if (fn2) { if (fn2) {
fp2 = xzopen(fn2, "r"); fp2 = xzopen(fn2, "r");
ks2 = kseq_init(fp2); ks2 = kseq_init(fp2);
@ -796,7 +792,7 @@ void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target, c
ks->name.l -= 2, ks->name.s[ks->name.l] = 0; ks->name.l -= 2, ks->name.s[ks->name.l] = 0;
if (_seq->n == _seq->max) { if (_seq->n == _seq->max) {
_seq->max = _seq->max? _seq->max<<1 : 1024; _seq->max = _seq->max? _seq->max<<1 : 1024;
_seq->seq = realloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t)); _seq->seq = xrealloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t));
} }
kseq_to_bsw2seq(ks, &_seq->seq[_seq->n++]); kseq_to_bsw2seq(ks, &_seq->seq[_seq->n++]);
size += ks->seq.l; size += ks->seq.l;

7
bwtsw2_chain.c
View File

@ -1,5 +1,6 @@
#include <stdio.h> #include <stdio.h>
#include "bwtsw2.h" #include "bwtsw2.h"
#include "utils.h"
typedef struct { typedef struct {
uint32_t tbeg, tend; uint32_t tbeg, tend;
@ -48,9 +49,9 @@ void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2])
char *flag; char *flag;
// initialization // initialization
n[0] = b[0]->n; n[1] = b[1]->n; n[0] = b[0]->n; n[1] = b[1]->n;
z[0] = calloc(n[0] + n[1], sizeof(hsaip_t)); z[0] = xcalloc(n[0] + n[1], sizeof(hsaip_t));
z[1] = z[0] + n[0]; z[1] = z[0] + n[0];
chain[0] = calloc(n[0] + n[1], sizeof(hsaip_t)); chain[0] = xcalloc(n[0] + n[1], sizeof(hsaip_t));
for (k = j = 0; k < 2; ++k) { for (k = j = 0; k < 2; ++k) {
for (i = 0; i < b[k]->n; ++i) { for (i = 0; i < b[k]->n; ++i) {
bsw2hit_t *p = b[k]->hits + i; bsw2hit_t *p = b[k]->hits + i;
@ -73,7 +74,7 @@ void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2])
} }
//for (k = 0; k < m[0]; ++k) printf("%d, [%d,%d), [%d,%d)\n", chain[0][k].chain, chain[0][k].tbeg, chain[0][k].tend, chain[0][k].qbeg, chain[0][k].qend); //for (k = 0; k < m[0]; ++k) printf("%d, [%d,%d), [%d,%d)\n", chain[0][k].chain, chain[0][k].tbeg, chain[0][k].tend, chain[0][k].qbeg, chain[0][k].qend);
// filtering // filtering
flag = calloc(m[0] + m[1], 1); flag = xcalloc(m[0] + m[1], 1);
ks_introsort(hsaip, m[0] + m[1], chain[0]); ks_introsort(hsaip, m[0] + m[1], chain[0]);
for (k = 1; k < m[0] + m[1]; ++k) { for (k = 1; k < m[0] + m[1]; ++k) {
hsaip_t *p = chain[0] + k; hsaip_t *p = chain[0] + k;

29
bwtsw2_core.c
View File

@ -7,6 +7,7 @@
#include "bwtsw2.h" #include "bwtsw2.h"
#include "bwt.h" #include "bwt.h"
#include "kvec.h" #include "kvec.h"
#include "utils.h"
typedef struct { typedef struct {
bwtint_t k, l; bwtint_t k, l;
@ -71,7 +72,7 @@ typedef struct __mempool_t {
inline static bsw2entry_p mp_alloc(mempool_t *mp) inline static bsw2entry_p mp_alloc(mempool_t *mp)
{ {
++mp->cnt; ++mp->cnt;
if (kv_size(mp->pool) == 0) return (bsw2entry_t*)calloc(1, sizeof(bsw2entry_t)); if (kv_size(mp->pool) == 0) return (bsw2entry_t*)xcalloc(1, sizeof(bsw2entry_t));
else return kv_pop(mp->pool); else return kv_pop(mp->pool);
} }
inline static void mp_free(mempool_t *mp, bsw2entry_p e) inline static void mp_free(mempool_t *mp, bsw2entry_p e)
@ -133,7 +134,7 @@ static void cut_tail(bsw2entry_t *u, int T, bsw2entry_t *aux)
if (u->n <= T) return; if (u->n <= T) return;
if (aux->max < u->n) { if (aux->max < u->n) {
aux->max = u->n; aux->max = u->n;
aux->array = (bsw2cell_t*)realloc(aux->array, aux->max * sizeof(bsw2cell_t)); aux->array = (bsw2cell_t*)xrealloc(aux->array, aux->max * sizeof(bsw2cell_t));
} }
a = (int*)aux->array; a = (int*)aux->array;
for (i = n = 0; i != u->n; ++i) for (i = n = 0; i != u->n; ++i)
@ -184,7 +185,7 @@ static void merge_entry(const bsw2opt_t * __restrict opt, bsw2entry_t *u, bsw2en
int i; int i;
if (u->n + v->n >= u->max) { if (u->n + v->n >= u->max) {
u->max = u->n + v->n; u->max = u->n + v->n;
u->array = (bsw2cell_t*)realloc(u->array, u->max * sizeof(bsw2cell_t)); u->array = (bsw2cell_t*)xrealloc(u->array, u->max * sizeof(bsw2cell_t));
} }
for (i = 0; i != v->n; ++i) { for (i = 0; i != v->n; ++i) {
bsw2cell_t *p = v->array + i; bsw2cell_t *p = v->array + i;
@ -202,7 +203,7 @@ static inline bsw2cell_t *push_array_p(bsw2entry_t *e)
{ {
if (e->n == e->max) { if (e->n == e->max) {
e->max = e->max? e->max<<1 : 256; e->max = e->max? e->max<<1 : 256;
e->array = (bsw2cell_t*)realloc(e->array, sizeof(bsw2cell_t) * e->max); e->array = (bsw2cell_t*)xrealloc(e->array, sizeof(bsw2cell_t) * e->max);
} }
return e->array + e->n; return e->array + e->n;
} }
@ -250,7 +251,7 @@ static void save_narrow_hits(const bwtl_t *bwtl, bsw2entry_t *u, bwtsw2_t *b1, i
if (p->G >= t && p->ql - p->qk + 1 <= IS) { // good narrow hit if (p->G >= t && p->ql - p->qk + 1 <= IS) { // good narrow hit
if (b1->max == b1->n) { if (b1->max == b1->n) {
b1->max = b1->max? b1->max<<1 : 4; b1->max = b1->max? b1->max<<1 : 4;
b1->hits = realloc(b1->hits, b1->max * sizeof(bsw2hit_t)); b1->hits = xrealloc(b1->hits, b1->max * sizeof(bsw2hit_t));
} }
q = &b1->hits[b1->n++]; q = &b1->hits[b1->n++];
q->k = p->qk; q->l = p->ql; q->k = p->qk; q->l = p->ql;
@ -279,7 +280,7 @@ int bsw2_resolve_duphits(const bntseq_t *bns, const bwt_t *bwt, bwtsw2_t *b, int
else if (p->G > 0) ++n; else if (p->G > 0) ++n;
} }
b->n = b->max = n; b->n = b->max = n;
b->hits = calloc(b->max, sizeof(bsw2hit_t)); b->hits = xcalloc(b->max, sizeof(bsw2hit_t));
for (i = j = 0; i < old_n; ++i) { for (i = j = 0; i < old_n; ++i) {
bsw2hit_t *p = old_hits + i; bsw2hit_t *p = old_hits + i;
if (p->l - p->k + 1 <= IS) { // the hit is no so repetitive if (p->l - p->k + 1 <= IS) { // the hit is no so repetitive
@ -399,9 +400,9 @@ bsw2global_t *bsw2_global_init()
{ {
bsw2global_t *pool; bsw2global_t *pool;
bsw2stack_t *stack; bsw2stack_t *stack;
pool = calloc(1, sizeof(bsw2global_t)); pool = xcalloc(1, sizeof(bsw2global_t));
stack = calloc(1, sizeof(bsw2stack_t)); stack = xcalloc(1, sizeof(bsw2stack_t));
stack->pool = (mempool_t*)calloc(1, sizeof(mempool_t)); stack->pool = (mempool_t*)xcalloc(1, sizeof(mempool_t));
pool->stack = (void*)stack; pool->stack = (void*)stack;
return pool; return pool;
} }
@ -461,13 +462,13 @@ bwtsw2_t **bsw2_core(const bntseq_t *bns, const bsw2opt_t *opt, const bwtl_t *ta
rhash = kh_init(qintv); rhash = kh_init(qintv);
init_bwtsw2(target, query, stack); init_bwtsw2(target, query, stack);
heap_size = opt->z; heap_size = opt->z;
heap = calloc(heap_size, sizeof(int)); heap = xcalloc(heap_size, sizeof(int));
// initialize the return struct // initialize the return struct
b = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t)); b = (bwtsw2_t*)xcalloc(1, sizeof(bwtsw2_t));
b->n = b->max = target->seq_len * 2; b->n = b->max = target->seq_len * 2;
b->hits = calloc(b->max, sizeof(bsw2hit_t)); b->hits = xcalloc(b->max, sizeof(bsw2hit_t));
b1 = (bwtsw2_t*)calloc(1, sizeof(bwtsw2_t)); b1 = (bwtsw2_t*)xcalloc(1, sizeof(bwtsw2_t));
b_ret = calloc(2, sizeof(void*)); b_ret = xcalloc(2, sizeof(void*));
b_ret[0] = b; b_ret[1] = b1; b_ret[0] = b; b_ret[1] = b1;
// initialize timer // initialize timer
getrusage(0, &last); getrusage(0, &last);

7
bwtsw2_pair.c
View File

@ -2,6 +2,7 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "utils.h"
#include "bwt.h" #include "bwt.h"
#include "bntseq.h" #include "bntseq.h"
#include "bwtsw2.h" #include "bwtsw2.h"
@ -30,7 +31,7 @@ bsw2pestat_t bsw2_stat(int n, bwtsw2_t **buf, kstring_t *msg, int max_ins)
bsw2pestat_t r; bsw2pestat_t r;
memset(&r, 0, sizeof(bsw2pestat_t)); memset(&r, 0, sizeof(bsw2pestat_t));
isize = calloc(n, 8); isize = xcalloc(n, 8);
for (i = k = 0; i < n; i += 2) { for (i = k = 0; i < n; i += 2) {
bsw2hit_t *t[2]; bsw2hit_t *t[2];
int l; int l;
@ -113,7 +114,7 @@ void bsw2_pair1(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, const b
if (end > l_pac) end = l_pac; if (end > l_pac) end = l_pac;
if (end - beg < l_mseq) return; if (end - beg < l_mseq) return;
// generate the sequence // generate the sequence
seq = malloc(l_mseq + (end - beg)); seq = xmalloc(l_mseq + (end - beg));
ref = seq + l_mseq; ref = seq + l_mseq;
for (k = beg; k < end; ++k) for (k = beg; k < end; ++k)
ref[k - beg] = pac[k>>2] >> ((~k&3)<<1) & 0x3; ref[k - beg] = pac[k>>2] >> ((~k&3)<<1) & 0x3;
@ -221,7 +222,7 @@ void bsw2_pair(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, int n, b
a[which].flag |= BSW2_FLAG_RESCUED; a[which].flag |= BSW2_FLAG_RESCUED;
if (p[1]->max == 0) { if (p[1]->max == 0) {
p[1]->max = 1; p[1]->max = 1;
p[1]->hits = malloc(sizeof(bsw2hit_t)); p[1]->hits = xmalloc(sizeof(bsw2hit_t));
} }
p[1]->hits[0] = a[which]; p[1]->hits[0] = a[which];
p[1]->n = 1; p[1]->n = 1;

3
cs2nt.c
View File

@ -3,6 +3,7 @@
#include <stdlib.h> #include <stdlib.h>
#include "bwtaln.h" #include "bwtaln.h"
#include "stdaln.h" #include "stdaln.h"
#include "utils.h"
/* /*
Here is a delicate example. ref_nt=ATTAAC(RBRBG), read_cs=RBBOG. If we Here is a delicate example. ref_nt=ATTAAC(RBRBG), read_cs=RBBOG. If we
@ -118,7 +119,7 @@ void bwa_cs2nt_core(bwa_seq_t *p, bwtint_t l_pac, ubyte_t *pac)
// set temporary arrays // set temporary arrays
if (p->type == BWA_TYPE_NO_MATCH) return; if (p->type == BWA_TYPE_NO_MATCH) return;
len = p->len + p->n_gapo + p->n_gape + 100; // leave enough space len = p->len + p->n_gapo + p->n_gape + 100; // leave enough space
ta = (uint8_t*)malloc(len * 7); ta = (uint8_t*)xmalloc(len * 7);
nt_ref = ta; nt_ref = ta;
cs_read = nt_ref + len; cs_read = nt_ref + len;
nt_read = cs_read + len; nt_read = cs_read + len;

13
fastmap.c
View File

@ -6,7 +6,8 @@
#include "bwt.h" #include "bwt.h"
#include "kvec.h" #include "kvec.h"
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) #include "utils.h"
KSEQ_INIT(gzFile, err_gzread)
extern unsigned char nst_nt4_table[256]; extern unsigned char nst_nt4_table[256];
@ -20,11 +21,11 @@ typedef struct {
smem_i *smem_iter_init(const bwt_t *bwt) smem_i *smem_iter_init(const bwt_t *bwt)
{ {
smem_i *iter; smem_i *iter;
iter = calloc(1, sizeof(smem_i)); iter = xcalloc(1, sizeof(smem_i));
iter->bwt = bwt; iter->bwt = bwt;
iter->tmpvec[0] = calloc(1, sizeof(bwtintv_v)); iter->tmpvec[0] = xcalloc(1, sizeof(bwtintv_v));
iter->tmpvec[1] = calloc(1, sizeof(bwtintv_v)); iter->tmpvec[1] = xcalloc(1, sizeof(bwtintv_v));
iter->matches = calloc(1, sizeof(bwtintv_v)); iter->matches = xcalloc(1, sizeof(bwtintv_v));
return iter; return iter;
} }
@ -78,7 +79,7 @@ int main_fastmap(int argc, char *argv[])
fp = gzopen(argv[optind + 1], "r"); fp = gzopen(argv[optind + 1], "r");
seq = kseq_init(fp); seq = kseq_init(fp);
{ // load the packed sequences, BWT and SA { // load the packed sequences, BWT and SA
char *tmp = calloc(strlen(argv[optind]) + 5, 1); char *tmp = xcalloc(strlen(argv[optind]) + 5, 1);
strcat(strcpy(tmp, argv[optind]), ".bwt"); strcat(strcpy(tmp, argv[optind]), ".bwt");
bwt = bwt_restore_bwt(tmp); bwt = bwt_restore_bwt(tmp);
strcat(strcpy(tmp, argv[optind]), ".sa"); strcat(strcpy(tmp, argv[optind]), ".sa");

6
is.c
View File

@ -25,6 +25,7 @@
*/ */
#include <stdlib.h> #include <stdlib.h>
#include "utils.h"
typedef unsigned char ubyte_t; typedef unsigned char ubyte_t;
#define chr(i) (cs == sizeof(int) ? ((const int *)T)[i]:((const unsigned char *)T)[i]) #define chr(i) (cs == sizeof(int) ? ((const int *)T)[i]:((const unsigned char *)T)[i])
@ -204,8 +205,9 @@ int is_sa(const ubyte_t *T, int *SA, int n)
int is_bwt(ubyte_t *T, int n) int is_bwt(ubyte_t *T, int n)
{ {
int *SA, i, primary = 0; int *SA, i, primary = 0;
SA = (int*)calloc(n+1, sizeof(int)); SA = (int*)xcalloc(n+1, sizeof(int));
is_sa(T, SA, n);
if (is_sa(T, SA, n)) err_fatal_simple("is_sa failed");
for (i = 0; i <= n; ++i) { for (i = 0; i <= n; ++i) {
if (SA[i] == 0) primary = i; if (SA[i] == 0) primary = i;

13
khash.h
View File

@ -95,6 +95,7 @@ int main() {
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <limits.h> #include <limits.h>
#include "utils.h"
/* compipler specific configuration */ /* compipler specific configuration */
@ -147,7 +148,7 @@ static const double __ac_HASH_UPPER = 0.77;
khval_t *vals; \ khval_t *vals; \
} kh_##name##_t; \ } kh_##name##_t; \
static inline kh_##name##_t *kh_init_##name() { \ static inline kh_##name##_t *kh_init_##name() { \
return (kh_##name##_t*)calloc(1, sizeof(kh_##name##_t)); \ return (kh_##name##_t*)xcalloc(1, sizeof(kh_##name##_t)); \
} \ } \
static inline void kh_destroy_##name(kh_##name##_t *h) \ static inline void kh_destroy_##name(kh_##name##_t *h) \
{ \ { \
@ -188,12 +189,12 @@ static const double __ac_HASH_UPPER = 0.77;
new_n_buckets = __ac_prime_list[t+1]; \ new_n_buckets = __ac_prime_list[t+1]; \
if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0; \ if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0; \
else { \ else { \
new_flags = (khint32_t*)malloc(((new_n_buckets>>4) + 1) * sizeof(khint32_t)); \ new_flags = (khint32_t*)xmalloc(((new_n_buckets>>4) + 1) * sizeof(khint32_t)); \
memset(new_flags, 0xaa, ((new_n_buckets>>4) + 1) * sizeof(khint32_t)); \ memset(new_flags, 0xaa, ((new_n_buckets>>4) + 1) * sizeof(khint32_t)); \
if (h->n_buckets < new_n_buckets) { \ if (h->n_buckets < new_n_buckets) { \
h->keys = (khkey_t*)realloc(h->keys, new_n_buckets * sizeof(khkey_t)); \ h->keys = (khkey_t*)xrealloc(h->keys, new_n_buckets * sizeof(khkey_t)); \
if (kh_is_map) \ if (kh_is_map) \
h->vals = (khval_t*)realloc(h->vals, new_n_buckets * sizeof(khval_t)); \ h->vals = (khval_t*)xrealloc(h->vals, new_n_buckets * sizeof(khval_t)); \
} \ } \
} \ } \
} \ } \
@ -227,9 +228,9 @@ static const double __ac_HASH_UPPER = 0.77;
} \ } \
} \ } \
if (h->n_buckets > new_n_buckets) { \ if (h->n_buckets > new_n_buckets) { \
h->keys = (khkey_t*)realloc(h->keys, new_n_buckets * sizeof(khkey_t)); \ h->keys = (khkey_t*)xrealloc(h->keys, new_n_buckets * sizeof(khkey_t)); \
if (kh_is_map) \ if (kh_is_map) \
h->vals = (khval_t*)realloc(h->vals, new_n_buckets * sizeof(khval_t)); \ h->vals = (khval_t*)xrealloc(h->vals, new_n_buckets * sizeof(khval_t)); \
} \ } \
free(h->flags); \ free(h->flags); \
h->flags = new_flags; \ h->flags = new_flags; \

13
kseq.h
View File

@ -29,6 +29,7 @@
#include <ctype.h> #include <ctype.h>
#include <string.h> #include <string.h>
#include <stdlib.h> #include <stdlib.h>
#include "utils.h"
#define __KS_TYPE(type_t) \ #define __KS_TYPE(type_t) \
typedef struct __kstream_t { \ typedef struct __kstream_t { \
@ -43,9 +44,9 @@
#define __KS_BASIC(type_t, __bufsize) \ #define __KS_BASIC(type_t, __bufsize) \
static inline kstream_t *ks_init(type_t f) \ static inline kstream_t *ks_init(type_t f) \
{ \ { \
kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t)); \ kstream_t *ks = (kstream_t*)xcalloc(1, sizeof(kstream_t)); \
ks->f = f; \ ks->f = f; \
ks->buf = (char*)malloc(__bufsize); \ ks->buf = (char*)xmalloc(__bufsize); \
return ks; \ return ks; \
} \ } \
static inline void ks_destroy(kstream_t *ks) \ static inline void ks_destroy(kstream_t *ks) \
@ -107,7 +108,7 @@ typedef struct __kstring_t {
if (str->m - str->l < i - ks->begin + 1) { \ if (str->m - str->l < i - ks->begin + 1) { \
str->m = str->l + (i - ks->begin) + 1; \ str->m = str->l + (i - ks->begin) + 1; \
kroundup32(str->m); \ kroundup32(str->m); \
str->s = (char*)realloc(str->s, str->m); \ str->s = (char*)xrealloc(str->s, str->m); \
} \ } \
memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \ memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \
str->l = str->l + (i - ks->begin); \ str->l = str->l + (i - ks->begin); \
@ -130,7 +131,7 @@ typedef struct __kstring_t {
#define __KSEQ_BASIC(type_t) \ #define __KSEQ_BASIC(type_t) \
static inline kseq_t *kseq_init(type_t fd) \ static inline kseq_t *kseq_init(type_t fd) \
{ \ { \
kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t)); \ kseq_t *s = (kseq_t*)xcalloc(1, sizeof(kseq_t)); \
s->f = ks_init(fd); \ s->f = ks_init(fd); \
return s; \ return s; \
} \ } \
@ -170,7 +171,7 @@ typedef struct __kstring_t {
if (seq->seq.l + 1 >= seq->seq.m) { /* double the memory */ \ if (seq->seq.l + 1 >= seq->seq.m) { /* double the memory */ \
seq->seq.m = seq->seq.l + 2; \ seq->seq.m = seq->seq.l + 2; \
kroundup32(seq->seq.m); /* rounded to next closest 2^k */ \ kroundup32(seq->seq.m); /* rounded to next closest 2^k */ \
seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \ seq->seq.s = (char*)xrealloc(seq->seq.s, seq->seq.m); \
} \ } \
seq->seq.s[seq->seq.l++] = (char)c; \ seq->seq.s[seq->seq.l++] = (char)c; \
} \ } \
@ -180,7 +181,7 @@ typedef struct __kstring_t {
if (c != '+') return seq->seq.l; /* FASTA */ \ if (c != '+') return seq->seq.l; /* FASTA */ \
if (seq->qual.m < seq->seq.m) { /* allocate enough memory */ \ if (seq->qual.m < seq->seq.m) { /* allocate enough memory */ \
seq->qual.m = seq->seq.m; \ seq->qual.m = seq->seq.m; \
seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \ seq->qual.s = (char*)xrealloc(seq->qual.s, seq->qual.m); \
} \ } \
while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \ while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \
if (c == -1) return -2; /* we should not stop here */ \ if (c == -1) return -2; /* we should not stop here */ \

5
ksort.h
View File

@ -57,6 +57,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "utils.h"
typedef struct { typedef struct {
void *left, *right; void *left, *right;
@ -72,7 +73,7 @@ typedef struct {
int curr, shift; \ int curr, shift; \
\ \
a2[0] = array; \ a2[0] = array; \
a2[1] = temp? temp : (type_t*)malloc(sizeof(type_t) * n); \ a2[1] = temp? temp : (type_t*)xmalloc(sizeof(type_t) * n); \
for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \ for (curr = 0, shift = 0; (1ul<<shift) < n; ++shift) { \
a = a2[curr]; b = a2[1-curr]; \ a = a2[curr]; b = a2[1-curr]; \
if (shift == 0) { \ if (shift == 0) { \
@ -182,7 +183,7 @@ typedef struct {
return; \ return; \
} \ } \
for (d = 2; 1ul<<d < n; ++d); \ for (d = 2; 1ul<<d < n; ++d); \
stack = (ks_isort_stack_t*)malloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \ stack = (ks_isort_stack_t*)xmalloc(sizeof(ks_isort_stack_t) * ((sizeof(size_t)*d)+2)); \
top = stack; s = a; t = a + (n-1); d <<= 1; \ top = stack; s = a; t = a + (n-1); d <<= 1; \
while (1) { \ while (1) { \
if (s < t) { \ if (s < t) { \

3
kstring.c
View File

@ -1,6 +1,7 @@
#include <stdarg.h> #include <stdarg.h>
#include <stdio.h> #include <stdio.h>
#include "kstring.h" #include "kstring.h"
#include "utils.h"
int ksprintf(kstring_t *s, const char *fmt, ...) int ksprintf(kstring_t *s, const char *fmt, ...)
{ {
@ -12,7 +13,7 @@ int ksprintf(kstring_t *s, const char *fmt, ...)
if (l + 1 > s->m - s->l) { if (l + 1 > s->m - s->l) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)realloc(s->s, s->m); s->s = (char*)xrealloc(s->s, s->m);
va_start(ap, fmt); va_start(ap, fmt);
l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap); l = vsnprintf(s->s + s->l, s->m - s->l, fmt, ap);
} }

5
kstring.h
View File

@ -3,6 +3,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include "utils.h"
#ifndef kroundup32 #ifndef kroundup32
#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
@ -22,7 +23,7 @@ static inline int kputs(const char *p, kstring_t *s)
if (s->l + l + 1 >= s->m) { if (s->l + l + 1 >= s->m) {
s->m = s->l + l + 2; s->m = s->l + l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)realloc(s->s, s->m); s->s = (char*)xrealloc(s->s, s->m);
} }
strcpy(s->s + s->l, p); strcpy(s->s + s->l, p);
s->l += l; s->l += l;
@ -34,7 +35,7 @@ static inline int kputc(int c, kstring_t *s)
if (s->l + 1 >= s->m) { if (s->l + 1 >= s->m) {
s->m = s->l + 2; s->m = s->l + 2;
kroundup32(s->m); kroundup32(s->m);
s->s = (char*)realloc(s->s, s->m); s->s = (char*)xrealloc(s->s, s->m);
} }
s->s[s->l++] = c; s->s[s->l++] = c;
s->s[s->l] = 0; s->s[s->l] = 0;

9
ksw.c
View File

@ -28,6 +28,7 @@
#include <stdint.h> #include <stdint.h>
#include <emmintrin.h> #include <emmintrin.h>
#include "ksw.h" #include "ksw.h"
#include "utils.h"
#ifdef __GNUC__ #ifdef __GNUC__
#define LIKELY(x) __builtin_expect((x),1) #define LIKELY(x) __builtin_expect((x),1)
@ -51,7 +52,7 @@ ksw_query_t *ksw_qinit(int size, int qlen, const uint8_t *query, int m, const in
size = size > 1? 2 : 1; size = size > 1? 2 : 1;
p = 8 * (3 - size); // # values per __m128i p = 8 * (3 - size); // # values per __m128i
slen = (qlen + p - 1) / p; // segmented length slen = (qlen + p - 1) / p; // segmented length
q = malloc(sizeof(ksw_query_t) + 256 + 16 * slen * (m + 4)); // a single block of memory q = xmalloc(sizeof(ksw_query_t) + 256 + 16 * slen * (m + 4)); // a single block of memory
q->qp = (__m128i*)(((size_t)q + sizeof(ksw_query_t) + 15) >> 4 << 4); // align memory q->qp = (__m128i*)(((size_t)q + sizeof(ksw_query_t) + 15) >> 4 << 4); // align memory
q->H0 = q->qp + slen * m; q->H0 = q->qp + slen * m;
q->H1 = q->H0 + slen; q->H1 = q->H0 + slen;
@ -169,7 +170,7 @@ end_loop16:
if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append
if (n_b == m_b) { if (n_b == m_b) {
m_b = m_b? m_b<<1 : 8; m_b = m_b? m_b<<1 : 8;
b = realloc(b, 8 * m_b); b = xrealloc(b, 8 * m_b);
} }
b[n_b++] = (uint64_t)imax<<32 | i; b[n_b++] = (uint64_t)imax<<32 | i;
} else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
@ -264,7 +265,7 @@ end_loop8:
if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) {
if (n_b == m_b) { if (n_b == m_b) {
m_b = m_b? m_b<<1 : 8; m_b = m_b? m_b<<1 : 8;
b = realloc(b, 8 * m_b); b = xrealloc(b, 8 * m_b);
} }
b[n_b++] = (uint64_t)imax<<32 | i; b[n_b++] = (uint64_t)imax<<32 | i;
} else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last } else if ((int)(b[n_b-1]>>32) < imax) b[n_b-1] = (uint64_t)imax<<32 | i; // modify the last
@ -310,7 +311,7 @@ int ksw_sse2(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a)
#include <stdio.h> #include <stdio.h>
#include <zlib.h> #include <zlib.h>
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) KSEQ_INIT(gzFile, err_gzread)
unsigned char seq_nt4_table[256] = { unsigned char seq_nt4_table[256] = {
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,

8
kvec.h
View File

@ -60,7 +60,7 @@ int main() {
#define kv_size(v) ((v).n) #define kv_size(v) ((v).n)
#define kv_max(v) ((v).m) #define kv_max(v) ((v).m)
#define kv_resize(type, v, s) ((v).m = (s), (v).a = (type*)realloc((v).a, sizeof(type) * (v).m)) #define kv_resize(type, v, s) ((v).m = (s), (v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m))
#define kv_copy(type, v1, v0) do { \ #define kv_copy(type, v1, v0) do { \
if ((v1).m < (v0).n) kv_resize(type, v1, (v0).n); \ if ((v1).m < (v0).n) kv_resize(type, v1, (v0).n); \
@ -71,19 +71,19 @@ int main() {
#define kv_push(type, v, x) do { \ #define kv_push(type, v, x) do { \
if ((v).n == (v).m) { \ if ((v).n == (v).m) { \
(v).m = (v).m? (v).m<<1 : 2; \ (v).m = (v).m? (v).m<<1 : 2; \
(v).a = (type*)realloc((v).a, sizeof(type) * (v).m); \ (v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m); \
} \ } \
(v).a[(v).n++] = (x); \ (v).a[(v).n++] = (x); \
} while (0) } while (0)
#define kv_pushp(type, v) (((v).n == (v).m)? \ #define kv_pushp(type, v) (((v).n == (v).m)? \
((v).m = ((v).m? (v).m<<1 : 2), \ ((v).m = ((v).m? (v).m<<1 : 2), \
(v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0) \ (v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m), 0) \
: 0), ((v).a + ((v).n++)) : 0), ((v).a + ((v).n++))
#define kv_a(type, v, i) ((v).m <= (size_t)(i)? \ #define kv_a(type, v, i) ((v).m <= (size_t)(i)? \
((v).m = (v).n = (i) + 1, kv_roundup32((v).m), \ ((v).m = (v).n = (i) + 1, kv_roundup32((v).m), \
(v).a = (type*)realloc((v).a, sizeof(type) * (v).m), 0) \ (v).a = (type*)xrealloc((v).a, sizeof(type) * (v).m), 0) \
: (v).n <= (size_t)(i)? (v).n = (i) \ : (v).n <= (size_t)(i)? (v).n = (i) \
: 0), (v).a[(i)] : 0), (v).a[(i)]

12
simple_dp.c
View File

@ -8,7 +8,7 @@
#include "utils.h" #include "utils.h"
#include "kseq.h" #include "kseq.h"
KSEQ_INIT(gzFile, gzread) KSEQ_INIT(gzFile, err_gzread)
typedef struct { typedef struct {
int l; int l;
@ -64,20 +64,20 @@ static seqs_t *load_seqs(const char *fn)
fp = xzopen(fn, "r"); fp = xzopen(fn, "r");
seq = kseq_init(fp); seq = kseq_init(fp);
s = (seqs_t*)calloc(1, sizeof(seqs_t)); s = (seqs_t*)xcalloc(1, sizeof(seqs_t));
s->m_seqs = 256; s->m_seqs = 256;
s->seqs = (seq1_t*)calloc(s->m_seqs, sizeof(seq1_t)); s->seqs = (seq1_t*)xcalloc(s->m_seqs, sizeof(seq1_t));
while ((l = kseq_read(seq)) >= 0) { while ((l = kseq_read(seq)) >= 0) {
if (s->n_seqs == s->m_seqs) { if (s->n_seqs == s->m_seqs) {
s->m_seqs <<= 1; s->m_seqs <<= 1;
s->seqs = (seq1_t*)realloc(s->seqs, s->m_seqs * sizeof(seq1_t)); s->seqs = (seq1_t*)xrealloc(s->seqs, s->m_seqs * sizeof(seq1_t));
} }
p = s->seqs + (s->n_seqs++); p = s->seqs + (s->n_seqs++);
p->l = seq->seq.l; p->l = seq->seq.l;
p->s = (unsigned char*)malloc(p->l + 1); p->s = (unsigned char*)xmalloc(p->l + 1);
memcpy(p->s, seq->seq.s, p->l); memcpy(p->s, seq->seq.s, p->l);
p->s[p->l] = 0; p->s[p->l] = 0;
p->n = strdup((const char*)seq->name.s); p->n = xstrdup((const char*)seq->name.s);
} }
kseq_destroy(seq); kseq_destroy(seq);
gzclose(fp); gzclose(fp);

37
stdaln.c
View File

@ -28,6 +28,7 @@
#include <string.h> #include <string.h>
#include <stdint.h> #include <stdint.h>
#include "stdaln.h" #include "stdaln.h"
#include "utils.h"
/* char -> 17 (=16+1) nucleotides */ /* char -> 17 (=16+1) nucleotides */
unsigned char aln_nt16_table[256] = { unsigned char aln_nt16_table[256] = {
@ -232,7 +233,7 @@ AlnParam aln_param_aa2aa = { 10, 2, 2, aln_sm_blosum62, 22, 50 };
AlnAln *aln_init_AlnAln() AlnAln *aln_init_AlnAln()
{ {
AlnAln *aa; AlnAln *aa;
aa = (AlnAln*)malloc(sizeof(AlnAln)); aa = (AlnAln*)xmalloc(sizeof(AlnAln));
aa->path = 0; aa->path = 0;
aa->out1 = aa->out2 = aa->outm = 0; aa->out1 = aa->out2 = aa->outm = 0;
aa->path_len = 0; aa->path_len = 0;
@ -382,13 +383,13 @@ int aln_global_core(unsigned char *seq1, int len1, unsigned char *seq2, int len2
/* allocate memory */ /* allocate memory */
end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1); end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1);
dpcell = (dpcell_t**)malloc(sizeof(dpcell_t*) * (len2 + 1)); dpcell = (dpcell_t**)xmalloc(sizeof(dpcell_t*) * (len2 + 1));
for (j = 0; j <= len2; ++j) for (j = 0; j <= len2; ++j)
dpcell[j] = (dpcell_t*)malloc(sizeof(dpcell_t) * end); dpcell[j] = (dpcell_t*)xmalloc(sizeof(dpcell_t) * end);
for (j = b2 + 1; j <= len2; ++j) for (j = b2 + 1; j <= len2; ++j)
dpcell[j] -= j - b2; dpcell[j] -= j - b2;
curr = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1)); curr = (dpscore_t*)xmalloc(sizeof(dpscore_t) * (len1 + 1));
last = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1)); last = (dpscore_t*)xmalloc(sizeof(dpscore_t) * (len1 + 1));
/* set first row */ /* set first row */
SET_INF(*curr); curr->M = 0; SET_INF(*curr); curr->M = 0;
@ -556,11 +557,11 @@ int aln_local_core(unsigned char *seq1, int len1, unsigned char *seq2, int len2,
if (len1 == 0 || len2 == 0) return -1; if (len1 == 0 || len2 == 0) return -1;
/* allocate memory */ /* allocate memory */
suba = (int*)malloc(sizeof(int) * (len2 + 1)); suba = (int*)xmalloc(sizeof(int) * (len2 + 1));
eh = (NT_LOCAL_SCORE*)malloc(sizeof(NT_LOCAL_SCORE) * (len1 + 1)); eh = (NT_LOCAL_SCORE*)xmalloc(sizeof(NT_LOCAL_SCORE) * (len1 + 1));
s_array = (int**)malloc(sizeof(int*) * N_MATRIX_ROW); s_array = (int**)xmalloc(sizeof(int*) * N_MATRIX_ROW);
for (i = 0; i != N_MATRIX_ROW; ++i) for (i = 0; i != N_MATRIX_ROW; ++i)
s_array[i] = (int*)malloc(sizeof(int) * len1); s_array[i] = (int*)xmalloc(sizeof(int) * len1);
/* initialization */ /* initialization */
aln_init_score_array(seq1, len1, N_MATRIX_ROW, score_matrix, s_array); aln_init_score_array(seq1, len1, N_MATRIX_ROW, score_matrix, s_array);
q = gap_open; q = gap_open;
@ -773,9 +774,9 @@ AlnAln *aln_stdaln_aux(const char *seq1, const char *seq2, const AlnParam *ap,
if (len2 < 0) len2 = strlen(seq2); if (len2 < 0) len2 = strlen(seq2);
aa = aln_init_AlnAln(); aa = aln_init_AlnAln();
seq11 = (unsigned char*)malloc(sizeof(unsigned char) * len1); seq11 = (unsigned char*)xmalloc(sizeof(unsigned char) * len1);
seq22 = (unsigned char*)malloc(sizeof(unsigned char) * len2); seq22 = (unsigned char*)xmalloc(sizeof(unsigned char) * len2);
aa->path = (path_t*)malloc(sizeof(path_t) * (len1 + len2 + 1)); aa->path = (path_t*)xmalloc(sizeof(path_t) * (len1 + len2 + 1));
if (ap->row < 10) { /* 4-nucleotide alignment */ if (ap->row < 10) { /* 4-nucleotide alignment */
for (i = 0; i < len1; ++i) for (i = 0; i < len1; ++i)
@ -805,9 +806,9 @@ AlnAln *aln_stdaln_aux(const char *seq1, const char *seq2, const AlnParam *ap,
aa->score = score; aa->score = score;
if (thres > 0) { if (thres > 0) {
out1 = aa->out1 = (char*)malloc(sizeof(char) * (aa->path_len + 1)); out1 = aa->out1 = (char*)xmalloc(sizeof(char) * (aa->path_len + 1));
out2 = aa->out2 = (char*)malloc(sizeof(char) * (aa->path_len + 1)); out2 = aa->out2 = (char*)xmalloc(sizeof(char) * (aa->path_len + 1));
outm = aa->outm = (char*)malloc(sizeof(char) * (aa->path_len + 1)); outm = aa->outm = (char*)xmalloc(sizeof(char) * (aa->path_len + 1));
--seq1; --seq2; --seq1; --seq2;
--seq11; --seq22; --seq11; --seq22;
@ -881,10 +882,10 @@ int aln_extend_core(unsigned char *seq1, int len1, unsigned char *seq2, int len2
if (len1 == 0 || len2 == 0) return -1; if (len1 == 0 || len2 == 0) return -1;
/* allocate memory */ /* allocate memory */
mem = _mem? _mem : calloc((len1 + 2) * (N_MATRIX_ROW + 1), 4); mem = _mem? _mem : xcalloc((len1 + 2) * (N_MATRIX_ROW + 1), 4);
_p = mem; _p = mem;
eh = (uint32_t*)_p, _p += 4 * (len1 + 2); eh = (uint32_t*)_p, _p += 4 * (len1 + 2);
s_array = calloc(N_MATRIX_ROW, sizeof(void*)); s_array = xcalloc(N_MATRIX_ROW, sizeof(void*));
for (i = 0; i != N_MATRIX_ROW; ++i) for (i = 0; i != N_MATRIX_ROW; ++i)
s_array[i] = (int32_t*)_p, _p += 4 * len1; s_array[i] = (int32_t*)_p, _p += 4 * len1;
/* initialization */ /* initialization */
@ -1024,7 +1025,7 @@ uint32_t *aln_path2cigar32(const path_t *path, int path_len, int *n_cigar)
last_type = path[i].ctype; last_type = path[i].ctype;
} }
*n_cigar = n; *n_cigar = n;
cigar = (uint32_t*)malloc(*n_cigar * 4); cigar = (uint32_t*)xmalloc(*n_cigar * 4);
cigar[0] = 1u << 4 | path[path_len-1].ctype; cigar[0] = 1u << 4 | path[path_len-1].ctype;
last_type = path[path_len-1].ctype; last_type = path[path_len-1].ctype;