1 2 /* $Id: smith_waterman_sse2.h 625 2011-03-23 17:21:38Z wrp $ */ 3 /* $Revision: 625 $ */ 4 5 /* The MIT License 6 Copyright (c) 2012-1015 Boston College. 7 Permission is hereby granted, free of charge, to any person obtaining 8 a copy of this software and associated documentation files (the 9 "Software"), to deal in the Software without restriction, including 10 without limitation the rights to use, copy, modify, merge, publish, 11 distribute, sublicense, and/or sell copies of the Software, and to 12 permit persons to whom the Software is furnished to do so, subject to 13 the following conditions: 14 The above copyright notice and this permission notice shall be 15 included in all copies or substantial portions of the Software. 16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS 20 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN 21 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 22 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 23 SOFTWARE. 24 */ 25 26 /* 27 Written by Michael Farrar, 2006. 28 Please send bug reports and/or suggestions to farrar.michael@gmail.com. 29 */ 30 31 #ifndef SMITH_WATERMAN_SSE2_H 32 #define SMITH_WATERMAN_SSE2_H 33 34 #include <climits> 35 36 #include <cstdio> 37 #include <cstdlib> 38 39 #if !defined(__APPLE__) && !defined(__llvm__) && !defined(__DragonFly__) 40 #include <malloc.h> 41 #endif 42 43 #include "simd.h" 44 #include "BaseMatrix.h" 45 46 #include "Sequence.h" 47 #include "EvalueComputation.h" 48 typedef struct { 49 short qStartPos; 50 short dbStartPos; 51 short qEndPos; 52 short dbEndPos; 53 } aln_t; 54 55 56 typedef struct { 57 uint32_t score1; 58 uint32_t score2; 59 int32_t dbStartPos1; 60 int32_t dbEndPos1; 61 int32_t qStartPos1; 62 int32_t qEndPos1; 63 int32_t ref_end2; 64 float qCov; 65 float tCov; 66 uint32_t* cigar; 67 int32_t cigarLen; 68 double evalue; 69 } s_align; 70 71 class SmithWaterman{ 72 public: 73 74 SmithWaterman(size_t maxSequenceLength, int aaSize, bool aaBiasCorrection); 75 ~SmithWaterman(); 76 77 // prints a __m128 vector containing 8 signed shorts 78 static void printVector (__m128i v); 79 80 // prints a __m128 vector containing 8 unsigned shorts, added 32768 81 static void printVectorUS (__m128i v); 82 83 static unsigned short sse2_extract_epi16(__m128i v, int pos); 84 85 // The dynamic programming matrix entries for the query and database sequences are stored sequentially (the order see the Farrar paper). 86 // This function calculates the index within the dynamic programming matrices for the given query and database sequence position. midx(int qpos,int dbpos,int iter)87 static inline int midx (int qpos, int dbpos, int iter){ 88 return dbpos * (8 * iter) + (qpos % iter) * 8 + (qpos / iter); 89 } 90 91 // @function ssw alignment. 92 /*! @function Do Striped Smith-Waterman alignment. 93 94 @param db_sequence pointer to the target sequence; the target sequence needs to be numbers and corresponding to the mat parameter of 95 function ssw_init 96 97 @param db_length length of the target sequence 98 99 @param gap_open the absolute value of gap open penalty 100 101 @param gap_extend the absolute value of gap extension penalty 102 103 @param alignmentMode bitwise FLAG; (from high to low) bit 5: when setted as 1, function ssw_align will return the best alignment 104 beginning position; bit 6: when setted as 1, if (ref_end1 - ref_begin1 < filterd && read_end1 - read_begin1 105 < filterd), (whatever bit 5 is setted) the function will return the best alignment beginning position and 106 cigar; bit 7: when setted as 1, if the best alignment score >= filters, (whatever bit 5 is setted) the function 107 will return the best alignment beginning position and cigar; bit 8: when setted as 1, (whatever bit 5, 6 or 7 is 108 setted) the function will always return the best alignment beginning position and cigar. When flag == 0, only 109 the optimal and sub-optimal scores and the optimal alignment ending position will be returned. 110 111 @param filters score filter: when bit 7 of flag is setted as 1 and bit 8 is setted as 0, filters will be used (Please check the 112 decription of the flag parameter for detailed usage.) 113 114 @param filterd distance filter: when bit 6 of flag is setted as 1 and bit 8 is setted as 0, filterd will be used (Please check 115 the decription of the flag parameter for detailed usage.) 116 117 @param maskLen The distance between the optimal and suboptimal alignment ending position >= maskLen. We suggest to use 118 readLen/2, if you don't have special concerns. Note: maskLen has to be >= 15, otherwise this function will NOT 119 return the suboptimal alignment information. Detailed description of maskLen: After locating the optimal 120 alignment ending position, the suboptimal alignment score can be heuristically found by checking the second 121 largest score in the array that contains the maximal score of each column of the SW matrix. In order to avoid 122 picking the scores that belong to the alignments sharing the partial best alignment, SSW C library masks the 123 reference loci nearby (mask length = maskLen) the best alignment ending position and locates the second largest 124 score from the unmasked elements. 125 126 @return pointer to the alignment result structure 127 128 @note Whatever the parameter flag is setted, this function will at least return the optimal and sub-optimal alignment score, 129 and the optimal alignment ending positions on target and query sequences. If both bit 6 and 7 of the flag are setted 130 while bit 8 is not, the function will return cigar only when both criteria are fulfilled. All returned positions are 131 0-based coordinate. 132 */ 133 s_align ssw_align (const unsigned char*db_sequence, 134 int32_t db_length, 135 const uint8_t gap_open, 136 const uint8_t gap_extend, 137 const uint8_t alignmentMode, // (from high to low) bit 5: return the best alignment beginning position; 6: if (ref_end1 - ref_begin1 <= filterd) && (read_end1 - read_begin1 <= filterd), return cigar; 7: if max score >= filters, return cigar; 8: always return cigar; if 6 & 7 are both setted, only return cigar when both filter fulfilled 138 const double filters, 139 EvalueComputation * filterd, 140 const int covMode, const float covThr, 141 const int32_t maskLen); 142 143 144 /*! @function computed ungapped alignment score 145 146 @param db_sequence pointer to the target sequence; the target sequence needs to be numbers and corresponding to the mat parameter of 147 function ssw_init 148 149 @param db_length length of the target sequence 150 @return max diagonal score 151 */ 152 int ungapped_alignment(const unsigned char *db_sequence, 153 int32_t db_length); 154 155 /*! @function Create the query profile using the query sequence. 156 @param read pointer to the query sequence; the query sequence needs to be numbers 157 @param readLen length of the query sequence 158 @param mat pointer to the substitution matrix; mat needs to be corresponding to the read sequence 159 @param n the square root of the number of elements in mat (mat has n*n elements) 160 @param score_size estimated Smith-Waterman score; if your estimated best alignment score is surely < 255 please set 0; if 161 your estimated best alignment score >= 255, please set 1; if you don't know, please set 2 162 @return pointer to the query profile structure 163 @note example for parameter read and mat: 164 If the query sequence is: ACGTATC, the sequence that read points to can be: 1234142 165 Then if the penalty for match is 2 and for mismatch is -2, the substitution matrix of parameter mat will be: 166 //A C G T 167 2 -2 -2 -2 //A 168 -2 2 -2 -2 //C 169 -2 -2 2 -2 //G 170 -2 -2 -2 2 //T 171 mat is the pointer to the array {2, -2, -2, -2, -2, 2, -2, -2, -2, -2, 2, -2, -2, -2, -2, 2} 172 */ 173 void ssw_init(const Sequence *q, const int8_t *mat, const BaseMatrix *m, const int8_t score_size); 174 175 176 static char cigar_int_to_op (uint32_t cigar_int); 177 178 static uint32_t cigar_int_to_len (uint32_t cigar_int); 179 180 181 static float computeCov(unsigned int startPos, unsigned int endPos, unsigned int len); 182 183 s_align scoreIdentical(unsigned char *dbSeq, int L, EvalueComputation * evaluer, int alignmentMode); 184 seq_reverse(int8_t * reverse,const int8_t * seq,int32_t end)185 static void seq_reverse(int8_t * reverse, const int8_t* seq, int32_t end) /* end is 0-based alignment ending position */ 186 { 187 int32_t start = 0; 188 while (LIKELY(start <= end)) { 189 reverse[start] = seq[end]; 190 reverse[end] = seq[start]; 191 ++start; 192 --end; 193 } 194 } 195 196 197 private: 198 199 struct s_profile{ 200 simd_int* profile_byte; // 0: none 201 simd_int* profile_word; // 0: none 202 simd_int* profile_rev_byte; // 0: none 203 simd_int* profile_rev_word; // 0: none 204 int8_t* query_sequence; 205 int8_t* query_rev_sequence; 206 int8_t* composition_bias; 207 int8_t* composition_bias_rev; 208 int8_t* mat; 209 // Memory layout of if mat + queryProfile is qL * AA 210 // Query length 211 // A -1 -3 -2 -1 -4 -2 -2 -3 -1 -3 -2 -2 7 -1 -2 -1 -1 -2 -5 -3 212 // C -1 -4 2 5 -3 -2 0 -3 1 -3 -2 0 -1 2 0 0 -1 -3 -4 -2 213 // ... 214 // Y -1 -3 -2 -1 -4 -2 -2 -3 -1 -3 -2 -2 7 -1 -2 -1 -1 -2 -5 -3 215 // Memory layout of if mat + sub is AA * AA 216 // A C ... Y 217 // A -1 -3 -2 -1 -4 -2 -2 -3 -1 -3 -2 -2 7 -1 -2 -1 -1 -2 -5 -3 218 // C -1 -4 2 5 -3 -2 0 -3 1 -3 -2 0 -1 2 0 0 -1 -3 -4 -2 219 // ... 220 // Y -1 -3 -2 -1 -4 -2 -2 -3 -1 -3 -2 -2 7 -1 -2 -1 -1 -2 -5 -3 221 int8_t* mat_rev; // needed for queryProfile 222 int32_t query_length; 223 int32_t sequence_type; 224 int32_t alphabetSize; 225 uint8_t bias; 226 short ** profile_word_linear; 227 }; 228 simd_int* vHStore; 229 simd_int* vHLoad; 230 simd_int* vE; 231 simd_int* vHmax; 232 uint8_t * maxColumn; 233 234 typedef struct { 235 uint16_t score; 236 int32_t ref; //0-based position 237 int32_t read; //alignment ending position on read, 0-based 238 } alignment_end; 239 240 241 typedef struct { 242 uint32_t* seq; 243 int32_t length; 244 } cigar; 245 246 /* Striped Smith-Waterman 247 Record the highest score of each reference position. 248 Return the alignment score and ending position of the best alignment, 2nd best alignment, etc. 249 Gap begin and gap extension are different. 250 wight_match > 0, all other weights < 0. 251 The returned positions are 0-based. 252 */ 253 std::pair<alignment_end, alignment_end> sw_sse2_byte (const unsigned char*db_sequence, 254 int8_t ref_dir, // 0: forward ref; 1: reverse ref 255 int32_t db_length, 256 int32_t query_length, 257 const uint8_t gap_open, /* will be used as - */ 258 const uint8_t gap_extend, /* will be used as - */ 259 const simd_int* query_profile_byte, 260 uint8_t terminate, /* the best alignment score: used to terminate 261 the matrix calculation when locating the 262 alignment beginning point. If this score 263 is set to 0, it will not be used */ 264 uint8_t bias, /* Shift 0 point to a positive value. */ 265 int32_t maskLen); 266 267 std::pair<alignment_end, alignment_end> sw_sse2_word (const unsigned char* db_sequence, 268 int8_t ref_dir, // 0: forward ref; 1: reverse ref 269 int32_t db_length, 270 int32_t query_length, 271 const uint8_t gap_open, /* will be used as - */ 272 const uint8_t gap_extend, /* will be used as - */ 273 const simd_int*query_profile_byte, 274 uint16_t terminate, 275 int32_t maskLen); 276 277 template <const unsigned int type> 278 SmithWaterman::cigar *banded_sw(const unsigned char *db_sequence, const int8_t *query_sequence, const int8_t * compositionBias, int32_t db_length, int32_t query_length, int32_t queryStart, int32_t score, const uint32_t gap_open, const uint32_t gap_extend, int32_t band_width, const int8_t *mat, int32_t n); 279 280 /*! @function Produce CIGAR 32-bit unsigned integer from CIGAR operation and CIGAR length 281 @param length length of CIGAR 282 @param op_letter CIGAR operation character ('M', 'I', etc) 283 @return 32-bit unsigned integer, representing encoded CIGAR operation and length 284 */ 285 inline uint32_t to_cigar_int (uint32_t length, char op_letter); 286 287 s_profile* profile; 288 289 290 const static unsigned int SUBSTITUTIONMATRIX = 1; 291 const static unsigned int PROFILE = 2; 292 293 template <typename T, size_t Elements, const unsigned int type> 294 void createQueryProfile(simd_int *profile, const int8_t *query_sequence, const int8_t * composition_bias, const int8_t *mat, const int32_t query_length, const int32_t aaSize, uint8_t bias, const int32_t offset, const int32_t entryLength); 295 296 float *tmp_composition_bias; 297 short * profile_word_linear_data; 298 bool aaBiasCorrection; 299 }; 300 #endif /* SMITH_WATERMAN_SSE2_H */ 301