diff --git a/CHANGES.md b/CHANGES.md index 2d907665..f9b67f88 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -21,6 +21,9 @@ * Include [ZStr](https://github.com/mateidavid/zstr/) in our own repository instead of downloading it at build time. This should make it possible to build strobealign without internet access. +* #357: Fix some suboptimal alignment ends. Sometimes an end was soft-clipped + although a better alignment with an insertion or deletion existed that + extends to the end of the read. ## v0.12.0 (2023-11-23) diff --git a/CMakeLists.txt b/CMakeLists.txt index 988b1235..d1088d22 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -63,6 +63,7 @@ add_library(salib STATIC ${SOURCES} ext/xxhash.c ext/ssw/ssw_cpp.cpp ext/ssw/ssw.c + ext/ksw2/ksw2_extz2_sse.c ) target_include_directories(salib PUBLIC src/ ext/ ${PROJECT_BINARY_DIR}) target_link_libraries(salib PUBLIC ZLIB::ZLIB Threads::Threads zstr::zstr) diff --git a/ext/README.md b/ext/README.md index 55e874b7..59e3503d 100644 --- a/ext/README.md +++ b/ext/README.md @@ -60,3 +60,10 @@ License: See xxhash.c Homepage: https://github.com/mateidavid/zstr Commit used: 755da7890ea22478a702e3139092e6c964fab1f5 License: See zstr/LICENSE + + +## ksw2 + +https://github.com/lh3/ksw2 +https://raw.githubusercontent.com/lh3/ksw2/06b2183b0f6646d82f2e3f5884008a1b4582f5b5/ksw2.h +https://raw.githubusercontent.com/lh3/ksw2/06b2183b0f6646d82f2e3f5884008a1b4582f5b5/ksw2_extz2_sse.c diff --git a/ext/ksw2.h b/ext/ksw2.h deleted file mode 100644 index 04eeda27..00000000 --- a/ext/ksw2.h +++ /dev/null @@ -1,290 +0,0 @@ -#ifndef KSW2_H_ -#define KSW2_H_ - -#include - -#define KSW_NEG_INF -0x40000000 - -#define KSW_EZ_SCORE_ONLY 0x01 // don't record alignment path/cigar -#define KSW_EZ_RIGHT 0x02 // right-align gaps -#define KSW_EZ_GENERIC_SC 0x04 // without this flag: match/mismatch only; last symbol is a wildcard -#define KSW_EZ_APPROX_MAX 0x08 // approximate max; this is faster with sse -#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse -#define KSW_EZ_EXTZ_ONLY 0x40 // only perform extension -#define KSW_EZ_REV_CIGAR 0x80 // reverse CIGAR in the output -#define KSW_EZ_SPLICE_FOR 0x100 -#define KSW_EZ_SPLICE_REV 0x200 -#define KSW_EZ_SPLICE_FLANK 0x400 - -#ifdef __cplusplus -extern "C" { -#endif - -typedef struct { - uint32_t max: 31, zdropped: 1; - int max_q, max_t; // max extension coordinate - int mqe, mqe_t; // max score when reaching the end of query - int mte, mte_q; // max score when reaching the end of target - int score; // max score reaching both ends; may be KSW_NEG_INF - int m_cigar, n_cigar; - int reach_end; - uint32_t *cigar; -} ksw_extz_t; - -/** - * NW-like extension - * - * @param km memory pool, when used with kalloc - * @param qlen query length - * @param query query sequence with 0 <= query[i] < m - * @param tlen target length - * @param target target sequence with 0 <= target[i] < m - * @param m number of residue types - * @param mat m*m scoring mattrix in one-dimension array - * @param gapo gap open penalty; a gap of length l cost "-(gapo+l*gape)" - * @param gape gap extension penalty - * @param w band width (<0 to disable) - * @param zdrop off-diagonal drop-off to stop extension (positive; <0 to disable) - * @param flag flag (see KSW_EZ_* macros) - * @param ez (out) scores and cigar - */ -void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez); - -void ksw_extz2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t q, - int8_t e, - int w, - int zdrop, - int end_bonus, - int flag, - ksw_extz_t *ez); - -void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez); - -void ksw_extd2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t gapo, - int8_t gape, - int8_t gapo2, - int8_t gape2, - int w, - int zdrop, - int end_bonus, - int flag, - ksw_extz_t *ez); - -void ksw_exts2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t gapo, - int8_t gape, - int8_t gapo2, - int8_t noncan, - int zdrop, - int flag, - ksw_extz_t *ez); - -void ksw_extf2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t mch, - int8_t mis, - int8_t e, - int w, - int xdrop, - ksw_extz_t *ez); - -/** - * Global alignment - * - * (first 10 parameters identical to ksw_extz_sse()) - * @param m_cigar (modified) max CIGAR length; feed 0 if cigar==0 - * @param n_cigar (out) number of CIGAR elements - * @param cigar (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, ) - * - * @return score of the alignment - */ -int ksw_gg(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t gapo, - int8_t gape, - int w, - int *m_cigar_, - int *n_cigar_, - uint32_t **cigar_); -int ksw_gg2(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t gapo, - int8_t gape, - int w, - int *m_cigar_, - int *n_cigar_, - uint32_t **cigar_); -int ksw_gg2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t gapo, - int8_t gape, - int w, - int *m_cigar_, - int *n_cigar_, - uint32_t **cigar_); - -void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat); -int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te); - -#ifdef __cplusplus -} -#endif - -/************************************ - *** Private macros and functions *** - ************************************/ - -#ifdef HAVE_KALLOC -#include "kalloc.h" -#else -#include -#define kmalloc(km, size) malloc((size)) -#define kcalloc(km, count, size) calloc((count), (size)) -#define krealloc(km, ptr, size) realloc((ptr), (size)) -#define kfree(km, ptr) free((ptr)) -#endif - -static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t 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 = (uint32_t *) krealloc(km, cigar, (*m_cigar) << 2); - } - cigar[(*n_cigar)++] = len << 4 | op; - } else - cigar[(*n_cigar) - 1] += len << 4; - return cigar; -} - -// In the backtrack matrix, value p[] has the following structure: -// bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F} -// bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E}) -// bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F}) -static inline void ksw_backtrack(void *km, - int is_rot, - int is_rev, - int min_intron_len, - const uint8_t *p, - const int *off, - const int *off_end, - int n_col, - int i0, - int j0, - int *m_cigar_, - int *n_cigar_, - uint32_t **cigar_) { // p[] - lower 3 bits: which type gets the max; bit - int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0; - uint32_t *cigar = *cigar_, tmp; - while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check - int force_state = -1; - if (is_rot) { - r = i + j; - if (i < off[r]) - force_state = 2; - if (off_end && i > off_end[r]) - force_state = 1; - tmp = force_state < 0 ? p[(size_t) r * n_col + i - off[r]] : 0; - } else { - if (j < off[i]) - force_state = 2; - if (off_end && j > off_end[i]) - force_state = 1; - tmp = force_state < 0 ? p[(size_t) i * n_col + j - off[i]] : 0; - } - if (state == 0) - state = tmp & 7; // if requesting the H state, find state one maximizes it. - else if (!(tmp >> (state + 2) & 1)) - state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H - if (state == 0) - state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure - if (force_state >= 0) - state = force_state; - if (state == 0) - cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match - else if (state == 1 || (state == 3 && min_intron_len <= 0)) - cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion - else if (state == 3 && min_intron_len > 0) - cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron - else - cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion - } - if (i >= 0) - cigar = ksw_push_cigar(km, - &n_cigar, - &m_cigar, - cigar, - min_intron_len > 0 && i >= min_intron_len ? 3 : 2, - i + 1); // first deletion - if (j >= 0) - cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion - if (!is_rev) - 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; - *m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar; -} - -static inline void ksw_reset_extz(ksw_extz_t *ez) { - ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1; - ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF; - ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0; -} - -static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e) { - int r, t; - if (is_rot) - r = a, t = b; - else - r = a + b, t = a; - if (H > (int32_t) ez->max) { - ez->max = H, ez->max_t = t, ez->max_q = r - t; - } else if (t >= ez->max_t && r - t >= ez->max_q) { - int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l; - l = tl > ql ? tl - ql : ql - tl; - if (zdrop >= 0 && ez->max - H > zdrop + l * e) { - ez->zdropped = 1; - return 1; - } - } - return 0; -} -#endif \ No newline at end of file diff --git a/ext/ksw2/LICENSE.txt b/ext/ksw2/LICENSE.txt new file mode 100644 index 00000000..1a06f649 --- /dev/null +++ b/ext/ksw2/LICENSE.txt @@ -0,0 +1,24 @@ +The MIT License + +Copyright (c) 2018- Dana-Farber Cancer Institute + 2017-2018 Broad Institute, Inc. + +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. diff --git a/ext/ksw2/ksw2.h b/ext/ksw2/ksw2.h new file mode 100644 index 00000000..01edd45e --- /dev/null +++ b/ext/ksw2/ksw2.h @@ -0,0 +1,208 @@ +#ifndef KSW2_H_ +#define KSW2_H_ + +#include + +#define KSW_NEG_INF -0x40000000 + +#define KSW_EZ_SCORE_ONLY 0x01 // don't record alignment path/cigar +#define KSW_EZ_RIGHT 0x02 // right-align gaps +#define KSW_EZ_GENERIC_SC 0x04 // without this flag: match/mismatch only; last symbol is a wildcard +#define KSW_EZ_APPROX_MAX 0x08 // approximate max; this is faster with sse +#define KSW_EZ_APPROX_DROP 0x10 // approximate Z-drop; faster with sse +#define KSW_EZ_EXTZ_ONLY 0x40 // only perform extension +#define KSW_EZ_REV_CIGAR 0x80 // reverse CIGAR in the output +#define KSW_EZ_SPLICE_FOR 0x100 +#define KSW_EZ_SPLICE_REV 0x200 +#define KSW_EZ_SPLICE_FLANK 0x400 +#define KSW_EZ_EQX 0x800 + +// The subset of CIGAR operators used by ksw code. +// Use MM_CIGAR_* from minimap.h if you need the full list. +#define KSW_CIGAR_MATCH 0 +#define KSW_CIGAR_INS 1 +#define KSW_CIGAR_DEL 2 +#define KSW_CIGAR_N_SKIP 3 +#define KSW_CIGAR_EQ 7 +#define KSW_CIGAR_X 8 + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct { + uint32_t max:31, zdropped:1; + int max_q, max_t; // max extension coordinate + int mqe, mqe_t; // max score when reaching the end of query + int mte, mte_q; // max score when reaching the end of target + int score; // max score reaching both ends; may be KSW_NEG_INF + int m_cigar, n_cigar; + int reach_end; + uint32_t *cigar; +} ksw_extz_t; + +/** + * NW-like extension + * + * @param km memory pool, when used with kalloc + * @param qlen query length + * @param query query sequence with 0 <= query[i] < m + * @param tlen target length + * @param target target sequence with 0 <= target[i] < m + * @param m number of residue types + * @param mat m*m scoring mattrix in one-dimension array + * @param gapo gap open penalty; a gap of length l cost "-(gapo+l*gape)" + * @param gape gap extension penalty + * @param w band width (<0 to disable) + * @param zdrop off-diagonal drop-off to stop extension (positive; <0 to disable) + * @param flag flag (see KSW_EZ_* macros) + * @param ez (out) scores and cigar + */ +void ksw_extz(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez); + +void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + +void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez); + +void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + +void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, + int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez); + +void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez); + +/** + * Global alignment + * + * (first 10 parameters identical to ksw_extz_sse()) + * @param m_cigar (modified) max CIGAR length; feed 0 if cigar==0 + * @param n_cigar (out) number of CIGAR elements + * @param cigar (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, ) + * + * @return score of the alignment + */ +int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); +int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); +int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); + +void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat); +int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te); + +#ifdef __cplusplus +} +#endif + +/************************************ + *** Private macros and functions *** + ************************************/ + +#ifdef HAVE_KALLOC +#include "kalloc.h" +#else +#include +#define kmalloc(km, size) malloc((size)) +#define kcalloc(km, count, size) calloc((count), (size)) +#define krealloc(km, ptr, size) realloc((ptr), (size)) +#define kfree(km, ptr) free((ptr)) +#endif + +static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uint32_t *cigar, uint32_t 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 = (uint32_t*)krealloc(km, cigar, (*m_cigar) << 2); + } + cigar[(*n_cigar)++] = len<<4 | op; + } else cigar[(*n_cigar)-1] += len<<4; + return cigar; +} + +// In the backtrack matrix, value p[] has the following structure: +// bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F} +// bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E}) +// bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F}) +static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0, + int *m_cigar_, int *n_cigar_, uint32_t **cigar_) +{ // p[] - lower 3 bits: which type gets the max; bit + int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0; + uint32_t *cigar = *cigar_, tmp; + while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check + int force_state = -1; + if (is_rot) { + r = i + j; + if (i < off[r]) force_state = 2; + if (off_end && i > off_end[r]) force_state = 1; + tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0; + } else { + if (j < off[i]) force_state = 2; + if (off_end && j > off_end[i]) force_state = 1; + tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0; + } + if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it. + else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H + if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure + if (force_state >= 0) state = force_state; + if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_MATCH, 1), --i, --j; + else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_DEL, 1), --i; + else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_N_SKIP, 1), --i; + else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, 1), --j; + } + if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? KSW_CIGAR_N_SKIP : KSW_CIGAR_DEL, i + 1); // first deletion + if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, j + 1); // first insertion + if (!is_rev) + 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; + *m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar; +} + +static inline void ksw_cigar2eqx(void *km, const uint8_t *query, const uint8_t *target, int nc0, const uint32_t *ci0, int *mc1, int *nc1, uint32_t **ci1) +{ + int i, k, x = 0, y = 0; + *nc1 = 0; + for (k = 0; k < nc0; ++k) { + int op = ci0[k]&0xf, len = ci0[k]>>4; + if (op == KSW_CIGAR_MATCH) { + for (i = 0; i < len; ++i) { + if (target[x + i] == query[y + i]) ksw_push_cigar(km, nc1, mc1, *ci1, KSW_CIGAR_EQ, 1); + else ksw_push_cigar(km, nc1, mc1, *ci1, KSW_CIGAR_X, 1); + } + x += len, y += len; + } else { + ksw_push_cigar(km, nc1, mc1, *ci1, op, len); + if (op == KSW_CIGAR_DEL || op == KSW_CIGAR_N_SKIP) x += len; + else if (op == KSW_CIGAR_INS) y += len; + else if (op == KSW_CIGAR_EQ || op == KSW_CIGAR_X) x += len, y += len; + } + } +} + +static inline void ksw_reset_extz(ksw_extz_t *ez) +{ + ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1; + ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF; + ez->n_cigar = 0, ez->zdropped = 0, ez->reach_end = 0; +} + +static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a, int b, int zdrop, int8_t e) +{ + int r, t; + if (is_rot) r = a, t = b; + else r = a + b, t = a; + if (H > (int32_t)ez->max) { + ez->max = H, ez->max_t = t, ez->max_q = r - t; + } else if (t >= ez->max_t && r - t >= ez->max_q) { + int tl = t - ez->max_t, ql = (r - t) - ez->max_q, l; + l = tl > ql? tl - ql : ql - tl; + if (zdrop >= 0 && ez->max - H > zdrop + l * e) { + ez->zdropped = 1; + return 1; + } + } + return 0; +} +#endif diff --git a/ext/ksw2/ksw2_extz2_sse.c b/ext/ksw2/ksw2_extz2_sse.c new file mode 100644 index 00000000..02bb4c2a --- /dev/null +++ b/ext/ksw2/ksw2_extz2_sse.c @@ -0,0 +1,305 @@ +#include +#include +#include "ksw2.h" + +#ifdef __SSE2__ +#include + +#ifdef KSW_SSE2_ONLY +#undef __SSE4_1__ +#endif + +#ifdef __SSE4_1__ +#include +#endif + +#ifdef KSW_CPU_DISPATCH +#ifdef __SSE4_1__ +void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +#else +void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +#endif +#else +void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) +#endif // ~KSW_CPU_DISPATCH +{ +#define __dp_code_block1 \ + z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \ + xt1 = _mm_load_si128(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ + tmp = _mm_srli_si128(xt1, 15); /* tmp <- x[r-1][t+15] */ \ + xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \ + x1_ = tmp; \ + vt1 = _mm_load_si128(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ + tmp = _mm_srli_si128(vt1, 15); /* tmp <- v[r-1][t+15] */ \ + vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \ + v1_ = tmp; \ + a = _mm_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ + ut = _mm_load_si128(&u[t]); /* ut <- u[t..t+15] */ \ + b = _mm_add_epi8(_mm_load_si128(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ + +#define __dp_code_block2 \ + z = _mm_max_epu8(z, b); /* z = max(z, b); this works because both are non-negative */ \ + z = _mm_min_epu8(z, max_sc_); \ + _mm_store_si128(&u[t], _mm_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ + _mm_store_si128(&v[t], _mm_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ + z = _mm_sub_epi8(z, q_); \ + a = _mm_sub_epi8(a, z); \ + b = _mm_sub_epi8(b, z); + + int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc; + int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); + int32_t *H = 0, H0 = 0, last_H0_t = 0; + uint8_t *qr, *sf, *mem, *mem2 = 0; + __m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_; + __m128i *u, *v, *x, *y, *s, *p = 0; + + ksw_reset_extz(ez); + if (m <= 0 || qlen <= 0 || tlen <= 0) return; + + zero_ = _mm_set1_epi8(0); + q_ = _mm_set1_epi8(q); + qe2_ = _mm_set1_epi8((q + e) * 2); + flag1_ = _mm_set1_epi8(1); + flag2_ = _mm_set1_epi8(2); + flag8_ = _mm_set1_epi8(0x08); + flag16_ = _mm_set1_epi8(0x10); + sc_mch_ = _mm_set1_epi8(mat[0]); + sc_mis_ = _mm_set1_epi8(mat[1]); + sc_N_ = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]); + m1_ = _mm_set1_epi8(m - 1); // wildcard + max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2); + + if (w < 0) w = tlen > qlen? tlen : qlen; + wl = wr = w; + tlen_ = (tlen + 15) / 16; + n_col_ = qlen < tlen? qlen : tlen; + n_col_ = ((n_col_ < w + 1? n_col_ : w + 1) + 15) / 16 + 1; + qlen_ = (qlen + 15) / 16; + for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { + max_sc = max_sc > mat[t]? max_sc : mat[t]; + min_sc = min_sc < mat[t]? min_sc : mat[t]; + } + if (-min_sc > 2 * (q + e)) return; // otherwise, we won't see any mismatches + + mem = (uint8_t*)kcalloc(km, tlen_ * 6 + qlen_ + 1, 16); + u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned + v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16; + if (!approx_max) { + H = (int32_t*)kmalloc(km, tlen_ * 16 * 4); + for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF; + } + if (with_cigar) { + mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16); + p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4); + off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); + off_end = off + qlen + tlen - 1; + } + + for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t]; + memcpy(sf, target, tlen); + + for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { + int st = 0, en = tlen - 1, st0, en0, st_, en_; + int8_t x1, v1; + uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v; + __m128i x1_, v1_; + // find the boundaries + if (st < r - qlen + 1) st = r - qlen + 1; + if (en > r) en = r; + if (st < (r-wr+1)>>1) st = (r-wr+1)>>1; // take the ceil + if (en > (r+wl)>>1) en = (r+wl)>>1; // take the floor + if (st > en) { + ez->zdropped = 1; + break; + } + st0 = st, en0 = en; + st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1; + // set boundary conditions + if (st > 0) { + if (st - 1 >= last_st && st - 1 <= last_en) + x1 = ((uint8_t*)x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round + else x1 = v1 = 0; // not calculated; set to zeros + } else x1 = 0, v1 = r? q : 0; + if (en >= r) ((uint8_t*)y)[r] = 0, u8[r] = r? q : 0; + // loop fission: set scores first + if (!(flag & KSW_EZ_GENERIC_SC)) { + for (t = st0; t <= en0; t += 16) { + __m128i sq, st, tmp, mask; + sq = _mm_loadu_si128((__m128i*)&sf[t]); + st = _mm_loadu_si128((__m128i*)&qrr[t]); + mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_)); + tmp = _mm_cmpeq_epi8(sq, st); +#ifdef __SSE4_1__ + tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp); + tmp = _mm_blendv_epi8(tmp, sc_N_, mask); +#else + tmp = _mm_or_si128(_mm_andnot_si128(tmp, sc_mis_), _mm_and_si128(tmp, sc_mch_)); + tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp), _mm_and_si128(mask, sc_N_)); +#endif + _mm_storeu_si128((__m128i*)((uint8_t*)s + t), tmp); + } + } else { + for (t = st0; t <= en0; ++t) + ((uint8_t*)s)[t] = mat[sf[t] * m + qrr[t]]; + } + // core loop + x1_ = _mm_cvtsi32_si128(x1); + v1_ = _mm_cvtsi32_si128(v1); + st_ = st / 16, en_ = en / 16; + assert(en_ - st_ + 1 <= n_col_); + if (!with_cigar) { // score only + for (t = st_; t <= en_; ++t) { + __m128i z, a, b, xt1, vt1, ut, tmp; + __dp_code_block1; +#ifdef __SSE4_1__ + z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) +#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() + z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; + z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative +#endif + __dp_code_block2; +#ifdef __SSE4_1__ + _mm_store_si128(&x[t], _mm_max_epi8(a, zero_)); + _mm_store_si128(&y[t], _mm_max_epi8(b, zero_)); +#else + tmp = _mm_cmpgt_epi8(a, zero_); + _mm_store_si128(&x[t], _mm_and_si128(a, tmp)); + tmp = _mm_cmpgt_epi8(b, zero_); + _mm_store_si128(&y[t], _mm_and_si128(b, tmp)); +#endif + } + } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment + __m128i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + for (t = st_; t <= en_; ++t) { + __m128i d, z, a, b, xt1, vt1, ut, tmp; + __dp_code_block1; + d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0 +#ifdef __SSE4_1__ + z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) + tmp = _mm_cmpgt_epi8(b, z); + d = _mm_blendv_epi8(d, flag2_, tmp); // d = b > z? 2 : d +#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8() + z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; + z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative + tmp = _mm_cmpgt_epi8(b, z); + d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv +#endif + __dp_code_block2; + tmp = _mm_cmpgt_epi8(a, zero_); + _mm_store_si128(&x[t], _mm_and_si128(tmp, a)); + d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_)); // d = a > 0? 0x08 : 0 + tmp = _mm_cmpgt_epi8(b, zero_); + _mm_store_si128(&y[t], _mm_and_si128(tmp, b)); + d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0 + _mm_store_si128(&pr[t], d); + } + } else { // gap right-alignment + __m128i *pr = p + (size_t)r * n_col_ - st_; + off[r] = st, off_end[r] = en; + for (t = st_; t <= en_; ++t) { + __m128i d, z, a, b, xt1, vt1, ut, tmp; + __dp_code_block1; + d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1 +#ifdef __SSE4_1__ + z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) + tmp = _mm_cmpgt_epi8(z, b); + d = _mm_blendv_epi8(flag2_, d, tmp); // d = z > b? d : 2 +#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8() + z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; + z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative + tmp = _mm_cmpgt_epi8(z, b); + d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv +#endif + __dp_code_block2; + tmp = _mm_cmpgt_epi8(zero_, a); + _mm_store_si128(&x[t], _mm_andnot_si128(tmp, a)); + d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_)); // d = 0 > a? 0 : 0x08 + tmp = _mm_cmpgt_epi8(zero_, b); + _mm_store_si128(&y[t], _mm_andnot_si128(tmp, b)); + d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10 + _mm_store_si128(&pr[t], d); + } + } + if (!approx_max) { // find the exact max with a 32-bit score array + int32_t max_H, max_t; + // compute H[], max_H and max_t + if (r > 0) { + int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i; + __m128i max_H_, max_t_, qe_; + max_H = H[en0] = en0 > 0? H[en0-1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element + max_t = en0; + max_H_ = _mm_set1_epi32(max_H); + max_t_ = _mm_set1_epi32(max_t); + qe_ = _mm_set1_epi32(q + e); + for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; + __m128i H1, tmp, t_; + H1 = _mm_loadu_si128((__m128i*)&H[t]); + t_ = _mm_setr_epi32(v8[t], v8[t+1], v8[t+2], v8[t+3]); + H1 = _mm_add_epi32(H1, t_); + H1 = _mm_sub_epi32(H1, qe_); + _mm_storeu_si128((__m128i*)&H[t], H1); + t_ = _mm_set1_epi32(t); + tmp = _mm_cmpgt_epi32(H1, max_H_); +#ifdef __SSE4_1__ + max_H_ = _mm_blendv_epi8(max_H_, H1, tmp); + max_t_ = _mm_blendv_epi8(max_t_, t_, tmp); +#else + max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_)); + max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_)); +#endif + } + _mm_storeu_si128((__m128i*)HH, max_H_); + _mm_storeu_si128((__m128i*)tt, max_t_); + for (i = 0; i < 4; ++i) + if (max_H < HH[i]) max_H = HH[i], max_t = tt[i] + i; + for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE + H[t] += (int32_t)v8[t] - qe; + if (H[t] > max_H) + max_H = H[t], max_t = t; + } + } else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0 + // update ez + if (en0 == tlen - 1 && H[en0] > ez->mte) + ez->mte = H[en0], ez->mte_q = r - en; + if (r - st0 == qlen - 1 && H[st0] > ez->mqe) + ez->mqe = H[st0], ez->mqe_t = st0; + if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break; + if (r == qlen + tlen - 2 && en0 == tlen - 1) + ez->score = H[tlen - 1]; + } else { // find approximate max; Z-drop might be inaccurate, too. + if (r > 0) { + if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { + int32_t d0 = v8[last_H0_t] - qe; + int32_t d1 = u8[last_H0_t + 1] - qe; + if (d0 > d1) H0 += d0; + else H0 += d1, ++last_H0_t; + } else if (last_H0_t >= st0 && last_H0_t <= en0) { + H0 += v8[last_H0_t] - qe; + } else { + ++last_H0_t, H0 += u8[last_H0_t] - qe; + } + if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) break; + } else H0 = v8[0] - qe - qe, last_H0_t = 0; + if (r == qlen + tlen - 2 && en0 == tlen - 1) + ez->score = H0; + } + last_st = st, last_en = en; + //for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging + } + kfree(km, mem); + if (!approx_max) kfree(km, H); + if (with_cigar) { // backtrack + int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); + if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) { + ez->reach_end = 1; + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } else if (ez->max_t >= 0 && ez->max_q >= 0) { + ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + } + kfree(km, mem2); kfree(km, off); + } +} +#endif // __SSE2__ diff --git a/ext/ksw2_extz2_sse.c b/ext/ksw2_extz2_sse.c deleted file mode 100644 index 767a749f..00000000 --- a/ext/ksw2_extz2_sse.c +++ /dev/null @@ -1,375 +0,0 @@ -#include -#include -#include "ksw2.h" - -#ifdef __SSE2__ -#include - -#ifdef KSW_SSE2_ONLY -#undef __SSE4_1__ -#endif - -#ifdef __SSE4_1__ -#include -#endif - -#ifdef KSW_CPU_DISPATCH -#ifdef __SSE4_1__ -void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) -#else -void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) -#endif -#else -void ksw_extz2_sse(void *km, - int qlen, - const uint8_t *query, - int tlen, - const uint8_t *target, - int8_t m, - const int8_t *mat, - int8_t q, - int8_t e, - int w, - int zdrop, - int end_bonus, - int flag, - ksw_extz_t *ez) -#endif // ~KSW_CPU_DISPATCH -{ -#define __dp_code_block1 \ - z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \ - xt1 = _mm_load_si128(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ - tmp = _mm_srli_si128(xt1, 15); /* tmp <- x[r-1][t+15] */ \ - xt1 = _mm_or_si128(_mm_slli_si128(xt1, 1), x1_); /* xt1 <- x[r-1][t-1..t+14] */ \ - x1_ = tmp; \ - vt1 = _mm_load_si128(&v[t]); /* vt1 <- v[r-1][t..t+15] */ \ - tmp = _mm_srli_si128(vt1, 15); /* tmp <- v[r-1][t+15] */ \ - vt1 = _mm_or_si128(_mm_slli_si128(vt1, 1), v1_); /* vt1 <- v[r-1][t-1..t+14] */ \ - v1_ = tmp; \ - a = _mm_add_epi8(xt1, vt1); /* a <- x[r-1][t-1..t+14] + v[r-1][t-1..t+14] */ \ - ut = _mm_load_si128(&u[t]); /* ut <- u[t..t+15] */ \ - b = _mm_add_epi8(_mm_load_si128(&y[t]), ut); /* b <- y[r-1][t..t+15] + u[r-1][t..t+15] */ - -#define __dp_code_block2 \ - z = _mm_max_epu8(z, b); /* z = max(z, b); this works because both are non-negative */ \ - z = _mm_min_epu8(z, max_sc_); \ - _mm_store_si128(&u[t], _mm_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ - _mm_store_si128(&v[t], _mm_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ - z = _mm_sub_epi8(z, q_); \ - a = _mm_sub_epi8(a, z); \ - b = _mm_sub_epi8(b, z); - - int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc; - int with_cigar = !(flag & KSW_EZ_SCORE_ONLY), approx_max = !!(flag & KSW_EZ_APPROX_MAX); - int32_t *H = 0, H0 = 0, last_H0_t = 0; - uint8_t *qr, *sf, *mem, *mem2 = 0; - __m128i q_, qe2_, zero_, flag1_, flag2_, flag8_, flag16_, sc_mch_, sc_mis_, sc_N_, m1_, max_sc_; - __m128i *u, *v, *x, *y, *s, *p = 0; - - ksw_reset_extz(ez); - if (m <= 0 || qlen <= 0 || tlen <= 0) - return; - - zero_ = _mm_set1_epi8(0); - q_ = _mm_set1_epi8(q); - qe2_ = _mm_set1_epi8((q + e) * 2); - flag1_ = _mm_set1_epi8(1); - flag2_ = _mm_set1_epi8(2); - flag8_ = _mm_set1_epi8(0x08); - flag16_ = _mm_set1_epi8(0x10); - sc_mch_ = _mm_set1_epi8(mat[0]); - sc_mis_ = _mm_set1_epi8(mat[1]); - sc_N_ = mat[m * m - 1] == 0 ? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m * m - 1]); - m1_ = _mm_set1_epi8(m - 1); // wildcard - max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2); - - if (w < 0) - w = tlen > qlen ? tlen : qlen; - wl = wr = w; - tlen_ = (tlen + 15) / 16; - n_col_ = qlen < tlen ? qlen : tlen; - n_col_ = ((n_col_ < w + 1 ? n_col_ : w + 1) + 15) / 16 + 1; - qlen_ = (qlen + 15) / 16; - for (t = 1, max_sc = mat[0], min_sc = mat[1]; t < m * m; ++t) { - max_sc = max_sc > mat[t] ? max_sc : mat[t]; - min_sc = min_sc < mat[t] ? min_sc : mat[t]; - } - if (-min_sc > 2 * (q + e)) - return; // otherwise, we won't see any mismatches - - mem = (uint8_t *) kcalloc(km, tlen_ * 6 + qlen_ + 1, 16); - u = (__m128i *) (((size_t) mem + 15) >> 4 << 4); // 16-byte aligned - v = u + tlen_, x = v + tlen_, y = x + tlen_, s = y + tlen_, sf = (uint8_t *) (s + tlen_), qr = sf + tlen_ * 16; - if (!approx_max) { - H = (int32_t *) kmalloc(km, tlen_ * 16 * 4); - for (t = 0; t < tlen_ * 16; ++t) - H[t] = KSW_NEG_INF; - } - if (with_cigar) { - mem2 = (uint8_t *) kmalloc(km, ((size_t) (qlen + tlen - 1) * n_col_ + 1) * 16); - p = (__m128i *) (((size_t) mem2 + 15) >> 4 << 4); - off = (int *) kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2); - off_end = off + qlen + tlen - 1; - } - - for (t = 0; t < qlen; ++t) - qr[t] = query[qlen - 1 - t]; - memcpy(sf, target, tlen); - - for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) { - int st = 0, en = tlen - 1, st0, en0, st_, en_; - int8_t x1, v1; - uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t *) u, *v8 = (uint8_t *) v; - __m128i x1_, v1_; - // find the boundaries - if (st < r - qlen + 1) - st = r - qlen + 1; - if (en > r) - en = r; - if (st < (r - wr + 1) >> 1) - st = (r - wr + 1) >> 1; // take the ceil - if (en > (r + wl) >> 1) - en = (r + wl) >> 1; // take the floor - if (st > en) { - ez->zdropped = 1; - break; - } - st0 = st, en0 = en; - st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1; - // set boundary conditions - if (st > 0) { - if (st - 1 >= last_st && st - 1 <= last_en) - x1 = ((uint8_t *) x)[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round - else - x1 = v1 = 0; // not calculated; set to zeros - } else - x1 = 0, v1 = r ? q : 0; - if (en >= r) - ((uint8_t *) y)[r] = 0, u8[r] = r ? q : 0; - // loop fission: set scores first - if (!(flag & KSW_EZ_GENERIC_SC)) { - for (t = st0; t <= en0; t += 16) { - __m128i sq, st, tmp, mask; - sq = _mm_loadu_si128((__m128i *) &sf[t]); - st = _mm_loadu_si128((__m128i *) &qrr[t]); - mask = _mm_or_si128(_mm_cmpeq_epi8(sq, m1_), _mm_cmpeq_epi8(st, m1_)); - tmp = _mm_cmpeq_epi8(sq, st); -#ifdef __SSE4_1__ - tmp = _mm_blendv_epi8(sc_mis_, sc_mch_, tmp); - tmp = _mm_blendv_epi8(tmp, sc_N_, mask); -#else - tmp = _mm_or_si128(_mm_andnot_si128(tmp, sc_mis_), _mm_and_si128(tmp, sc_mch_)); - tmp = _mm_or_si128(_mm_andnot_si128(mask, tmp), _mm_and_si128(mask, sc_N_)); -#endif - _mm_storeu_si128((__m128i *) ((uint8_t *) s + t), tmp); - } - } else { - for (t = st0; t <= en0; ++t) - ((uint8_t *) s)[t] = mat[sf[t] * m + qrr[t]]; - } - // core loop - x1_ = _mm_cvtsi32_si128(x1); - v1_ = _mm_cvtsi32_si128(v1); - st_ = st / 16, en_ = en / 16; - assert(en_ - st_ + 1 <= n_col_); - if (!with_cigar) { // score only - for (t = st_; t <= en_; ++t) { - __m128i z, a, b, xt1, vt1, ut, tmp; - __dp_code_block1; -#ifdef __SSE4_1__ - z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) -#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() - z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; - z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative -#endif - __dp_code_block2; -#ifdef __SSE4_1__ - _mm_store_si128(&x[t], _mm_max_epi8(a, zero_)); - _mm_store_si128(&y[t], _mm_max_epi8(b, zero_)); -#else - tmp = _mm_cmpgt_epi8(a, zero_); - _mm_store_si128(&x[t], _mm_and_si128(a, tmp)); - tmp = _mm_cmpgt_epi8(b, zero_); - _mm_store_si128(&y[t], _mm_and_si128(b, tmp)); -#endif - } - } else if (!(flag & KSW_EZ_RIGHT)) { // gap left-alignment - __m128i *pr = p + (size_t) r * n_col_ - st_; - off[r] = st, off_end[r] = en; - for (t = st_; t <= en_; ++t) { - __m128i d, z, a, b, xt1, vt1, ut, tmp; - __dp_code_block1; - d = _mm_and_si128(_mm_cmpgt_epi8(a, z), flag1_); // d = a > z? 1 : 0 -#ifdef __SSE4_1__ - z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) - tmp = _mm_cmpgt_epi8(b, z); - d = _mm_blendv_epi8(d, flag2_, tmp); // d = b > z? 2 : d -#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8() - z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; - z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative - tmp = _mm_cmpgt_epi8(b, z); - d = _mm_or_si128(_mm_andnot_si128(tmp, d), _mm_and_si128(tmp, flag2_)); // d = b > z? 2 : d; emulating blendv -#endif - __dp_code_block2; - tmp = _mm_cmpgt_epi8(a, zero_); - _mm_store_si128(&x[t], _mm_and_si128(tmp, a)); - d = _mm_or_si128(d, _mm_and_si128(tmp, flag8_)); // d = a > 0? 0x08 : 0 - tmp = _mm_cmpgt_epi8(b, zero_); - _mm_store_si128(&y[t], _mm_and_si128(tmp, b)); - d = _mm_or_si128(d, _mm_and_si128(tmp, flag16_)); // d = b > 0? 0x10 : 0 - _mm_store_si128(&pr[t], d); - } - } else { // gap right-alignment - __m128i *pr = p + (size_t) r * n_col_ - st_; - off[r] = st, off_end[r] = en; - for (t = st_; t <= en_; ++t) { - __m128i d, z, a, b, xt1, vt1, ut, tmp; - __dp_code_block1; - d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), flag1_); // d = z > a? 0 : 1 -#ifdef __SSE4_1__ - z = _mm_max_epi8(z, a); // z = z > a? z : a (signed) - tmp = _mm_cmpgt_epi8(z, b); - d = _mm_blendv_epi8(flag2_, d, tmp); // d = z > b? d : 2 -#else // we need to emulate SSE4.1 intrinsics _mm_max_epi8() and _mm_blendv_epi8() - z = _mm_and_si128(z, _mm_cmpgt_epi8(z, zero_)); // z = z > 0? z : 0; - z = _mm_max_epu8(z, a); // z = max(z, a); this works because both are non-negative - tmp = _mm_cmpgt_epi8(z, b); - d = _mm_or_si128(_mm_andnot_si128(tmp, flag2_), _mm_and_si128(tmp, d)); // d = z > b? d : 2; emulating blendv -#endif - __dp_code_block2; - tmp = _mm_cmpgt_epi8(zero_, a); - _mm_store_si128(&x[t], _mm_andnot_si128(tmp, a)); - d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag8_)); // d = 0 > a? 0 : 0x08 - tmp = _mm_cmpgt_epi8(zero_, b); - _mm_store_si128(&y[t], _mm_andnot_si128(tmp, b)); - d = _mm_or_si128(d, _mm_andnot_si128(tmp, flag16_)); // d = 0 > b? 0 : 0x10 - _mm_store_si128(&pr[t], d); - } - } - if (!approx_max) { // find the exact max with a 32-bit score array - int32_t max_H, max_t; - // compute H[], max_H and max_t - if (r > 0) { - int32_t HH[4], tt[4], en1 = st0 + (en0 - st0) / 4 * 4, i; - __m128i max_H_, max_t_, qe_; - max_H = H[en0] = en0 > 0 ? H[en0 - 1] + u8[en0] - qe : H[en0] + v8[en0] - qe; // special casing the last element - max_t = en0; - max_H_ = _mm_set1_epi32(max_H); - max_t_ = _mm_set1_epi32(max_t); - qe_ = _mm_set1_epi32(q + e); - for (t = st0; t < en1; t += 4) { // this implements: H[t]+=v8[t]-qe; if(H[t]>max_H) max_H=H[t],max_t=t; - __m128i H1, tmp, t_; - H1 = _mm_loadu_si128((__m128i *) &H[t]); - t_ = _mm_setr_epi32(v8[t], v8[t + 1], v8[t + 2], v8[t + 3]); - H1 = _mm_add_epi32(H1, t_); - H1 = _mm_sub_epi32(H1, qe_); - _mm_storeu_si128((__m128i *) &H[t], H1); - t_ = _mm_set1_epi32(t); - tmp = _mm_cmpgt_epi32(H1, max_H_); -#ifdef __SSE4_1__ - max_H_ = _mm_blendv_epi8(max_H_, H1, tmp); - max_t_ = _mm_blendv_epi8(max_t_, t_, tmp); -#else - max_H_ = _mm_or_si128(_mm_and_si128(tmp, H1), _mm_andnot_si128(tmp, max_H_)); - max_t_ = _mm_or_si128(_mm_and_si128(tmp, t_), _mm_andnot_si128(tmp, max_t_)); -#endif - } - _mm_storeu_si128((__m128i *) HH, max_H_); - _mm_storeu_si128((__m128i *) tt, max_t_); - for (i = 0; i < 4; ++i) - if (max_H < HH[i]) - max_H = HH[i], max_t = tt[i] + i; - for (; t < en0; ++t) { // for the rest of values that haven't been computed with SSE - H[t] += (int32_t) v8[t] - qe; - if (H[t] > max_H) - max_H = H[t], max_t = t; - } - } else - H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0 - // update ez - if (en0 == tlen - 1 && H[en0] > ez->mte) - ez->mte = H[en0], ez->mte_q = r - en; - if (r - st0 == qlen - 1 && H[st0] > ez->mqe) - ez->mqe = H[st0], ez->mqe_t = st0; - if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) - break; - if (r == qlen + tlen - 2 && en0 == tlen - 1) - ez->score = H[tlen - 1]; - } else { // find approximate max; Z-drop might be inaccurate, too. - if (r > 0) { - if (last_H0_t >= st0 && last_H0_t <= en0 && last_H0_t + 1 >= st0 && last_H0_t + 1 <= en0) { - int32_t d0 = v8[last_H0_t] - qe; - int32_t d1 = u8[last_H0_t + 1] - qe; - if (d0 > d1) - H0 += d0; - else - H0 += d1, ++last_H0_t; - } else if (last_H0_t >= st0 && last_H0_t <= en0) { - H0 += v8[last_H0_t] - qe; - } else { - ++last_H0_t, H0 += u8[last_H0_t] - qe; - } - if ((flag & KSW_EZ_APPROX_DROP) && ksw_apply_zdrop(ez, 1, H0, r, last_H0_t, zdrop, e)) - break; - } else - H0 = v8[0] - qe - qe, last_H0_t = 0; - if (r == qlen + tlen - 2 && en0 == tlen - 1) - ez->score = H0; - } - last_st = st, last_en = en; - //for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging - } - kfree(km, mem); - if (!approx_max) - kfree(km, H); - if (with_cigar) { // backtrack - int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR); - if (!ez->zdropped && !(flag & KSW_EZ_EXTZ_ONLY)) { - ksw_backtrack(km, - 1, - rev_cigar, - 0, - (uint8_t *) p, - off, - off_end, - n_col_ * 16, - tlen - 1, - qlen - 1, - &ez->m_cigar, - &ez->n_cigar, - &ez->cigar); - } else if (!ez->zdropped && (flag & KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int) ez->max) { - ez->reach_end = 1; - ksw_backtrack(km, - 1, - rev_cigar, - 0, - (uint8_t *) p, - off, - off_end, - n_col_ * 16, - ez->mqe_t, - qlen - 1, - &ez->m_cigar, - &ez->n_cigar, - &ez->cigar); - } else if (ez->max_t >= 0 && ez->max_q >= 0) { - ksw_backtrack(km, - 1, - rev_cigar, - 0, - (uint8_t *) p, - off, - off_end, - n_col_ * 16, - ez->max_t, - ez->max_q, - &ez->m_cigar, - &ez->n_cigar, - &ez->cigar); - } - kfree(km, mem2); - kfree(km, off); - } -} -#endif // __SSE2__ \ No newline at end of file diff --git a/src/aligner.cpp b/src/aligner.cpp index 6d775739..93e8396c 100644 --- a/src/aligner.cpp +++ b/src/aligner.cpp @@ -3,14 +3,19 @@ * * This is for anything that returns an aln_info object, currently * Aligner::align and hamming_align. + * + * ksw_extend code is based on https://github.com/lh3/ksw2/blob/master/cli.c */ #include #include #include #include +#include // memset +#include +#include "ksw2/ksw2.h" #include "aligner.hpp" -AlignmentInfo Aligner::align(const std::string &query, const std::string &ref) const { +AlignmentInfo Aligner::align(const std::string& query, const std::string_view ref) const { m_align_calls++; AlignmentInfo aln; int32_t maskLen = query.length() / 2; @@ -25,8 +30,8 @@ AlignmentInfo Aligner::align(const std::string &query, const std::string &ref) c StripedSmithWaterman::Alignment alignment_ssw; - // query must be NULL-terminated - auto flag = ssw_aligner.Align(query.c_str(), ref.c_str(), ref.size(), filter, &alignment_ssw, maskLen); + // only query must be NULL-terminated + auto flag = ssw_aligner.Align(query.c_str(), ref.begin(), ref.size(), filter, &alignment_ssw, maskLen); if (flag != 0) { aln.edit_distance = 100000; aln.ref_start = 0; @@ -116,7 +121,7 @@ AlignmentInfo Aligner::align(const std::string &query, const std::string &ref) c * of the query, once for each end. */ std::tuple highest_scoring_segment( - const std::string& query, const std::string& ref, int match, int mismatch, int end_bonus + const std::string& query, const std::string_view ref, int match, int mismatch, int end_bonus ) { size_t n = query.length(); @@ -151,7 +156,7 @@ std::tuple highest_scoring_segment( } AlignmentInfo hamming_align( - const std::string &query, const std::string &ref, int match, int mismatch, int end_bonus + const std::string &query, const std::string_view ref, int match, int mismatch, int end_bonus ) { AlignmentInfo aln; if (query.length() != ref.length()) { @@ -199,3 +204,133 @@ AlignmentInfo hamming_align( aln.query_end = segment_end; return aln; } + +namespace { + +unsigned char seq_nt4_table[256] = { + 0, 1, 2, 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, 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 +}; + +} // namespace + +void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t wildcard_score) +{ + int i, j; + a = a < 0? -a : a; + b = b > 0? -b : b; + wildcard_score = wildcard_score > 0 ? -wildcard_score : wildcard_score; + for (i = 0; i < m - 1; ++i) { + for (j = 0; j < m - 1; ++j) + mat[i * m + j] = i == j? a : b; + mat[i * m + m - 1] = wildcard_score; + } + for (j = 0; j < m; ++j) + mat[(m - 1) * m + j] = wildcard_score; +} + +std::ostream& operator<<(std::ostream& os, const ksw_extz_t& ez) { + os << "ksw_extz_t(" + // + << "\n max: " << ez.max // max overall score + << "\n coord max_q: " << ez.max_q // max extension coordinate + << "\n coord max_t: " << ez.max_t // max extension coordinate + + << "\n score mqe: " << ez.mqe // max score when reaching the end of query + << "\n mqe_t: " << ez.mqe_t // coordinate in target corresponding to mqe + + << "\n score mte: " << ez.mte // max score when reaching the end of target + << "\n mte_q: " << ez.mte_q // coordinate in query corresponding to mte + + << "\n score both ends: " << ez.score // max score reaching both ends + << "\n cigar: " << Cigar(ez.cigar, ez.n_cigar) + << "\n zdropped: " << ez.zdropped + << "\n reach_end: " << ez.reach_end + << "\n)"; + return os; +} + +AlignmentInfo Aligner::ksw_extend(const std::string& query, const std::string_view ref, bool right_align) const { + AlignmentInfo info; + if (query.size() == 0) { + info.cigar = Cigar(); + info.edit_distance = 0; + info.ref_start = 0; + info.query_start = 0; + info.ref_end = 0; + info.query_end = 0; + info.sw_score = parameters.end_bonus; + return info; + } + + int w = -1; // band width; -1 is inf + int zdrop = -1; // -1 to disable + int flag = KSW_EZ_EXTZ_ONLY | KSW_EZ_GENERIC_SC; + if (right_align) { + flag |= KSW_EZ_RIGHT; + } + ksw_extz_t ez; + memset(&ez, 0, sizeof(ksw_extz_t)); + + ez.max_q = ez.max_t = ez.mqe_t = ez.mte_q = -1; + ez.max = 0; ez.mqe = ez.mte = KSW_NEG_INF; + ez.n_cigar = 0; + int qlen = query.length(); + int tlen = ref.length(); + uint8_t *qseq = (uint8_t*)calloc(qlen + 33, 1); + uint8_t *tseq = (uint8_t*)calloc(tlen + 33, 1); + for (int i = 0; i < qlen; ++i) + qseq[i] = seq_nt4_table[(uint8_t)query[i]]; + for (int i = 0; i < tlen; ++i) + tseq[i] = seq_nt4_table[(uint8_t)ref[i]]; + + ksw_extz2_sse( + nullptr, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, ksw_matrix_m, ksw_matrix, parameters.gap_open, parameters.gap_extend, w, zdrop, parameters.end_bonus, flag, &ez + ); + free(qseq); + free(tseq); + + auto cigar = Cigar(ez.cigar, ez.n_cigar).to_eqx(query, ref); + info.edit_distance = cigar.edit_distance(); + info.cigar = std::move(cigar); + info.ref_start = 0; + info.query_start = 0; + if (ez.reach_end) { + info.ref_end = ez.mqe_t + 1; + info.query_end = query.size(); + info.sw_score = ez.mqe + parameters.end_bonus; + } else { + info.ref_end = ez.max_t + 1; + info.query_end = ez.max_q + 1; + info.sw_score = ez.max; + info.cigar.push(CIGAR_SOFTCLIP, query.length() - info.query_end); + } + assert(info.cigar.derived_sequence_length() == query.length()); + + kfree(km, ez.cigar); + return info; +} + +std::ostream& operator<<(std::ostream& os, const AlignmentInfo& info) { + os << "AlignmentInfo(cigar='" << info.cigar + << "', ref=" << info.ref_start << ".." << info.ref_end + << ", query=" << info.query_start << ".." << info.query_end + << ", NM=" << info.edit_distance + << ", AS=" << info.sw_score + << ")"; + return os; +} diff --git a/src/aligner.hpp b/src/aligner.hpp index 0238885f..77f83c23 100644 --- a/src/aligner.hpp +++ b/src/aligner.hpp @@ -26,16 +26,25 @@ struct AlignmentInfo { int sw_score{0}; int ref_span() const { return ref_end - ref_start; } + int query_span() const { return query_end - query_start; } }; +std::ostream& operator<<(std::ostream& os, const AlignmentInfo& info); + +void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t wildcard_score); + struct Aligner { public: Aligner(AlignmentParameters parameters) : parameters(parameters) , ssw_aligner(StripedSmithWaterman::Aligner(parameters.match, parameters.mismatch, parameters.gap_open, parameters.gap_extend)) - { } + { + int8_t wildcard_score = -parameters.mismatch; + ksw_gen_simple_mat(ksw_matrix_m, ksw_matrix, parameters.match, -parameters.mismatch, wildcard_score); + } - AlignmentInfo align(const std::string &query, const std::string &ref) const; + AlignmentInfo align(const std::string& query, const std::string_view ref) const; + AlignmentInfo ksw_extend(const std::string& query, const std::string_view ref, bool right_align) const; AlignmentParameters parameters; @@ -47,9 +56,11 @@ struct Aligner { const StripedSmithWaterman::Aligner ssw_aligner; const StripedSmithWaterman::Filter filter; mutable unsigned m_align_calls{0}; // no. of calls to the align() method + const int8_t ksw_matrix_m{5}; + int8_t ksw_matrix[25]; }; -inline int hamming_distance(const std::string &s, const std::string &t) { +inline int hamming_distance(const std::string& s, const std::string_view t) { if (s.length() != t.length()){ return -1; } @@ -65,11 +76,11 @@ inline int hamming_distance(const std::string &s, const std::string &t) { } std::tuple highest_scoring_segment( - const std::string& query, const std::string& ref, int match, int mismatch, int end_bonus + const std::string& query, const std::string_view ref, int match, int mismatch, int end_bonus ); AlignmentInfo hamming_align( - const std::string &query, const std::string &ref, int match, int mismatch, int end_bonus + const std::string& query, const std::string_view ref, int match, int mismatch, int end_bonus ); #endif diff --git a/src/aln.cpp b/src/aln.cpp index b5105bb1..a8bae9c3 100644 --- a/src/aln.cpp +++ b/src/aln.cpp @@ -213,7 +213,7 @@ inline Alignment extend_seed( bool consistent_nam ) { const std::string query = nam.is_rc ? read.rc : read.seq; - const std::string& ref = references.sequences[nam.ref_id]; + const std::string_view ref = references.sequences[nam.ref_id]; const auto projected_ref_start = nam.projected_ref_start(); const auto projected_ref_end = std::min(nam.ref_end + query.size() - nam.query_end, ref.size()); @@ -222,11 +222,48 @@ inline Alignment extend_seed( int result_ref_start; bool gapped = true; if (projected_ref_end - projected_ref_start == query.size() && consistent_nam) { - std::string ref_segm_ham = ref.substr(projected_ref_start, query.size()); - auto hamming_dist = hamming_distance(query, ref_segm_ham); + const std::string_view projected_ref = ref.substr(projected_ref_start, query.size()); + info = hamming_align(query, projected_ref, aligner.parameters.match, aligner.parameters.mismatch, aligner.parameters.end_bonus); + + if (info.edit_distance + (query.length() - info.query_span()) < 0.05 * query.length()) { + if (info.query_end < query.length()) { + // Right end is soft clipped, do gapped alignment on it + const std::string query_right = query.substr(info.query_end); + const int ext_right = std::min(ref.size() - projected_ref_end, size_t(50)); + const std::string_view ref_right = ref.substr(projected_ref_start + info.ref_end, ext_right); + auto right = aligner.ksw_extend(query_right, ref_right, false); + info.query_end += right.query_end; + info.ref_end += right.ref_end; + info.edit_distance += right.edit_distance; + info.sw_score += right.sw_score; + assert(!info.cigar.empty()); + info.cigar.pop_oplen(); + info.cigar += right.cigar; + } + + if (info.query_start > 0) { + // Left end is soft clipped, do gapped alignment on it + std::string query_left = query.substr(0, info.query_start); + const int ext_left = std::min(50, projected_ref_start); + const int ref_start = projected_ref_start - ext_left; + std::string ref_left{ref.substr(ref_start, ext_left + info.ref_start)}; + std::reverse(query_left.begin(), query_left.end()); + std::reverse(ref_left.begin(), ref_left.end()); + auto left = aligner.ksw_extend(query_left, ref_left, true); + info.query_start -= left.query_end; + info.ref_start -= left.ref_end; + info.edit_distance += left.edit_distance; + info.sw_score += left.sw_score; + + // TODO this just removes the soft-clipping from the beginning, + // a bit too complicated + info.cigar.reverse(); + info.cigar.pop_oplen(); + info.cigar += left.cigar; + info.cigar.reverse(); + } - if (hamming_dist >= 0 && (((float) hamming_dist / query.size()) < 0.05) ) { //Hamming distance worked fine, no need to ksw align - info = hamming_align(query, ref_segm_ham, aligner.parameters.match, aligner.parameters.mismatch, aligner.parameters.end_bonus); + assert(info.cigar.derived_sequence_length() == query.length()); result_ref_start = projected_ref_start + info.ref_start; gapped = false; } diff --git a/src/cigar.cpp b/src/cigar.cpp index 88676d53..c2799a1c 100644 --- a/src/cigar.cpp +++ b/src/cigar.cpp @@ -17,7 +17,7 @@ Cigar Cigar::to_m() const { return cigar; } -Cigar Cigar::to_eqx(const std::string& query, const std::string& ref) const { +Cigar Cigar::to_eqx(const std::string& query, const std::string_view ref) const { size_t i = 0, j = 0; Cigar cigar; for (auto op_len : m_ops) { @@ -89,6 +89,11 @@ Cigar::Cigar(const std::string& cig) { } } +std::ostream& operator<<(std::ostream& os, const Cigar& cigar) { + os << cigar.to_string(); + return os; +} + std::string compress_cigar(const std::string& ops) { char prev = 0; int count = 0; diff --git a/src/cigar.hpp b/src/cigar.hpp index ced193b1..55bda51d 100644 --- a/src/cigar.hpp +++ b/src/cigar.hpp @@ -77,21 +77,39 @@ class Cigar { return dist; } + size_t derived_sequence_length() const { + size_t length = 0; + for (auto op_len : m_ops) { + auto op = op_len & 0xf; + auto len = op_len >> 4; + if (op == CIGAR_MATCH || op == CIGAR_EQ || op == CIGAR_X || op == CIGAR_INS || op == CIGAR_SOFTCLIP) { + length += len; + } + } + return length; + } + void reverse() { std::reverse(m_ops.begin(), m_ops.end()); } + void pop_oplen() { + m_ops.pop_back(); + } + /* Return a new Cigar that uses M instead of =/X */ Cigar to_m() const; /* Return a new Cigar that uses =/X instead of M */ - Cigar to_eqx(const std::string& query, const std::string& ref) const; + Cigar to_eqx(const std::string& query, const std::string_view ref) const; std::string to_string() const; std::vector m_ops; }; +std::ostream& operator<<(std::ostream& os, const Cigar& cigar); + std::string compress_cigar(const std::string& ops); #endif diff --git a/tests/baseline-commit.txt b/tests/baseline-commit.txt index 471fbbbb..12f828e5 100644 --- a/tests/baseline-commit.txt +++ b/tests/baseline-commit.txt @@ -1 +1 @@ -2e4ff9500e68d6e465735dd276d362cf71851dcd +26b472e7a748aabab5e46328622b1cd46a0c0841 diff --git a/tests/phix.1.fastq b/tests/phix.1.fastq index 148f0248..1ed483a6 100644 --- a/tests/phix.1.fastq +++ b/tests/phix.1.fastq @@ -31,9 +31,9 @@ NTCATTTGGCGAGAAAGCTCAGTCTCAGGAGGAAGCGGAGCAGGCCAAATGTTTTTGAGATGGCAGCAACGGAAACCATA + #88ABCFG=>+@+3,3CFBFFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? +#8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDF>>9>9? @SRR1377138.14 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT + diff --git a/tests/phix.pe.sam b/tests/phix.pe.sam index ef4157d7..ccf0dbca 100644 --- a/tests/phix.pe.sam +++ b/tests/phix.pe.sam @@ -17,7 +17,7 @@ SRR1377138.7 83 NC_001422.1 256 60 300=1X = 141 -416 GGCACGTTCGTCAAGGACTGGTTTAGA SRR1377138.7 163 NC_001422.1 141 60 1X299=1X = 256 416 NTTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTTCAACCACTAATAGGTAAGAAATCATGAGTCAAGTTACTGAACAATCCGTACGTTTCCAGACCGCTTTGGCCTCTT #8AF6BFFCFFF4AA>A4A4>:E>7?-55>;(>?E38*878CB6@EC7C;6C>C<)6FFF@E>2)(4179:BB?FB#### NM:i:2 AS:i:602 RG:Z:1 SRR1377138.8 83 NC_001422.1 1254 60 124=1X132=1X42=1X = 1179 -376 TCATGAAGGATGGTGTTAATGCCACTCCTCTCCCGACTGTTAACACTACTGGTTATATTGACCATGCCGCTTTTCTTGGCACGATTAACCCTGATACCAATAAAATCCCTAAGCATTTGTTTCATGGTTATTTGAATATCTATAACAACTATTTTAAAGCGCCGTGGATGCCTGACCGTACCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGCCTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGAN #################################?<9<464:AC;CFAF;AFFA77AFAGFGFACB73DC4CD64GGGGCGEC=>+FGGFEFGE8GGGGCE=EC7GEFEFF5>A@FFE=<)1*-..:@29@29@0:1;3=6@=(544?4)9;C>5))7).474CB)99?###### NM:i:5 AS:i:572 RG:Z:1 -SRR1377138.9 83 NC_001422.1 4731 60 300=1X = 4650 -382 TGTGTGCCTGAGTATGGTACAGCTAATGGCCGTCTTCATTTCCATGCGGTGCACTTTATGCGGACACTTCCTACAGGTAGCGTTGACCCTAATTTTGGTCGTCGGGTACGCAATCGCCGCCAGTTAAATAGCTTGCAAAATACGTGGCCTTATGGTTACAGTATGCCCATCGCAGTTCGCTACACGCAGGACGCTTTTTCACGTTCTGGTTGGTTGTGGCCTGTTGATGCTAAAGGTGAGCCGCTTAAAGCTACCAGTTATATGGCTGTTGGTTTCTATGTGGCTAAATACGTTAACAAAN ?BFF?<:?FCF@FFF?EFFF@EFEFFFFEEEB?FD<<)8083:?>FFDFFECC49FFD?3<@FFFFEDEC?3E?C?0EFFFFFFF>@5FECFFGGGFGD:FGGCDGFGGGGGG>GGGGGGGDGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:1 AS:i:612 RG:Z:1 @@ -25,7 +25,7 @@ SRR1377138.11 99 NC_001422.1 5117 60 1X268=32S = 5220 270 NAAGCTGTTCAGAATCAGAATG SRR1377138.11 147 NC_001422.1 5220 60 167=134S = 5117 -270 GACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAAN ########BFFBA>82><9B>4<2BAA2A<<:?:4,.(B?C?.37;)EE<;)C?6(((;FFFFFEE?7)?9/);FDFB>A=FDAFFFFFFFDEDDCCFAGFFGGGGGGF:GGGFGGGGGFGGGGFGCGGEECDEFCCGGFGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCA8# NM:i:0 AS:i:344 RG:Z:1 SRR1377138.12 99 NC_001422.1 1 60 57S238=1X5= = 28 328 NAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGATGTGGC #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGDCGGGGGGGGGGFGGGGGGGFGFGGFGFGFFFFFFFFFFFAFFFFFFFFFFFA3.1.48?A<:?FFFFFFBDFFFF?1GGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:2 AS:i:602 RG:Z:1 -SRR1377138.13 99 NC_001422.1 5006 60 1X300= = 5049 344 NTATGTGGCTAAATACGTTAACAAAAAGTCAGATATGGACCTTGCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTT #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 RG:Z:1 +SRR1377138.13 99 NC_001422.1 5006 60 1X293=1D6= = 5049 344 NTATGTGGCTAAATACGTTAACAAAAAGTCAGATATGGACCTTGCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCGACGTT #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDF>>9>9? NM:i:2 AS:i:597 RG:Z:1 SRR1377138.13 147 NC_001422.1 5049 60 300=1X = 5006 -344 GCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGN #FFFFBFFFBD>A<::64>AAFFFB><?FFBABBAFFFFFFBAFFFFFFFFFFFFFE;@F@FFFFFFFD@FE?FFFFFFFFFFFFFEBFFFFFFFFGFFGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCC8# NM:i:1 AS:i:612 RG:Z:1 SRR1377138.14 99 NC_001422.1 66 60 1X300= = 156 391 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT #8@CCGGGGGGFFGGGGGGGGGGGGGFGGGGGGGGGGGGGFGGEEFGGGGGGGGGGGGGFGDCFEFGGGGGFGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGGGGGGGGGDFGGGGGGGFGGEGGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGDGGGGGGGGDGGGGGGGFEGGGGGGGFGGFFGFGGGGGGGGGFGGGFG7FGGFGDGGGGGFFCGGGGGGGGGGGGGFGFGGC>FGCF=A@FFFFFFFFFFEFFFFFFE<=/04(34:501:AB>4 NM:i:1 AS:i:612 RG:Z:1 SRR1377138.14 147 NC_001422.1 156 60 300=1X = 66 -391 GCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTTCAACCACTAATAGGTAAGAAATCATGAGTCAAGTTACTGAACAATCCGTACGTTTCCAGACCGCTTTGGCCTCTATTAAGCTCATTCAGN ############?>BFAAA45FFFFB:7/(9>9;B@=63?FFED>52FE?>5*>852BCEEFFAA@A>8@8C>=7EFCFFEAFFFF:FAFFFFFFFFD?9*DDGGGGGDGGGGGGFGFGGFFCFFGDGFGGGGGFGFEF?GGFCGGGGGGGGGGGGFEEGFFGGFE9GGGFFGFAGGGGGGGGFAGGGFGFFFGGGGGGGEFAFFCFGGGF@GF@AE@FFCEFEEFFFFFFFFFFFFFFF0.044:6B0>BF>003>>:(:<>B?FFFA>:5)9FFFFDAAB<4<60.,AB>FFFFDFFEFA6FFFFFFA8=FGFCFCGFGFF7CGGGGFGGGGFECCGGGGGGGGGGGGGGGGGGGEGGGGGGGFFFGGGGGFEFGGCGGGGGGGGFGEAFGGGFGGGGGGGGGGFCGFFGEGGGGGGGGFGGGGGGGGGDFCFFF<9GGGGGGGGGGGGGGGFFF=8GFGGGGGGGGGGGGGGGGGGEGGGEFFEFGGGGGGGGGGGEFGFFEFGGEFDFCFGGFCCGGDFGGGGGFCGGGGGEGEFGEDEGGGEGGGGGGBCA8# NM:i:1 AS:i:612 RG:Z:1 SRR1377138.8 16 NC_001422.1 1254 60 124=1X132=1X42=1X * 0 0 TCATGAAGGATGGTGTTAATGCCACTCCTCTCCCGACTGTTAACACTACTGGTTATATTGACCATGCCGCTTTTCTTGGCACGATTAACCCTGATACCAATAAAATCCCTAAGCATTTGTTTCATGGTTATTTGAATATCTATAACAACTATTTTAAAGCGCCGTGGATGCCTGACCGTACCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGCCTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGAN #################################?<9<464:AC;CFAF;AFFA77AFAGFGFACB73DC4CD64GGGGCGEC=>+FGGFEFGE8GGGGCE=EC7GEF*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 RG:Z:1 +SRR1377138.13 0 NC_001422.1 5006 60 1X293=1D6= * 0 0 NTATGTGGCTAAATACGTTAACAAAAAGTCAGATATGGACCTTGCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCGACGTT #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDF>>9>9? NM:i:2 AS:i:597 RG:Z:1 SRR1377138.14 0 NC_001422.1 66 60 1X300= * 0 0 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT #8@CCGGGGGGFFGGGGGGGGGGGGGFGGGGGGGGGGGGGFGGEEFGGGGGGGGGGGGGFGDCFEFGGGGGFGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGGGGGGGGGDFGGGGGGGFGGEGGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGDGGGGGGGGDGGGGGGGFEGGGGGGGFGGFFGFGGGGGGGGGFGGGFG7FGGFGDGGGGGFFCGGGGGGGGGGGGGFGFGGC>FGCF=A@FFFFFFFFFFEFFFFFFE<=/04(34:501:AB>4 NM:i:1 AS:i:612 RG:Z:1 SRR1377138.15 0 NC_001422.1 1907 60 1X185=115S * 0 0 NTTAACGCTAAAGGTGCTTTGACTTATACCGATATTGCTGGCGACCCTGTTTTGTATGGCAACTTGCCGCCGCGTGAAATTTCTATGAAGGATGTTTTCCGTTCTGGTGATTCGTCTAAGAAGTTTAAGATTGCTGAGGGTCAGTGGTATCGTTATGCGCCTTCGTATGTTTCTCCTGCTTATCACAAGAACGGAAAACATCCTTCATAGAAATTTCACGCAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAAAGATA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFGFFFFFFFFFFFFF>73?FFFFFFF3-048:BDFFAFFFFFF06>BFF######### NM:i:1 AS:i:372 RG:Z:1 SRR1377138.16 16 NC_001422.1 2072 60 27=1X20=1X83=1X50=1X116=1X * 0 0 TATGTTTCTCCTGCTTATCACCTTCTTTAAGGCTTCCCATTCATTCAGCAACCGCCTTCTGGTGATTTGCAAGAACGCGTACTTATTCGCCACCATGATTATGACCAGTGTTTCCAGTCCGTTCAGTTGTTGAAGTGGAATAGTCAGGTTAAATTTAATGTGACCGTTTATCGCAATCTGCCGTCCACTCGCGATTCAATCATGACTTCGTGATAAAAGATTGAGTGTGAGGTTATAACGCCGAAGCGGTAAAAATTTTAATTTTTGCCGCTGAGGGGTTGACCAAGCGAAGCGCGGTAGN ############################################################################################################GC;D>,5E>B*,,;98F>FE,9E4,@F@7,;CF>3,,BF7=@,3CF>,E,,F@83++3>@3++CBC,5B+68+FB,7=>=F=F,