Merge branch 'master' into master_fixes

Merged to master version b621d3a

Conflicts:
	Makefile
	bntseq.c
	bwa.c
	bwase.c
	bwaseqio.c
	bwtaln.c
	bwtindex.c
	bwtio.c
	bwtmisc.c
	bwtsw2_aux.c
	cs2nt.c
	fastmap.c
	khash.h
	kseq.h
	ksw.c
	kvec.h
	simple_dp.c
	utils.c
	utils.h

.gitignore

@@ -1,4 +1,5 @@
 *.[oa]
 bwa
 test
+test64
 .*.swp

Makefile

@@ -1,14 +1,11 @@
 CC=			gcc
-CXX=		g++
-CFLAGS=		-g -Wall -O2
+CFLAGS=		-g -Wall -O2 -msse2
 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 stdaln.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 \
+LOBJS=		utils.o kstring.o ksw.o kopen.o bwt.o bntseq.o bwa.o bwamem.o bwamem_pair.o
+AOBJS=		QSufSort.o bwt_gen.o stdaln.o bwase.o bwaseqio.o bwtgap.o bwtaln.o bamlite.o \
+			is.o bwtindex.o bwape.o \
 			bwtsw2_core.o bwtsw2_main.o bwtsw2_aux.o bwt_lite.o \
 			bwtsw2_chain.o fastmap.o bwtsw2_pair.o
 PROG=		bwa
@@ -26,43 +23,48 @@ SUBDIRS=	.
 all:$(PROG)
 
 bwa:libbwa.a $(AOBJS) main.o
-		$(CC) $(CFLAGS) $(DFLAGS) $(AOBJS) main.o -o $@ -L. -lbwa $(LIBS)
+		$(CC) $(CFLAGS) $(DFLAGS) $(AOBJS) main.o -o $@ $(LIBS) -L. -lbwa
 
 libbwa.a:$(LOBJS)
 		$(AR) -csru $@ $(LOBJS)
 
 clean:
 		rm -f gmon.out *.o a.out $(PROG) *~ *.a
 
+depend:
+	( LC_ALL=C ; export LC_ALL; makedepend -Y -- $(CFLAGS) -- *.c )
+
+# DO NOT DELETE THIS LINE -- make depend depends on it.
+
 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
+bamlite.o: utils.h bamlite.h
+bntseq.o: bntseq.h main.h utils.h kseq.h
+bwa.o: bntseq.h bwa.h bwt.h ksw.h utils.h kseq.h
+bwamem.o: kstring.h utils.h bwamem.h bwt.h bntseq.h bwa.h ksw.h kvec.h
+bwamem.o: ksort.h kbtree.h
+bwamem_pair.o: kstring.h utils.h bwamem.h bwt.h bntseq.h bwa.h kvec.h ksw.h
+bwape.o: bwtaln.h bwt.h stdaln.h kvec.h bntseq.h utils.h bwase.h bwa.h
+bwape.o: khash.h
+bwase.o: stdaln.h bwase.h bntseq.h bwt.h bwtaln.h utils.h kstring.h bwa.h
+bwaseqio.o: bwtaln.h bwt.h stdaln.h utils.h bamlite.h kseq.h
+bwt.o: utils.h bwt.h kvec.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
+bwtaln.o: bwtaln.h bwt.h stdaln.h bwtgap.h utils.h bwa.h bntseq.h
+bwtgap.o: bwtgap.h bwt.h bwtaln.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
+bwtsw2_aux.o: bntseq.h bwt_lite.h utils.h bwtsw2.h bwt.h stdaln.h kstring.h
+bwtsw2_aux.o: bwa.h kseq.h ksort.h
+bwtsw2_chain.o: bwtsw2.h bntseq.h bwt_lite.h bwt.h utils.h ksort.h
+bwtsw2_core.o: bwt_lite.h bwtsw2.h bntseq.h bwt.h kvec.h utils.h khash.h
+bwtsw2_core.o: ksort.h
+bwtsw2_main.o: bwt.h bwtsw2.h bntseq.h bwt_lite.h utils.h bwa.h
+bwtsw2_pair.o: utils.h bwt.h bntseq.h bwtsw2.h bwt_lite.h kstring.h ksw.h
+fastmap.o: bwa.h bntseq.h bwt.h bwamem.h kvec.h utils.h kseq.h
 is.o: utils.h
+kopen.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
+utils.o: utils.h ksort.h kseq.h

QSufSort.c

@@ -59,12 +59,9 @@ void QSufSortSuffixSort(qsint_t* __restrict V, qsint_t* __restrict I, const qsin
 	qsint_t i, j;
 	qsint_t s, negatedSortedGroupLength;
 	qsint_t numSymbolAggregated;
-	qsint_t maxNumInputSymbol;
 	qsint_t numSortedPos = 1;
 	qsint_t newAlphabetSize;
 
-	maxNumInputSymbol = largestInputSymbol - smallestInputSymbol + 1;
-
 	if (!skipTransform) {
 		/* bucketing possible*/
 		newAlphabetSize = QSufSortTransform(V, I, numChar, largestInputSymbol, smallestInputSymbol, 

bntseq.c

@@ -36,7 +36,7 @@
 #include "utils.h"
 
 #include "kseq.h"
-KSEQ_INIT(gzFile, err_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, 
@@ -310,21 +310,26 @@ int bwa_fa2pac(int argc, char *argv[])
 	return 0;
 }
 
+int bns_pos2rid(const bntseq_t *bns, int64_t pos_f)
+{
+	int left, mid, right;
+	if (pos_f >= bns->l_pac) return -1;
+	left = 0; mid = 0; right = bns->n_seqs;
+	while (left < right) { // binary search
+		mid = (left + right) >> 1;
+		if (pos_f >= bns->anns[mid].offset) {
+			if (mid == bns->n_seqs - 1) break;
+			if (pos_f < bns->anns[mid+1].offset) break; // bracketed
+			left = mid + 1;
+		} else right = mid;
+	}
+	return mid;
+}
+
 int bns_cnt_ambi(const bntseq_t *bns, int64_t pos_f, int len, int *ref_id)
 {
 	int left, mid, right, nn;
-	if (ref_id) {
-		left = 0; mid = 0; right = bns->n_seqs;
-		while (left < right) {
-			mid = (left + right) >> 1;
-			if (pos_f >= bns->anns[mid].offset) {
-				if (mid == bns->n_seqs - 1) break;
-				if (pos_f < bns->anns[mid+1].offset) break; // bracketed
-				left = mid + 1;
-			} else right = mid;
-		}
-		*ref_id = mid;
-	}
+	if (ref_id) *ref_id = bns_pos2rid(bns, pos_f);
 	left = 0; right = bns->n_holes; nn = 0;
 	while (left < right) {
 		mid = (left + right) >> 1;
@@ -343,3 +348,26 @@ int bns_cnt_ambi(const bntseq_t *bns, int64_t pos_f, int len, int *ref_id)
 	}
 	return nn;
 }
+
+uint8_t *bns_get_seq(int64_t l_pac, const uint8_t *pac, int64_t beg, int64_t end, int64_t *len)
+{
+	uint8_t *seq = 0;
+	if (end < beg) end ^= beg, beg ^= end, end ^= beg; // if end is smaller, swap
+	if (end > l_pac<<1) end = l_pac<<1;
+	if (beg < 0) beg = 0;
+	if (beg >= l_pac || end <= l_pac) {
+		int64_t k, l = 0;
+		*len = end - beg;
+		seq = xmalloc(end - beg);
+		if (beg >= l_pac) { // reverse strand
+			int64_t beg_f = (l_pac<<1) - 1 - end;
+			int64_t end_f = (l_pac<<1) - 1 - beg;
+			for (k = end_f; k > beg_f; --k)
+				seq[l++] = 3 - _get_pac(pac, k);
+		} else { // forward strand
+			for (k = beg; k < end; ++k)
+				seq[l++] = _get_pac(pac, k);
+		}
+	} else *len = 0; // if bridging the forward-reverse boundary, return nothing
+	return seq;
+}

bntseq.h

@@ -29,6 +29,7 @@
 #define BWT_BNTSEQ_H
 
 #include <stdint.h>
+#include <stdio.h>
 #include <zlib.h>
 
 #ifndef BWA_UBYTE
@@ -71,7 +72,9 @@ extern "C" {
 	bntseq_t *bns_restore_core(const char *ann_filename, const char* amb_filename, const char* pac_filename);
 	void bns_destroy(bntseq_t *bns);
 	int64_t bns_fasta2bntseq(gzFile fp_fa, const char *prefix, int for_only);
+	int bns_pos2rid(const bntseq_t *bns, int64_t pos_f);
 	int bns_cnt_ambi(const bntseq_t *bns, int64_t pos_f, int len, int *ref_id);
+	uint8_t *bns_get_seq(int64_t l_pac, const uint8_t *pac, int64_t beg, int64_t end, int64_t *len);
 
 #ifdef __cplusplus
 }