From c96f16e8103783f6b98d9f844ea3a7f2195e4834 Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Thu, 30 May 2024 04:49:40 +0200 Subject: Adding debian version 1.9.4-1. Signed-off-by: Daniel Baumann --- debian/vendor-h2o/deps/klib/ksw.c | 633 -------------------------------------- 1 file changed, 633 deletions(-) delete mode 100644 debian/vendor-h2o/deps/klib/ksw.c (limited to 'debian/vendor-h2o/deps/klib/ksw.c') diff --git a/debian/vendor-h2o/deps/klib/ksw.c b/debian/vendor-h2o/deps/klib/ksw.c deleted file mode 100644 index 742fec9..0000000 --- a/debian/vendor-h2o/deps/klib/ksw.c +++ /dev/null @@ -1,633 +0,0 @@ -/* The MIT License - - Copyright (c) 2011 by Attractive Chaos - - Permission is hereby granted, free of charge, to any person obtaining - a copy of this software and associated documentation files (the - "Software"), to deal in the Software without restriction, including - without limitation the rights to use, copy, modify, merge, publish, - distribute, sublicense, and/or sell copies of the Software, and to - permit persons to whom the Software is furnished to do so, subject to - the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS - BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN - ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN - CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -#include -#include -#include -#include "ksw.h" - -#ifdef __GNUC__ -#define LIKELY(x) __builtin_expect((x),1) -#define UNLIKELY(x) __builtin_expect((x),0) -#else -#define LIKELY(x) (x) -#define UNLIKELY(x) (x) -#endif - -const kswr_t g_defr = { 0, -1, -1, -1, -1, -1, -1 }; - -struct _kswq_t { - int qlen, slen; - uint8_t shift, mdiff, max, size; - __m128i *qp, *H0, *H1, *E, *Hmax; -}; - -/** - * Initialize the query data structure - * - * @param size Number of bytes used to store a score; valid valures are 1 or 2 - * @param qlen Length of the query sequence - * @param query Query sequence - * @param m Size of the alphabet - * @param mat Scoring matrix in a one-dimension array - * - * @return Query data structure - */ -kswq_t *ksw_qinit(int size, int qlen, const uint8_t *query, int m, const int8_t *mat) -{ - kswq_t *q; - int slen, a, tmp, p; - - size = size > 1? 2 : 1; - p = 8 * (3 - size); // # values per __m128i - slen = (qlen + p - 1) / p; // segmented length - q = (kswq_t*)malloc(sizeof(kswq_t) + 256 + 16 * slen * (m + 4)); // a single block of memory - q->qp = (__m128i*)(((size_t)q + sizeof(kswq_t) + 15) >> 4 << 4); // align memory - q->H0 = q->qp + slen * m; - q->H1 = q->H0 + slen; - q->E = q->H1 + slen; - q->Hmax = q->E + slen; - q->slen = slen; q->qlen = qlen; q->size = size; - // compute shift - tmp = m * m; - for (a = 0, q->shift = 127, q->mdiff = 0; a < tmp; ++a) { // find the minimum and maximum score - if (mat[a] < (int8_t)q->shift) q->shift = mat[a]; - if (mat[a] > (int8_t)q->mdiff) q->mdiff = mat[a]; - } - q->max = q->mdiff; - q->shift = 256 - q->shift; // NB: q->shift is uint8_t - q->mdiff += q->shift; // this is the difference between the min and max scores - // An example: p=8, qlen=19, slen=3 and segmentation: - // {{0,3,6,9,12,15,18,-1},{1,4,7,10,13,16,-1,-1},{2,5,8,11,14,17,-1,-1}} - if (size == 1) { - int8_t *t = (int8_t*)q->qp; - for (a = 0; a < m; ++a) { - int i, k, nlen = slen * p; - const int8_t *ma = mat + a * m; - for (i = 0; i < slen; ++i) - for (k = i; k < nlen; k += slen) // p iterations - *t++ = (k >= qlen? 0 : ma[query[k]]) + q->shift; - } - } else { - int16_t *t = (int16_t*)q->qp; - for (a = 0; a < m; ++a) { - int i, k, nlen = slen * p; - const int8_t *ma = mat + a * m; - for (i = 0; i < slen; ++i) - for (k = i; k < nlen; k += slen) // p iterations - *t++ = (k >= qlen? 0 : ma[query[k]]); - } - } - return q; -} - -kswr_t ksw_u8(kswq_t *q, int tlen, const uint8_t *target, int _gapo, int _gape, int xtra) // the first gap costs -(_o+_e) -{ - int slen, i, m_b, n_b, te = -1, gmax = 0, minsc, endsc; - uint64_t *b; - __m128i zero, gapoe, gape, shift, *H0, *H1, *E, *Hmax; - kswr_t r; - -#define __max_16(ret, xx) do { \ - (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 8)); \ - (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 4)); \ - (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 2)); \ - (xx) = _mm_max_epu8((xx), _mm_srli_si128((xx), 1)); \ - (ret) = _mm_extract_epi16((xx), 0) & 0x00ff; \ - } while (0) - - // initialization - r = g_defr; - minsc = (xtra&KSW_XSUBO)? xtra&0xffff : 0x10000; - endsc = (xtra&KSW_XSTOP)? xtra&0xffff : 0x10000; - m_b = n_b = 0; b = 0; - zero = _mm_set1_epi32(0); - gapoe = _mm_set1_epi8(_gapo + _gape); - gape = _mm_set1_epi8(_gape); - shift = _mm_set1_epi8(q->shift); - H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax; - slen = q->slen; - for (i = 0; i < slen; ++i) { - _mm_store_si128(E + i, zero); - _mm_store_si128(H0 + i, zero); - _mm_store_si128(Hmax + i, zero); - } - // the core loop - for (i = 0; i < tlen; ++i) { - int j, k, cmp, imax; - __m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector - h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example - h = _mm_slli_si128(h, 1); // h=H(i-1,-1); << instead of >> because x64 is little-endian - for (j = 0; LIKELY(j < slen); ++j) { - /* SW cells are computed in the following order: - * H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)} - * E(i+1,j) = max{H(i,j)-q, E(i,j)-r} - * F(i,j+1) = max{H(i,j)-q, F(i,j)-r} - */ - // compute H'(i,j); note that at the beginning, h=H'(i-1,j-1) - h = _mm_adds_epu8(h, _mm_load_si128(S + j)); - h = _mm_subs_epu8(h, shift); // h=H'(i-1,j-1)+S(i,j) - e = _mm_load_si128(E + j); // e=E'(i,j) - h = _mm_max_epu8(h, e); - h = _mm_max_epu8(h, f); // h=H'(i,j) - max = _mm_max_epu8(max, h); // set max - _mm_store_si128(H1 + j, h); // save to H'(i,j) - // now compute E'(i+1,j) - h = _mm_subs_epu8(h, gapoe); // h=H'(i,j)-gapo - e = _mm_subs_epu8(e, gape); // e=E'(i,j)-gape - e = _mm_max_epu8(e, h); // e=E'(i+1,j) - _mm_store_si128(E + j, e); // save to E'(i+1,j) - // now compute F'(i,j+1) - f = _mm_subs_epu8(f, gape); - f = _mm_max_epu8(f, h); - // get H'(i-1,j) and prepare for the next j - h = _mm_load_si128(H0 + j); // h=H'(i-1,j) - } - // NB: we do not need to set E(i,j) as we disallow adjecent insertion and then deletion - for (k = 0; LIKELY(k < 16); ++k) { // this block mimics SWPS3; NB: H(i,j) updated in the lazy-F loop cannot exceed max - f = _mm_slli_si128(f, 1); - for (j = 0; LIKELY(j < slen); ++j) { - h = _mm_load_si128(H1 + j); - h = _mm_max_epu8(h, f); // h=H'(i,j) - _mm_store_si128(H1 + j, h); - h = _mm_subs_epu8(h, gapoe); - f = _mm_subs_epu8(f, gape); - cmp = _mm_movemask_epi8(_mm_cmpeq_epi8(_mm_subs_epu8(f, h), zero)); - if (UNLIKELY(cmp == 0xffff)) goto end_loop16; - } - } -end_loop16: - //int k;for (k=0;k<16;++k)printf("%d ", ((uint8_t*)&max)[k]);printf("\n"); - __max_16(imax, max); // imax is the maximum number in max - if (imax >= minsc) { // write the b array; this condition adds branching unfornately - if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { // then append - if (n_b == m_b) { - m_b = m_b? m_b<<1 : 8; - b = (uint64_t*)realloc(b, 8 * m_b); - } - 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 - } - if (imax > gmax) { - gmax = imax; te = i; // te is the end position on the target - for (j = 0; LIKELY(j < slen); ++j) // keep the H1 vector - _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j)); - if (gmax + q->shift >= 255 || gmax >= endsc) break; - } - S = H1; H1 = H0; H0 = S; // swap H0 and H1 - } - r.score = gmax + q->shift < 255? gmax : 255; - r.te = te; - if (r.score != 255) { // get a->qe, the end of query match; find the 2nd best score - int max = -1, low, high, qlen = slen * 16; - uint8_t *t = (uint8_t*)Hmax; - for (i = 0; i < qlen; ++i, ++t) - if ((int)*t > max) max = *t, r.qe = i / 16 + i % 16 * slen; - //printf("%d,%d\n", max, gmax); - if (b) { - i = (r.score + q->max - 1) / q->max; - low = te - i; high = te + i; - for (i = 0; i < n_b; ++i) { - int e = (int32_t)b[i]; - if ((e < low || e > high) && (int)(b[i]>>32) > r.score2) - r.score2 = b[i]>>32, r.te2 = e; - } - } - } - free(b); - return r; -} - -kswr_t ksw_i16(kswq_t *q, int tlen, const uint8_t *target, int _gapo, int _gape, int xtra) // the first gap costs -(_o+_e) -{ - int slen, i, m_b, n_b, te = -1, gmax = 0, minsc, endsc; - uint64_t *b; - __m128i zero, gapoe, gape, *H0, *H1, *E, *Hmax; - kswr_t r; - -#define __max_8(ret, xx) do { \ - (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 8)); \ - (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 4)); \ - (xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 2)); \ - (ret) = _mm_extract_epi16((xx), 0); \ - } while (0) - - // initialization - r = g_defr; - minsc = (xtra&KSW_XSUBO)? xtra&0xffff : 0x10000; - endsc = (xtra&KSW_XSTOP)? xtra&0xffff : 0x10000; - m_b = n_b = 0; b = 0; - zero = _mm_set1_epi32(0); - gapoe = _mm_set1_epi16(_gapo + _gape); - gape = _mm_set1_epi16(_gape); - H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax; - slen = q->slen; - for (i = 0; i < slen; ++i) { - _mm_store_si128(E + i, zero); - _mm_store_si128(H0 + i, zero); - _mm_store_si128(Hmax + i, zero); - } - // the core loop - for (i = 0; i < tlen; ++i) { - int j, k, imax; - __m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector - h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example - h = _mm_slli_si128(h, 2); - for (j = 0; LIKELY(j < slen); ++j) { - h = _mm_adds_epi16(h, *S++); - e = _mm_load_si128(E + j); - h = _mm_max_epi16(h, e); - h = _mm_max_epi16(h, f); - max = _mm_max_epi16(max, h); - _mm_store_si128(H1 + j, h); - h = _mm_subs_epu16(h, gapoe); - e = _mm_subs_epu16(e, gape); - e = _mm_max_epi16(e, h); - _mm_store_si128(E + j, e); - f = _mm_subs_epu16(f, gape); - f = _mm_max_epi16(f, h); - h = _mm_load_si128(H0 + j); - } - for (k = 0; LIKELY(k < 16); ++k) { - f = _mm_slli_si128(f, 2); - for (j = 0; LIKELY(j < slen); ++j) { - h = _mm_load_si128(H1 + j); - h = _mm_max_epi16(h, f); - _mm_store_si128(H1 + j, h); - h = _mm_subs_epu16(h, gapoe); - f = _mm_subs_epu16(f, gape); - if(UNLIKELY(!_mm_movemask_epi8(_mm_cmpgt_epi16(f, h)))) goto end_loop8; - } - } -end_loop8: - __max_8(imax, max); - if (imax >= minsc) { - if (n_b == 0 || (int32_t)b[n_b-1] + 1 != i) { - if (n_b == m_b) { - m_b = m_b? m_b<<1 : 8; - b = (uint64_t*)realloc(b, 8 * m_b); - } - 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 - } - if (imax > gmax) { - gmax = imax; te = i; - for (j = 0; LIKELY(j < slen); ++j) - _mm_store_si128(Hmax + j, _mm_load_si128(H1 + j)); - if (gmax >= endsc) break; - } - S = H1; H1 = H0; H0 = S; - } - r.score = gmax; r.te = te; - { - int max = -1, low, high, qlen = slen * 8; - uint16_t *t = (uint16_t*)Hmax; - for (i = 0, r.qe = -1; i < qlen; ++i, ++t) - if ((int)*t > max) max = *t, r.qe = i / 8 + i % 8 * slen; - if (b) { - i = (r.score + q->max - 1) / q->max; - low = te - i; high = te + i; - for (i = 0; i < n_b; ++i) { - int e = (int32_t)b[i]; - if ((e < low || e > high) && (int)(b[i]>>32) > r.score2) - r.score2 = b[i]>>32, r.te2 = e; - } - } - } - free(b); - return r; -} - -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; -} - -kswr_t ksw_align(int qlen, uint8_t *query, int tlen, uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int xtra, kswq_t **qry) -{ - int size; - kswq_t *q; - kswr_t r, rr; - kswr_t (*func)(kswq_t*, int, const uint8_t*, int, int, int); - - q = (qry && *qry)? *qry : ksw_qinit((xtra&KSW_XBYTE)? 1 : 2, qlen, query, m, mat); - if (qry && *qry == 0) *qry = q; - func = q->size == 2? ksw_i16 : ksw_u8; - size = q->size; - r = func(q, tlen, target, gapo, gape, xtra); - if (qry == 0) free(q); - if ((xtra&KSW_XSTART) == 0 || ((xtra&KSW_XSUBO) && r.score < (xtra&0xffff))) return r; - revseq(r.qe + 1, query); revseq(r.te + 1, target); // +1 because qe/te points to the exact end, not the position after the end - q = ksw_qinit(size, r.qe + 1, query, m, mat); - rr = func(q, tlen, target, gapo, gape, KSW_XSTOP | r.score); - revseq(r.qe + 1, query); revseq(r.te + 1, target); - free(q); - if (r.score == rr.score) - r.tb = r.te - rr.te, r.qb = r.qe - rr.qe; - return r; -} - -/******************** - *** SW extension *** - ********************/ - -typedef struct { - int32_t h, e; -} eh_t; - -int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, int *_qle, int *_tle) -{ - eh_t *eh; // score array - int8_t *qp; // query profile - int i, j, k, gapoe = gapo + gape, beg, end, max, max_i, max_j, max_gap; - if (h0 < 0) h0 = 0; - // allocate memory - qp = malloc(qlen * m); - eh = calloc(qlen + 1, 8); - // generate the query profile - for (k = i = 0; k < m; ++k) { - const int8_t *p = &mat[k * m]; - for (j = 0; j < qlen; ++j) qp[i++] = p[query[j]]; - } - // fill the first row - eh[0].h = h0; eh[1].h = h0 > gapoe? h0 - gapoe : 0; - for (j = 2; j <= qlen && eh[j-1].h > gape; ++j) - eh[j].h = eh[j-1].h - gape; - // adjust $w if it is too large - k = m * m; - for (i = 0, max = 0; i < k; ++i) // get the max score - max = max > mat[i]? max : mat[i]; - max_gap = (int)((double)(qlen * max - gapo) / gape + 1.); - max_gap = max_gap > 1? max_gap : 1; - w = w < max_gap? w : max_gap; - // DP loop - max = h0, max_i = max_j = -1; - beg = 0, end = qlen; - for (i = 0; LIKELY(i < tlen); ++i) { - int f = 0, h1, m = 0, mj = -1; - int8_t *q = &qp[target[i] * qlen]; - // compute the first column - h1 = h0 - (gapo + gape * (i + 1)); - if (h1 < 0) h1 = 0; - // apply the band and the constraint (if provided) - if (beg < i - w) beg = i - w; - if (end > i + w + 1) end = i + w + 1; - if (end > qlen) end = qlen; - for (j = beg; LIKELY(j < end); ++j) { - // At the beginning of the loop: eh[j] = { H(i-1,j-1), E(i,j) }, f = F(i,j) and h1 = H(i,j-1) - // Similar to SSE2-SW, cells are computed in the following order: - // H(i,j) = max{H(i-1,j-1)+S(i,j), E(i,j), F(i,j)} - // E(i+1,j) = max{H(i,j)-gapo, E(i,j)} - gape - // F(i,j+1) = max{H(i,j)-gapo, F(i,j)} - gape - eh_t *p = &eh[j]; - int h = p->h, e = p->e; // get H(i-1,j-1) and E(i-1,j) - p->h = h1; // set H(i,j-1) for the next row - h += q[j]; - h = h > e? h : e; - h = h > f? h : f; - h1 = h; // save H(i,j) to h1 for the next column - mj = m > h? mj : j; - m = m > h? m : h; // m is stored at eh[mj+1] - h -= gapoe; - h = h > 0? h : 0; - e -= gape; - e = e > h? e : h; // computed E(i+1,j) - p->e = e; // save E(i+1,j) for the next row - f -= gape; - f = f > h? f : h; // computed F(i,j+1) - } - eh[end].h = h1; eh[end].e = 0; - if (m == 0) break; - if (m > max) max = m, max_i = i, max_j = mj; - // update beg and end for the next round - for (j = mj; j >= beg && eh[j].h; --j); - beg = j + 1; - for (j = mj + 2; j <= end && eh[j].h; ++j); - end = j; - //beg = 0; end = qlen; // uncomment this line for debugging - } - free(eh); free(qp); - if (_qle) *_qle = max_j + 1; - if (_tle) *_tle = max_i + 1; - return max; -} - -/******************** - * Global alignment * - ********************/ - -#define MINUS_INF -0x40000000 - -static inline uint32_t *push_cigar(int *n_cigar, int *m_cigar, uint32_t *cigar, int op, int len) -{ - if (*n_cigar == 0 || op != (cigar[(*n_cigar) - 1]&0xf)) { - if (*n_cigar == *m_cigar) { - *m_cigar = *m_cigar? (*m_cigar)<<1 : 4; - cigar = realloc(cigar, (*m_cigar) << 2); - } - cigar[(*n_cigar)++] = len<<4 | op; - } else cigar[(*n_cigar)-1] += len<<4; - return cigar; -} - -int ksw_global(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int *n_cigar_, uint32_t **cigar_) -{ - eh_t *eh; - int8_t *qp; // query profile - int i, j, k, gapoe = gapo + gape, score, n_col; - uint8_t *z; // backtrack matrix; in each cell: f<<4|e<<2|h; in principle, we can halve the memory, but backtrack will be a little more complex - if (n_cigar_) *n_cigar_ = 0; - // allocate memory - n_col = qlen < 2*w+1? qlen : 2*w+1; // maximum #columns of the backtrack matrix - z = malloc(n_col * tlen); - qp = malloc(qlen * m); - eh = calloc(qlen + 1, 8); - // generate the query profile - for (k = i = 0; k < m; ++k) { - const int8_t *p = &mat[k * m]; - for (j = 0; j < qlen; ++j) qp[i++] = p[query[j]]; - } - // fill the first row - eh[0].h = 0; eh[0].e = MINUS_INF; - for (j = 1; j <= qlen && j <= w; ++j) - eh[j].h = -(gapo + gape * j), eh[j].e = MINUS_INF; - for (; j <= qlen; ++j) eh[j].h = eh[j].e = MINUS_INF; // everything is -inf outside the band - // DP loop - for (i = 0; LIKELY(i < tlen); ++i) { // target sequence is in the outer loop - int32_t f = MINUS_INF, h1, beg, end; - int8_t *q = &qp[target[i] * qlen]; - uint8_t *zi = &z[i * n_col]; - beg = i > w? i - w : 0; - end = i + w + 1 < qlen? i + w + 1 : qlen; // only loop through [beg,end) of the query sequence - h1 = beg == 0? -(gapo + gape * (i + 1)) : MINUS_INF; - for (j = beg; LIKELY(j < end); ++j) { - // This loop is organized in a similar way to ksw_extend() and ksw_sse2(), except: - // 1) not checking h>0; 2) recording direction for backtracking - eh_t *p = &eh[j]; - int32_t h = p->h, e = p->e; - uint8_t d; // direction - p->h = h1; - h += q[j]; - d = h > e? 0 : 1; - h = h > e? h : e; - d = h > f? d : 2; - h = h > f? h : f; - h1 = h; - h -= gapoe; - e -= gape; - d |= e > h? 1<<2 : 0; - e = e > h? e : h; - p->e = e; - f -= gape; - d |= f > h? 2<<4 : 0; // if we want to halve the memory, use one bit only, instead of two - f = f > h? f : h; - zi[j - beg] = d; // z[i,j] keeps h for the current cell and e/f for the next cell - } - eh[end].h = h1; eh[end].e = MINUS_INF; - } - score = eh[qlen].h; - if (n_cigar_ && cigar_) { // backtrack - int n_cigar = 0, m_cigar = 0, which = 0; - uint32_t *cigar = 0, tmp; - i = tlen - 1; k = (i + w + 1 < qlen? i + w + 1 : qlen) - 1; // (i,k) points to the last cell - while (i >= 0 && k >= 0) { - which = z[i * n_col + (k - (i > w? i - w : 0))] >> (which<<1) & 3; - if (which == 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 0, 1), --i, --k; - else if (which == 1) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 2, 1), --i; - else cigar = push_cigar(&n_cigar, &m_cigar, cigar, 1, 1), --k; - } - if (i >= 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 2, i + 1); - if (k >= 0) cigar = push_cigar(&n_cigar, &m_cigar, cigar, 1, k + 1); - for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR - tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp; - *n_cigar_ = n_cigar, *cigar_ = cigar; - } - free(eh); free(qp); free(z); - return score; -} - -/******************************************* - * Main function (not compiled by default) * - *******************************************/ - -#ifdef _KSW_MAIN - -#include -#include -#include -#include "kseq.h" -KSEQ_INIT(gzFile, gzread) - -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, - 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, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 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, 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, 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, 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, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 -}; - -int main(int argc, char *argv[]) -{ - int c, sa = 1, sb = 3, i, j, k, forward_only = 0, max_rseq = 0; - int8_t mat[25]; - int gapo = 5, gape = 2, minsc = 0, xtra = KSW_XSTART; - uint8_t *rseq = 0; - gzFile fpt, fpq; - kseq_t *kst, *ksq; - - // parse command line - while ((c = getopt(argc, argv, "a:b:q:r:ft:1")) >= 0) { - switch (c) { - case 'a': sa = atoi(optarg); break; - case 'b': sb = atoi(optarg); break; - case 'q': gapo = atoi(optarg); break; - case 'r': gape = atoi(optarg); break; - case 't': minsc = atoi(optarg); break; - case 'f': forward_only = 1; break; - case '1': xtra |= KSW_XBYTE; break; - } - } - if (optind + 2 > argc) { - fprintf(stderr, "Usage: ksw [-1] [-f] [-a%d] [-b%d] [-q%d] [-r%d] [-t%d] \n", sa, sb, gapo, gape, minsc); - return 1; - } - if (minsc > 0xffff) minsc = 0xffff; - xtra |= KSW_XSUBO | minsc; - // initialize scoring matrix - for (i = k = 0; i < 4; ++i) { - for (j = 0; j < 4; ++j) - mat[k++] = i == j? sa : -sb; - mat[k++] = 0; // ambiguous base - } - for (j = 0; j < 5; ++j) mat[k++] = 0; - // open file - fpt = gzopen(argv[optind], "r"); kst = kseq_init(fpt); - fpq = gzopen(argv[optind+1], "r"); ksq = kseq_init(fpq); - // all-pair alignment - while (kseq_read(ksq) > 0) { - kswq_t *q[2] = {0, 0}; - kswr_t r; - for (i = 0; i < (int)ksq->seq.l; ++i) ksq->seq.s[i] = seq_nt4_table[(int)ksq->seq.s[i]]; - if (!forward_only) { // reverse - if ((int)ksq->seq.m > max_rseq) { - max_rseq = ksq->seq.m; - rseq = (uint8_t*)realloc(rseq, max_rseq); - } - for (i = 0, j = ksq->seq.l - 1; i < (int)ksq->seq.l; ++i, --j) - rseq[j] = ksq->seq.s[i] == 4? 4 : 3 - ksq->seq.s[i]; - } - gzrewind(fpt); kseq_rewind(kst); - while (kseq_read(kst) > 0) { - for (i = 0; i < (int)kst->seq.l; ++i) kst->seq.s[i] = seq_nt4_table[(int)kst->seq.s[i]]; - r = ksw_align(ksq->seq.l, (uint8_t*)ksq->seq.s, kst->seq.l, (uint8_t*)kst->seq.s, 5, mat, gapo, gape, xtra, &q[0]); - if (r.score >= minsc) - printf("%s\t%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n", kst->name.s, r.tb, r.te+1, ksq->name.s, r.qb, r.qe+1, r.score, r.score2, r.te2); - if (rseq) { - r = ksw_align(ksq->seq.l, rseq, kst->seq.l, (uint8_t*)kst->seq.s, 5, mat, gapo, gape, xtra, &q[1]); - if (r.score >= minsc) - printf("%s\t%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n", kst->name.s, r.tb, r.te+1, ksq->name.s, (int)ksq->seq.l - r.qb, (int)ksq->seq.l - 1 - r.qe, r.score, r.score2, r.te2); - } - } - free(q[0]); free(q[1]); - } - free(rseq); - kseq_destroy(kst); gzclose(fpt); - kseq_destroy(ksq); gzclose(fpq); - return 0; -} -#endif -- cgit v1.2.3