From: Heng Li <lh3@me.com>
Date: Fri, 2 Mar 2012 20:52:44 +0000 (-0500)
Subject: find the start positions; fixed a bug in example
X-Git-Tag: spawn-final~69
X-Git-Url: http://www.kaiwu.me/postgresql/commit/static/gitweb.js?a=commitdiff_plain;h=509017b53a5adb1cc58f3051c9720eb67ed79099;p=klib.git

find the start positions; fixed a bug in example
---

diff --git a/ksw.c b/ksw.c
index c2b5f9c..b7436ff 100644
--- a/ksw.c
+++ b/ksw.c
@@ -200,7 +200,7 @@ end_loop16:
 	return a->score + q->shift >= 255? 255 : a->score;
 }
 
-int ksw_sse2_8(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a) // the first gap costs -(_o+_e)
+int ksw_sse2_8(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a, int cutsc) // the first gap costs -(_o+_e)
 {
 	int slen, i, m_b, n_b, te = -1, gmax = 0;
 	uint64_t *b;
@@ -272,6 +272,7 @@ end_loop8:
 			gmax = imax; te = i;
 			for (j = 0; LIKELY(j < slen); ++j)
 				_mm_store_si128(Hmax + j, _mm_load_si128(H1 + j));
+			if (gmax == cutsc) break;
 		}
 		S = H1; H1 = H0; H0 = S;
 	}
@@ -295,8 +296,37 @@ end_loop8:
 
 int ksw_sse2(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a)
 {
+	a->qb = a->tb = -1;
 	if (q->size == 1) return ksw_sse2_16(q, tlen, target, a);
-	else return ksw_sse2_8(q, tlen, target, a);
+	else return ksw_sse2_8(q, tlen, target, a, 0);
+}
+
+static void revseq(int l, uint8_t *s)
+{
+	int i, t;
+	for (i = 0; i < l>>1; ++i)
+		t = s[i], s[i] = s[l - 1 - i], s[l - 1 - i] = t;
+}
+
+int ksw_align_short(int qlen, uint8_t *query, int tlen, uint8_t *target, int m, const int8_t *mat, ksw_aux_t *a)
+{
+	ksw_query_t *q;
+	ksw_aux_t alt;
+
+	q = ksw_qinit(2, qlen, query, m, mat);
+	ksw_sse2_8(q, tlen, target, a, 0);
+	free(q);
+	if (a->score == 0) return 0;
+	alt = *a;
+	revseq(a->qe + 1, query); revseq(a->te + 1, target); // +1 because qe/te points to the exact end, not the position after the end
+	q = ksw_qinit(2, a->qe + 1, query, m, mat);
+	ksw_sse2_8(q, a->te + 1, target, &alt, a->score);
+	revseq(a->qe + 1, query); revseq(a->te + 1, target);
+	free(q);
+	if (alt.score == a->score)
+		a->tb = a->te - alt.te, a->qb = a->qe - alt.qe;
+	else a->tb = a->qb = -1;
+	return a->score;
 }
 
 /*******************************************
@@ -355,7 +385,7 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 	// initialize scoring matrix
-	for (i = k = 0; i < 5; ++i) {
+	for (i = k = 0; i < 4; ++i) {
 		for (j = 0; j < 4; ++j)
 			mat[k++] = i == j? sa : -sb;
 		mat[k++] = 0; // ambiguous base
diff --git a/ksw.h b/ksw.h
index d93d6a9..ce8945e 100644
--- a/ksw.h
+++ b/ksw.h
@@ -10,6 +10,7 @@ typedef struct {
 	unsigned T; // threshold
 	// output
 	int score, te, qe, score2, te2;
+	int tb, qb; // tb and qb are only generated when calling ksw_align_16()
 } ksw_aux_t;
 
 #ifdef __cplusplus
@@ -39,7 +40,7 @@ extern "C" {
 	 *
 	 * @return        The maximum local score; if the returned value equals 255, the SW may not be finished
 	 */
-	int ksw_sse2_8(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a);
+	int ksw_sse2_8(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a, int cutsc);
 
 	/** Compute the maximum local score for queries initialized with ksw_qinit(2, ...) */
 	int ksw_sse2_16(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a);
@@ -47,6 +48,8 @@ extern "C" {
 	/** Unified interface for ksw_sse2_8() and ksw_sse2_16() */
 	int ksw_sse2(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a);
 
+	int ksw_align_short(int qlen, uint8_t *query, int tlen, uint8_t *target, int m, const int8_t *mat, ksw_aux_t *a);
+
 #ifdef __cplusplus
 }
 #endif