updated to the latest kseq.h

Makefile

@@ -5,7 +5,7 @@ CXXFLAGS=	$(CFLAGS)
 AR=			ar
 DFLAGS=		-DHAVE_PTHREAD #-D_NO_SSE2 #-D_FILE_OFFSET_BITS=64
 LOBJS=		bwa.o bamlite.o utils.o bwt.o bwtio.o bwtaln.o bwtgap.o bntseq.o bwamem.o stdaln.o \
-			bwaseqio.o bwase.o kstring.o
+			bseq.o bwaseqio.o bwase.o kstring.o
 AOBJS=		QSufSort.o bwt_gen.o \
 			is.o bwtmisc.o bwtindex.o ksw.o simple_dp.o \
 			bwape.o cs2nt.o \

bntseq.c

@@ -35,7 +35,7 @@
 #include "utils.h"
 
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_DECLARE(gzFile)
 
 unsigned char nst_nt4_table[256] = {
 	4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 

bseq.c

@@ -0,0 +1,4 @@
+#include <zlib.h>
+#include "bseq.h"
+#include "kseq.h"
+KSEQ_INIT2(, gzFile, gzread)

bseq.h

@@ -0,0 +1,15 @@
+#ifndef BATCHSEQ_H_
+#define BATCHSEQ_H_
+
+typedef struct {
+	char *name, *comment, *seq, *qual;
+} bseq1_t;
+
+typedef struct {
+	int n, m;
+	bseq1_t *seqs;
+} bseq_t;
+
+int bseq_read(int chunk_size, bseq_t *bs);
+
+#endif

bwamem.c

@@ -274,6 +274,32 @@ int mem_chain_flt(const mem_opt_t *opt, int n_chn, mem_chain1_t *chains)
 	return n;
 }
 
+#define alnreg_lt(a, b) ((a).score > (b).score)
+KSORT_INIT(mem_ar, mem_alnreg_t, alnreg_lt)
+
+int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a)
+{ // similar to the loop in mem_chain_flt()
+	int i, j, m;
+	if (n <= 1) return n;
+	ks_introsort(mem_ar, n, a);
+	for (i = 0; i < n; ++i) a[i].sub = 0;
+	for (i = 1, m = 1; i < n; ++i) {
+		for (j = 0; j < m; ++j) {
+			int b_max = a[j].qb > a[i].qb? a[j].qb : a[i].qb;
+			int e_min = a[j].qe < a[i].qe? a[j].qe : a[i].qe;
+			if (e_min > b_max) { // have overlap
+				int min_l = a[i].qe - a[i].qb < a[j].qe - a[j].qb? a[i].qe - a[i].qb : a[j].qe - a[j].qb;
+				if (e_min - b_max >= min_l * opt->mask_level) { // significant overlap
+					if (a[j].sub == 0) a[j].sub = a[i].score;
+					break;
+				}
+			}
+		}
+		if (j == m) a[m++] = a[i];
+	}
+	return m;
+}
+
 /****************************************
  * Construct the alignment from a chain *
  ****************************************/
@@ -388,3 +414,7 @@ ret_gen_cigar:
 	free(rseq);
 	return cigar;
 }
+
+/****************
+ * Sequence I/O *
+ ****************/

bwamem.h

@@ -31,7 +31,7 @@ typedef struct {
 
 typedef struct {
 	int64_t rb, re;
-	int score, qb, qe;
+	int score, qb, qe, sub;
 } mem_alnreg_t;
 
 #ifdef __cplusplus

bwaseqio.c

@@ -5,7 +5,7 @@
 #include "bamlite.h"
 
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_DECLARE(gzFile)
 
 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 };

bwtsw2_aux.c

@@ -15,7 +15,7 @@
 #include "kstring.h"
 
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_DECLARE(gzFile)
 
 #include "ksort.h"
 #define __left_lt(a, b) ((a).end > (b).end)

fastmap.c

@@ -7,7 +7,7 @@
 #include "bwamem.h"
 #include "kvec.h"
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_DECLARE(gzFile)
 
 extern unsigned char nst_nt4_table[256];
 

kseq.h

@@ -1,6 +1,6 @@
 /* The MIT License
 
-   Copyright (c) 2008, by Heng Li <lh3@sanger.ac.uk>
+   Copyright (c) 2008, 2009, 2011 Attractive Chaos <attractor@live.co.uk>
 
    Permission is hereby granted, free of charge, to any person obtaining
    a copy of this software and associated documentation files (the
@@ -23,6 +23,8 @@
    SOFTWARE.
 */
 
+/* Last Modified: 05MAR2012 */
+
 #ifndef AC_KSEQ_H
 #define AC_KSEQ_H
 
@@ -30,9 +32,14 @@
 #include <string.h>
 #include <stdlib.h>
 
+#define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r
+#define KS_SEP_TAB   1 // isspace() && !' '
+#define KS_SEP_LINE  2 // line separator: "\n" (Unix) or "\r\n" (Windows)
+#define KS_SEP_MAX   2
+
 #define __KS_TYPE(type_t)						\
 	typedef struct __kstream_t {				\
-		char *buf;								\
+		unsigned char *buf;						\
 		int begin, end, is_eof;					\
 		type_t f;								\
 	} kstream_t;
@@ -45,7 +52,7 @@
 	{																\
 		kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t));	\
 		ks->f = f;													\
-		ks->buf = (char*)malloc(__bufsize);							\
+		ks->buf = (unsigned char*)malloc(__bufsize);				\
 		return ks;													\
 	}																\
 	static inline void ks_destroy(kstream_t *ks)					\
@@ -82,10 +89,10 @@ typedef struct __kstring_t {
 #endif
 
 #define __KS_GETUNTIL(__read, __bufsize)								\
-	static int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \
+	static int ks_getuntil2(kstream_t *ks, int delimiter, kstring_t *str, int *dret, int append) \
 	{																	\
 		if (dret) *dret = 0;											\
-		str->l = 0;														\
+		str->l = append? str->l : 0;									\
 		if (ks->begin >= ks->end && ks->is_eof) return -1;				\
 		for (;;) {														\
 			int i;														\
@@ -97,14 +104,20 @@ typedef struct __kstring_t {
 					if (ks->end == 0) break;							\
 				} else break;											\
 			}															\
-			if (delimiter) {											\
+			if (delimiter == KS_SEP_LINE) { \
+				for (i = ks->begin; i < ks->end; ++i) \
+					if (ks->buf[i] == '\n') break; \
+			} else if (delimiter > KS_SEP_MAX) {						\
 				for (i = ks->begin; i < ks->end; ++i)					\
 					if (ks->buf[i] == delimiter) break;					\
-			} else {													\
+			} else if (delimiter == KS_SEP_SPACE) {						\
 				for (i = ks->begin; i < ks->end; ++i)					\
 					if (isspace(ks->buf[i])) break;						\
-			}															\
+			} else if (delimiter == KS_SEP_TAB) {						\
-			if (str->m - str->l < i - ks->begin + 1) {					\
+				for (i = ks->begin; i < ks->end; ++i)					\
+					if (isspace(ks->buf[i]) && ks->buf[i] != ' ') break; \
+			} else i = 0; /* never come to here! */						\
+			if (str->m - str->l < (size_t)(i - ks->begin + 1)) {		\
 				str->m = str->l + (i - ks->begin) + 1;					\
 				kroundup32(str->m);										\
 				str->s = (char*)realloc(str->s, str->m);				\
@@ -117,9 +130,15 @@ typedef struct __kstring_t {
 				break;													\
 			}															\
 		}																\
+		if (str->s == 0) {												\
+			str->m = 1;													\
+			str->s = (char*)calloc(1, 1);								\
+		} else if (delimiter == KS_SEP_LINE && str->l > 1 && str->s[str->l-1] == '\r') --str->l; \
 		str->s[str->l] = '\0';											\
 		return str->l;													\
-	}
+	} \
+	static inline int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \
+	{ return ks_getuntil2(ks, delimiter, str, dret, 0); }
 
 #define KSTREAM_INIT(type_t, __read, __bufsize) \
 	__KS_TYPE(type_t)							\
@@ -127,19 +146,16 @@ typedef struct __kstring_t {
 	__KS_GETC(__read, __bufsize)				\
 	__KS_GETUNTIL(__read, __bufsize)
 
-#define __KSEQ_BASIC(type_t)											\
+#define kseq_rewind(ks) ((ks)->last_char = (ks)->f->is_eof = (ks)->f->begin = (ks)->f->end = 0)
-	static inline kseq_t *kseq_init(type_t fd)							\
+
+#define __KSEQ_BASIC(SCOPE, type_t)										\
+	SCOPE kseq_t *kseq_init(type_t fd)									\
 	{																	\
 		kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t));					\
 		s->f = ks_init(fd);												\
 		return s;														\
 	}																	\
-	static inline void kseq_rewind(kseq_t *ks)							\
+	SCOPE void kseq_destroy(kseq_t *ks)									\
-	{																	\
-		ks->last_char = 0;												\
-		ks->f->is_eof = ks->f->begin = ks->f->end = 0;					\
-	}																	\
-	static inline void kseq_destroy(kseq_t *ks)							\
 	{																	\
 		if (!ks) return;												\
 		free(ks->name.s); free(ks->comment.s); free(ks->seq.s);	free(ks->qual.s); \
@@ -152,44 +168,46 @@ typedef struct __kstring_t {
    -1   end-of-file
    -2   truncated quality string
  */
-#define __KSEQ_READ														\
+#define __KSEQ_READ(SCOPE) \
-	static int kseq_read(kseq_t *seq)									\
+	SCOPE int kseq_read(kseq_t *seq) \
-	{																	\
+	{ \
-		int c;															\
+		int c; \
-		kstream_t *ks = seq->f;											\
+		kstream_t *ks = seq->f; \
 		if (seq->last_char == 0) { /* then jump to the next header line */ \
-			while ((c = ks_getc(ks)) != -1 && c != '>' && c != '@');	\
+			while ((c = ks_getc(ks)) != -1 && c != '>' && c != '@'); \
-			if (c == -1) return -1; /* end of file */					\
+			if (c == -1) return -1; /* end of file */ \
-			seq->last_char = c;											\
+			seq->last_char = c; \
-		} /* the first header char has been read */						\
+		} /* else: the first header char has been read in the previous call */ \
-		seq->comment.l = seq->seq.l = seq->qual.l = 0;					\
+		seq->comment.l = seq->seq.l = seq->qual.l = 0; /* reset all members */ \
-		if (ks_getuntil(ks, 0, &seq->name, &c) < 0) return -1;			\
+		if (ks_getuntil(ks, 0, &seq->name, &c) < 0) return -1; /* normal exit: EOF */ \
-		if (c != '\n') ks_getuntil(ks, '\n', &seq->comment, 0);			\
+		if (c != '\n') ks_getuntil(ks, KS_SEP_LINE, &seq->comment, 0); /* read FASTA/Q comment */ \
+		if (seq->seq.s == 0) { /* we can do this in the loop below, but that is slower */ \
+			seq->seq.m = 256; \
+			seq->seq.s = (char*)malloc(seq->seq.m); \
+		} \
 		while ((c = ks_getc(ks)) != -1 && c != '>' && c != '+' && c != '@') { \
-			if (isgraph(c)) { /* printable non-space character */		\
+			if (c == '\n') continue; /* skip empty lines */ \
-				if (seq->seq.l + 1 >= seq->seq.m) { /* double the memory */ \
+			seq->seq.s[seq->seq.l++] = c; /* this is safe: we always have enough space for 1 char */ \
-					seq->seq.m = seq->seq.l + 2;						\
+			ks_getuntil2(ks, KS_SEP_LINE, &seq->seq, 0, 1); /* read the rest of the line */ \
-					kroundup32(seq->seq.m); /* rounded to next closest 2^k */ \
+		} \
-					seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \
-				}														\
-				seq->seq.s[seq->seq.l++] = (char)c;						\
-			}															\
-		}																\
 		if (c == '>' || c == '@') seq->last_char = c; /* the first header char has been read */	\
-		seq->seq.s[seq->seq.l] = 0;	/* null terminated string */		\
+		if (seq->seq.l + 1 >= seq->seq.m) { /* seq->seq.s[seq->seq.l] below may be out of boundary */ \
-		if (c != '+') return seq->seq.l; /* FASTA */					\
+			seq->seq.m = seq->seq.l + 2; \
-		if (seq->qual.m < seq->seq.m) {	/* allocate enough memory */	\
+			kroundup32(seq->seq.m); /* rounded to the next closest 2^k */ \
-			seq->qual.m = seq->seq.m;									\
+			seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \
-			seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m);		\
+		} \
-		}																\
+		seq->seq.s[seq->seq.l] = 0;	/* null terminated string */ \
+		if (c != '+') return seq->seq.l; /* FASTA */ \
+		if (seq->qual.m < seq->seq.m) {	/* allocate memory for qual in case insufficient */ \
+			seq->qual.m = seq->seq.m; \
+			seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \
+		} \
 		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; /* error: no quality string */ \
-		while ((c = ks_getc(ks)) != -1 && seq->qual.l < seq->seq.l)		\
+		while (ks_getuntil2(ks, KS_SEP_LINE, &seq->qual, 0, 1) >= 0 && seq->qual.l < seq->seq.l); \
-			if (c >= 33 && c <= 127) seq->qual.s[seq->qual.l++] = (unsigned char)c;	\
-		seq->qual.s[seq->qual.l] = 0; /* null terminated string */		\
 		seq->last_char = 0;	/* we have not come to the next header line */ \
-		if (seq->seq.l != seq->qual.l) return -2; /* qual string is shorter than seq string */ \
+		if (seq->seq.l != seq->qual.l) return -2; /* error: qual string is of a different length */ \
-		return seq->seq.l;												\
+		return seq->seq.l; \
 	}
 
 #define __KSEQ_TYPE(type_t)						\
@@ -199,10 +217,19 @@ typedef struct __kstring_t {
 		kstream_t *f;							\
 	} kseq_t;
 
-#define KSEQ_INIT(type_t, __read)				\
+#define KSEQ_INIT2(SCOPE, type_t, __read)		\
-	KSTREAM_INIT(type_t, __read, 4096)			\
+	KSTREAM_INIT(type_t, __read, 16384)			\
 	__KSEQ_TYPE(type_t)							\
-	__KSEQ_BASIC(type_t)						\
+	__KSEQ_BASIC(SCOPE, type_t)					\
-	__KSEQ_READ
+	__KSEQ_READ(SCOPE)
+
+#define KSEQ_INIT(type_t, __read) KSEQ_INIT2(static, type_t, __read)
+
+#define KSEQ_DECLARE(type_t) \
+	__KS_TYPE(type_t) \
+	__KSEQ_TYPE(type_t) \
+	extern kseq_t *kseq_init(type_t fd); \
+	void kseq_destroy(kseq_t *ks); \
+	int kseq_read(kseq_t *seq);
 
 #endif

simple_dp.c

@@ -8,7 +8,7 @@
 #include "utils.h"
 
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_DECLARE(gzFile)
 
 typedef struct {
 	int l;