From 182cb2e89ca57d1d916d7a967a19414e08751e25 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 19 Nov 2011 19:38:21 -0500 Subject: [PATCH] use standard SW when no SSE2 --- Makefile | 2 +- bwtsw2_pair.c | 80 +++++++++++++++++++++++++++++++++------------------ ksw.c | 4 ++- main.c | 2 +- stdaln.c | 2 +- 5 files changed, 58 insertions(+), 32 deletions(-) diff --git a/Makefile b/Makefile index 65ce6dd..1fe8462 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CC= gcc CXX= g++ CFLAGS= -g -Wall -O2 CXXFLAGS= $(CFLAGS) -DFLAGS= -DHAVE_PTHREAD #-D_FILE_OFFSET_BITS=64 +DFLAGS= -DHAVE_PTHREAD #-D_NO_SSE2 #-D_FILE_OFFSET_BITS=64 OBJS= QSufSort.o bwt_gen.o utils.o bwt.o bwtio.o bwtaln.o bwtgap.o \ is.o bntseq.o bwtmisc.o bwtindex.o ksw.o stdaln.o simple_dp.o \ bwaseqio.o bwase.o bwape.o kstring.o cs2nt.o \ diff --git a/bwtsw2_pair.c b/bwtsw2_pair.c index 581e19c..ad1bb3f 100644 --- a/bwtsw2_pair.c +++ b/bwtsw2_pair.c @@ -5,8 +5,12 @@ #include "bwt.h" #include "bntseq.h" #include "bwtsw2.h" -#include "ksw.h" #include "kstring.h" +#ifndef _NO_SSE2 +#include "ksw.h" +#else +#include "stdaln.h" +#endif #define MAX_INS 20000 #define MIN_RATIO 0.8 @@ -92,8 +96,6 @@ void bsw2_pair1(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, const b int64_t k, beg, end; uint8_t *seq, *ref; int i; - ksw_query_t *q; - ksw_aux_t aux[2]; // compute the region start and end a->n_seeds = 1; a->flag |= BSW2_FLAG_MATESW; // before calling this routine, *a has been cleared with memset(0); the flag is set with 1<<6/7 if (h->is_rev == 0) { @@ -123,32 +125,54 @@ void bsw2_pair1(const bsw2opt_t *opt, int64_t l_pac, const uint8_t *pac, const b for (i = 0; i < l_mseq; ++i) // on the forward strand seq[i] = nst_nt4_table[(int)mseq[i]]; } - // forward Smith-Waterman - aux[0].T = opt->t; aux[0].gapo = opt->q; aux[0].gape = opt->r; aux[1] = aux[0]; - q = ksw_qinit(l_mseq * g_mat[0] < 250? 1 : 2, l_mseq, seq, 5, g_mat); - ksw_sse2(q, end - beg, ref, &aux[0]); - free(q); - if (aux[0].score < opt->t) { - free(seq); - return; +#ifndef _NO_SSE2 + { + ksw_query_t *q; + ksw_aux_t aux[2]; + // forward Smith-Waterman + aux[0].T = opt->t; aux[0].gapo = opt->q; aux[0].gape = opt->r; aux[1] = aux[0]; + q = ksw_qinit(l_mseq * g_mat[0] < 250? 1 : 2, l_mseq, seq, 5, g_mat); + ksw_sse2(q, end - beg, ref, &aux[0]); + free(q); + if (aux[0].score < opt->t) { + free(seq); + return; + } + ++aux[0].qe; ++aux[0].te; + // reverse Smith-Waterman + seq_reverse(aux[0].qe, seq, 0); + seq_reverse(aux[0].te, ref, 0); + q = ksw_qinit(aux[0].qe * g_mat[0] < 250? 1 : 2, aux[0].qe, seq, 5, g_mat); + ksw_sse2(q, aux[0].te, ref, &aux[1]); + free(q); + ++aux[1].qe; ++aux[1].te; + // write output + a->G = aux[0].score; + a->G2 = aux[0].score2 > aux[1].score2? aux[0].score2 : aux[1].score2; + if (a->G2 < opt->t) a->G2 = 0; + if (a->G2) a->flag |= BSW2_FLAG_TANDEM; + a->k = beg + (aux[0].te - aux[1].te); + a->len = aux[1].te; + a->beg = aux[0].qe - aux[1].qe; + a->end = aux[0].qe; } - ++aux[0].qe; ++aux[0].te; - // reverse Smith-Waterman - seq_reverse(aux[0].qe, seq, 0); - seq_reverse(aux[0].te, ref, 0); - q = ksw_qinit(aux[0].qe * g_mat[0] < 250? 1 : 2, aux[0].qe, seq, 5, g_mat); - ksw_sse2(q, aux[0].te, ref, &aux[1]); - free(q); - ++aux[1].qe; ++aux[1].te; - // write output - a->G = aux[0].score; - a->G2 = aux[0].score2 > aux[1].score2? aux[0].score2 : aux[1].score2; - if (a->G2 < opt->t) a->G2 = 0; - if (a->G2) a->flag |= BSW2_FLAG_TANDEM; - a->k = beg + (aux[0].te - aux[1].te); - a->len = aux[1].te; - a->beg = aux[0].qe - aux[1].qe; - a->end = aux[0].qe; +#else + { + AlnParam ap; + path_t path[2]; + int matrix[25]; + for (i = 0; i < 25; ++i) matrix[i] = g_mat[i]; + ap.gap_open = opt->q; ap.gap_ext = opt->r; ap.gap_end = opt->r; + ap.matrix = matrix; ap.row = 5; ap.band_width = 50; + a->G = aln_local_core(ref, end - beg, seq, l_mseq, &ap, path, 0, opt->t, &a->G2); + if (a->G < opt->t) a->G = 0; + if (a->G2 < opt->t) a->G2 = 0; + a->k = beg + path[0].i - 1; + a->len = path[1].i - path[0].i + 1; + a->beg = path[0].j - 1; + a->end = path[1].j; + } +#endif if (a->is_rev) i = a->beg, a->beg = l_mseq - a->end, a->end = l_mseq - i; free(seq); } diff --git a/ksw.c b/ksw.c index c2b5f9c..bd29e96 100644 --- a/ksw.c +++ b/ksw.c @@ -23,6 +23,7 @@ SOFTWARE. */ +#ifndef _NO_SSE2 #include #include #include @@ -396,4 +397,5 @@ int main(int argc, char *argv[]) kseq_destroy(ksq); gzclose(fpq); return 0; } -#endif +#endif // _KSW_MAIN +#endif // _NO_SSE2 diff --git a/main.c b/main.c index 4d4402f..10e3314 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.6.0-r90-dev" +#define PACKAGE_VERSION "0.6.0-r91-dev" #endif static int usage() diff --git a/stdaln.c b/stdaln.c index 7b55b2e..eb41882 100644 --- a/stdaln.c +++ b/stdaln.c @@ -631,7 +631,7 @@ int aln_local_core(unsigned char *seq1, int len1, unsigned char *seq2, int len2, score_f += of_base; if (score_f < thres) { /* no matching residue at all, 090218 */ - *path_len = 0; + if (path_len) *path_len = 0; goto end_func; } if (path == 0) goto end_func; /* skip path-filling */