-
Notifications
You must be signed in to change notification settings - Fork 19
/
aligner_swsse_ee_u8.cpp
1901 lines (1764 loc) · 61.2 KB
/
aligner_swsse_ee_u8.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* Copyright 2011, Ben Langmead <[email protected]>
*
* This file is part of Bowtie 2.
*
* Bowtie 2 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Bowtie 2 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Bowtie 2. If not, see <http://www.gnu.org/licenses/>.
*/
/**
* aligner_sw_sse.cpp
*
* Versions of key alignment functions that use vector instructions to
* accelerate dynamic programming. Based chiefly on the striped Smith-Waterman
* paper and implementation by Michael Farrar. See:
*
* Farrar M. Striped Smith-Waterman speeds database searches six times over
* other SIMD implementations. Bioinformatics. 2007 Jan 15;23(2):156-61.
* http://sites.google.com/site/farrarmichael/smith-waterman
*
* While the paper describes an implementation of Smith-Waterman, we extend it
* do end-to-end read alignment as well as local alignment. The change
* required for this is minor: we simply let vmax be the maximum element in the
* score domain rather than the minimum.
*
* The vectorized dynamic programming implementation lacks some features that
* make it hard to adapt to solving the entire dynamic-programming alignment
* problem. For instance:
*
* - It doesn't respect gap barriers on either end of the read
* - It just gives a maximum; not enough information to backtrace without
* redoing some alignment
* - It's a little difficult to handle st_ and en_, especially st_.
* - The query profile mechanism makes handling of ambiguous reference bases a
* little tricky (16 cols in query profile lookup table instead of 5)
*
* Given the drawbacks, it is tempting to use SSE dynamic programming as a
* filter rather than as an aligner per se. Here are a few ideas for how it
* can be extended to handle more of the alignment problem:
*
* - Save calculated scores to a big array as we go. We return to this array
* to find and backtrace from good solutions.
*/
#include <limits>
#include "aligner_sw.h"
static const size_t NBYTES_PER_REG = 16;
static const size_t NWORDS_PER_REG = 16;
// static const size_t NBITS_PER_WORD = 8;
static const size_t NBYTES_PER_WORD = 1;
// In end-to-end mode, we start high (255) and go low (0). Factoring in
// a query profile involves unsigned saturating subtraction, so all the
// query profile elements should be expressed as a positive penalty rather
// than a negative score.
typedef uint8_t TCScore;
/**
* Build query profile look up tables for the read. The query profile look
* up table is organized as a 1D array indexed by [i][j] where i is the
* reference character in the current DP column (0=A, 1=C, etc), and j is
* the segment of the query we're currently working on.
*/
void SwAligner::buildQueryProfileEnd2EndSseU8(bool fw) {
bool& done = fw ? sseU8fwBuilt_ : sseU8rcBuilt_;
if(done) {
return;
}
done = true;
const BTDnaString* rd = fw ? rdfw_ : rdrc_;
const BTString* qu = fw ? qufw_ : qurc_;
// daehwan - allows to align a portion of a read, not the whole.
// const size_t len = rd->length();
const size_t len = dpRows();
const size_t seglen = (len + (NWORDS_PER_REG-1)) / NWORDS_PER_REG;
// How many __m128i's are needed
size_t n128s =
64 + // slack bytes, for alignment?
(seglen * ALPHA_SIZE) // query profile data
* 2; // & gap barrier data
assert_gt(n128s, 0);
SSEData& d = fw ? sseU8fw_ : sseU8rc_;
d.profbuf_.resizeNoCopy(n128s);
assert(!d.profbuf_.empty());
d.maxPen_ = d.maxBonus_ = 0;
d.lastIter_ = d.lastWord_ = 0;
d.qprofStride_ = d.gbarStride_ = 2;
d.bias_ = 0; // no bias needed for end-to-end alignment; just use subtraction
// For each reference character A, C, G, T, N ...
for(size_t refc = 0; refc < ALPHA_SIZE; refc++) {
// For each segment ...
for(size_t i = 0; i < seglen; i++) {
size_t j = i;
uint8_t *qprofWords =
reinterpret_cast<uint8_t*>(d.profbuf_.ptr() + (refc * seglen * 2) + (i * 2));
uint8_t *gbarWords =
reinterpret_cast<uint8_t*>(d.profbuf_.ptr() + (refc * seglen * 2) + (i * 2) + 1);
// For each sub-word (byte) ...
for(size_t k = 0; k < NWORDS_PER_REG; k++) {
int sc = 0;
*gbarWords = 0;
if(j < len) {
int readc = (*rd)[j];
int readq = (*qu)[j];
sc = sc_->score(readc, (int)(1 << refc), readq - 33);
// Make score positive, to fit in an unsigned
sc = -sc;
assert_range(0, 255, sc);
size_t j_from_end = len - j - 1;
if(j < (size_t)sc_->gapbar ||
j_from_end < (size_t)sc_->gapbar)
{
// Inside the gap barrier
*gbarWords = 0xff;
}
}
if(refc == 0 && j == len-1) {
// Remember which 128-bit word and which smaller word has
// the final row
d.lastIter_ = i;
d.lastWord_ = k;
}
if((size_t)sc > d.maxPen_) {
d.maxPen_ = (size_t)sc;
}
*qprofWords = (uint8_t)sc;
gbarWords++;
qprofWords++;
j += seglen; // update offset into query
}
}
}
}
#ifndef NDEBUG
/**
* Return true iff the cell has sane E/F/H values w/r/t its predecessors.
*/
static bool cellOkEnd2EndU8(
SSEData& d,
size_t row,
size_t col,
int refc,
int readc,
int readq,
const Scoring& sc) // scoring scheme
{
TCScore floorsc = 0;
TAlScore ceilsc = MAX_I64;
TAlScore offsetsc = -0xff;
TAlScore sc_h_cur = (TAlScore)d.mat_.helt(row, col);
TAlScore sc_e_cur = (TAlScore)d.mat_.eelt(row, col);
TAlScore sc_f_cur = (TAlScore)d.mat_.felt(row, col);
if(sc_h_cur > floorsc) {
sc_h_cur += offsetsc;
}
if(sc_e_cur > floorsc) {
sc_e_cur += offsetsc;
}
if(sc_f_cur > floorsc) {
sc_f_cur += offsetsc;
}
bool gapsAllowed = true;
size_t rowFromEnd = d.mat_.nrow() - row - 1;
if(row < (size_t)sc.gapbar || rowFromEnd < (size_t)sc.gapbar) {
gapsAllowed = false;
}
bool e_left_trans = false, h_left_trans = false;
bool f_up_trans = false, h_up_trans = false;
bool h_diag_trans = false;
if(gapsAllowed) {
TAlScore sc_h_left = floorsc;
TAlScore sc_e_left = floorsc;
TAlScore sc_h_up = floorsc;
TAlScore sc_f_up = floorsc;
if(col > 0 && sc_e_cur > floorsc && sc_e_cur <= ceilsc) {
sc_h_left = d.mat_.helt(row, col-1) + offsetsc;
sc_e_left = d.mat_.eelt(row, col-1) + offsetsc;
e_left_trans = (sc_e_left > floorsc && sc_e_cur == sc_e_left - sc.readGapExtend());
h_left_trans = (sc_h_left > floorsc && sc_e_cur == sc_h_left - sc.readGapOpen());
assert(e_left_trans || h_left_trans);
// Check that we couldn't have got a better E score
assert_geq(sc_e_cur, sc_e_left - sc.readGapExtend());
assert_geq(sc_e_cur, sc_h_left - sc.readGapOpen());
}
if(row > 0 && sc_f_cur > floorsc && sc_f_cur <= ceilsc) {
sc_h_up = d.mat_.helt(row-1, col) + offsetsc;
sc_f_up = d.mat_.felt(row-1, col) + offsetsc;
f_up_trans = (sc_f_up > floorsc && sc_f_cur == sc_f_up - sc.refGapExtend());
h_up_trans = (sc_h_up > floorsc && sc_f_cur == sc_h_up - sc.refGapOpen());
assert(f_up_trans || h_up_trans);
// Check that we couldn't have got a better F score
assert_geq(sc_f_cur, sc_f_up - sc.refGapExtend());
assert_geq(sc_f_cur, sc_h_up - sc.refGapOpen());
}
} else {
assert_geq(floorsc, sc_e_cur);
assert_geq(floorsc, sc_f_cur);
}
if(col > 0 && row > 0 && sc_h_cur > floorsc && sc_h_cur <= ceilsc) {
TAlScore sc_h_upleft = d.mat_.helt(row-1, col-1) + offsetsc;
TAlScore sc_diag = sc.score(readc, (int)refc, readq - 33);
h_diag_trans = sc_h_cur == sc_h_upleft + sc_diag;
}
assert(
sc_h_cur <= floorsc ||
e_left_trans ||
h_left_trans ||
f_up_trans ||
h_up_trans ||
h_diag_trans ||
sc_h_cur > ceilsc ||
row == 0 ||
col == 0);
return true;
}
#endif /*ndef NDEBUG*/
#ifdef NDEBUG
#define assert_all_eq0(x)
#define assert_all_gt(x, y)
#define assert_all_gt_lo(x)
#define assert_all_lt(x, y)
#define assert_all_lt_hi(x)
#else
#define assert_all_eq0(x) { \
__m128i z = _mm_setzero_si128(); \
__m128i tmp = _mm_setzero_si128(); \
z = _mm_xor_si128(z, z); \
tmp = _mm_cmpeq_epi16(x, z); \
assert_eq(0xffff, _mm_movemask_epi8(tmp)); \
}
#define assert_all_gt(x, y) { \
__m128i tmp = _mm_cmpgt_epu8(x, y); \
assert_eq(0xffff, _mm_movemask_epi8(tmp)); \
}
#define assert_all_gt_lo(x) { \
__m128i z = _mm_setzero_si128(); \
__m128i tmp = _mm_setzero_si128(); \
z = _mm_xor_si128(z, z); \
tmp = _mm_cmpgt_epu8(x, z); \
assert_eq(0xffff, _mm_movemask_epi8(tmp)); \
}
#define assert_all_lt(x, y) { \
__m128i z = _mm_setzero_si128(); \
z = _mm_xor_si128(z, z); \
__m128i tmp = _mm_subs_epu8(y, x); \
tmp = _mm_cmpeq_epi16(tmp, z); \
assert_eq(0x0000, _mm_movemask_epi8(tmp)); \
}
#define assert_all_lt_hi(x) { \
__m128i z = _mm_setzero_si128(); \
__m128i tmp = _mm_setzero_si128(); \
z = _mm_cmpeq_epu8(z, z); \
z = _mm_srli_epu8(z, 1); \
tmp = _mm_cmplt_epu8(x, z); \
assert_eq(0xffff, _mm_movemask_epi8(tmp)); \
}
#endif
/**
* Aligns by filling a dynamic programming matrix with the SSE-accelerated,
* banded DP approach of Farrar. As it goes, it determines which cells we
* might backtrace from and tallies the best (highest-scoring) N backtrace
* candidate cells per diagonal. Also returns the alignment score of the best
* alignment in the matrix.
*
* This routine does *not* maintain a matrix holding the entire matrix worth of
* scores, nor does it maintain any other dense O(mn) data structure, as this
* would quickly exhaust memory for queries longer than about 10,000 kb.
* Instead, in the fill stage it maintains two columns worth of scores at a
* time (current/previous, or right/left) - these take O(m) space. When
* finished with the current column, it determines which cells from the
* previous column, if any, are candidates we might backtrace from to find a
* full alignment. A candidate cell has a score that rises above the threshold
* and isn't improved upon by a match in the next column. The best N
* candidates per diagonal are stored in a O(m + n) data structure.
*/
TAlScore SwAligner::alignGatherEE8(int& flag, bool debug) {
assert_leq(rdf_, rd_->length());
assert_leq(rdf_, qu_->length());
assert_lt(rfi_, rff_);
assert_lt(rdi_, rdf_);
assert_eq(rd_->length(), qu_->length());
assert_geq(sc_->gapbar, 1);
assert(repOk());
#ifndef NDEBUG
for(size_t i = (size_t)rfi_; i < (size_t)rff_; i++) {
assert_range(0, 16, (int)rf_[i]);
}
#endif
SSEData& d = fw_ ? sseU8fw_ : sseU8rc_;
SSEMetrics& met = extend_ ? sseU8ExtendMet_ : sseU8MateMet_;
if(!debug) met.dp++;
buildQueryProfileEnd2EndSseU8(fw_);
assert(!d.profbuf_.empty());
assert_eq(0, d.maxBonus_);
size_t iter =
(dpRows() + (NWORDS_PER_REG-1)) / NWORDS_PER_REG; // iter = segLen
int dup;
// Now set up the score vectors. We just need two columns worth, which
// we'll call "left" and "right".
d.vecbuf_.resize(4 * 2 * iter);
d.vecbuf_.zero();
__m128i *vbuf_l = d.vecbuf_.ptr();
__m128i *vbuf_r = d.vecbuf_.ptr() + (4 * iter);
// This is the data structure that holds candidate cells per diagonal.
const size_t ndiags = rff_ - rfi_ + dpRows() - 1;
if(!debug) {
btdiag_.init(ndiags, 2);
}
// Data structure that holds checkpointed anti-diagonals
TAlScore perfectScore = sc_->perfectScore(dpRows());
bool checkpoint = true;
bool cpdebug = false;
#ifndef NDEBUG
cpdebug = dpRows() < 1000;
#endif
cper_.init(
dpRows(), // # rows
rff_ - rfi_, // # columns
cperPerPow2_, // checkpoint every 1 << perpow2 diags (& next)
perfectScore, // perfect score (for sanity checks)
true, // matrix cells have 8-bit scores?
cperTri_, // triangular mini-fills?
false, // alignment is local?
cpdebug); // save all cells for debugging?
// Many thanks to Michael Farrar for releasing his striped Smith-Waterman
// implementation:
//
// http://sites.google.com/site/farrarmichael/smith-waterman
//
// Much of the implmentation below is adapted from Michael's code.
// Set all elts to reference gap open penalty
__m128i rfgapo = _mm_setzero_si128();
__m128i rfgape = _mm_setzero_si128();
__m128i rdgapo = _mm_setzero_si128();
__m128i rdgape = _mm_setzero_si128();
__m128i vlo = _mm_setzero_si128();
__m128i vhi = _mm_setzero_si128();
__m128i ve = _mm_setzero_si128();
__m128i vf = _mm_setzero_si128();
__m128i vh = _mm_setzero_si128();
__m128i vhd = _mm_setzero_si128();
__m128i vhdtmp = _mm_setzero_si128();
__m128i vtmp = _mm_setzero_si128();
__m128i vzero = _mm_setzero_si128();
__m128i vhilsw = _mm_setzero_si128();
assert_gt(sc_->refGapOpen(), 0);
assert_leq(sc_->refGapOpen(), MAX_U8);
dup = (sc_->refGapOpen() << 8) | (sc_->refGapOpen() & 0x00ff);
rfgapo = _mm_insert_epi16(rfgapo, dup, 0);
rfgapo = _mm_shufflelo_epi16(rfgapo, 0);
rfgapo = _mm_shuffle_epi32(rfgapo, 0);
// Set all elts to reference gap extension penalty
assert_gt(sc_->refGapExtend(), 0);
assert_leq(sc_->refGapExtend(), MAX_U8);
assert_leq(sc_->refGapExtend(), sc_->refGapOpen());
dup = (sc_->refGapExtend() << 8) | (sc_->refGapExtend() & 0x00ff);
rfgape = _mm_insert_epi16(rfgape, dup, 0);
rfgape = _mm_shufflelo_epi16(rfgape, 0);
rfgape = _mm_shuffle_epi32(rfgape, 0);
// Set all elts to read gap open penalty
assert_gt(sc_->readGapOpen(), 0);
assert_leq(sc_->readGapOpen(), MAX_U8);
dup = (sc_->readGapOpen() << 8) | (sc_->readGapOpen() & 0x00ff);
rdgapo = _mm_insert_epi16(rdgapo, dup, 0);
rdgapo = _mm_shufflelo_epi16(rdgapo, 0);
rdgapo = _mm_shuffle_epi32(rdgapo, 0);
// Set all elts to read gap extension penalty
assert_gt(sc_->readGapExtend(), 0);
assert_leq(sc_->readGapExtend(), MAX_U8);
assert_leq(sc_->readGapExtend(), sc_->readGapOpen());
dup = (sc_->readGapExtend() << 8) | (sc_->readGapExtend() & 0x00ff);
rdgape = _mm_insert_epi16(rdgape, dup, 0);
rdgape = _mm_shufflelo_epi16(rdgape, 0);
rdgape = _mm_shuffle_epi32(rdgape, 0);
vhi = _mm_cmpeq_epi16(vhi, vhi); // all elts = 0xffff
vlo = _mm_xor_si128(vlo, vlo); // all elts = 0
// vhilsw: topmost (least sig) word set to 0x7fff, all other words=0
vhilsw = _mm_shuffle_epi32(vhi, 0);
vhilsw = _mm_srli_si128(vhilsw, NBYTES_PER_REG - NBYTES_PER_WORD);
// Points to a long vector of __m128i where each element is a block of
// contiguous cells in the E, F or H matrix. If the index % 3 == 0, then
// the block of cells is from the E matrix. If index % 3 == 1, they're
// from the F matrix. If index % 3 == 2, then they're from the H matrix.
// Blocks of cells are organized in the same interleaved manner as they are
// calculated by the Farrar algorithm.
const __m128i *pvScore; // points into the query profile
const size_t colstride = ROWSTRIDE_2COL * iter;
// Initialize the H and E vectors in the first matrix column
__m128i *pvELeft = vbuf_l + 0; __m128i *pvERight = vbuf_r + 0;
/* __m128i *pvFLeft = vbuf_l + 1; */ __m128i *pvFRight = vbuf_r + 1;
__m128i *pvHLeft = vbuf_l + 2; __m128i *pvHRight = vbuf_r + 2;
// Maximum score in final row
bool found = false;
TCScore lrmax = MIN_U8;
for(size_t i = 0; i < iter; i++) {
_mm_store_si128(pvERight, vlo); pvERight += ROWSTRIDE_2COL;
// Could initialize Hs to high or low. If high, cells in the lower
// triangle will have somewhat more legitiate scores, but still won't
// be exhaustively scored.
_mm_store_si128(pvHRight, vlo); pvHRight += ROWSTRIDE_2COL;
}
assert_gt(sc_->gapbar, 0);
size_t nfixup = 0;
// Fill in the table as usual but instead of using the same gap-penalty
// vector for each iteration of the inner loop, load words out of a
// pre-calculated gap vector parallel to the query profile. The pre-
// calculated gap vectors enforce the gap barrier constraint by making it
// infinitely costly to introduce a gap in barrier rows.
//
// AND use a separate loop to fill in the first row of the table, enforcing
// the st_ constraints in the process. This is awkward because it
// separates the processing of the first row from the others and might make
// it difficult to use the first-row results in the next row, but it might
// be the simplest and least disruptive way to deal with the st_ constraint.
for(size_t i = (size_t)rfi_; i < (size_t)rff_; i++) {
// Swap left and right; vbuf_l is the vector on the left, which we
// generally load from, and vbuf_r is the vector on the right, which we
// generally store to.
swap(vbuf_l, vbuf_r);
pvELeft = vbuf_l + 0; pvERight = vbuf_r + 0;
/* pvFLeft = vbuf_l + 1; */ pvFRight = vbuf_r + 1;
pvHLeft = vbuf_l + 2; pvHRight = vbuf_r + 2;
// Fetch the appropriate query profile. Note that elements of rf_ must
// be numbers, not masks.
const int refc = (int)rf_[i];
// Fetch the appropriate query profile
size_t off = (size_t)firsts5[refc] * iter * 2;
pvScore = d.profbuf_.ptr() + off; // even elts = query profile, odd = gap barrier
// Set all cells to low value
vf = _mm_xor_si128(vf, vf);
// Load H vector from the final row of the previous column
vh = _mm_load_si128(pvHLeft + colstride - ROWSTRIDE_2COL);
// Shift 2 bytes down so that topmost (least sig) cell gets 0
vh = _mm_slli_si128(vh, NBYTES_PER_WORD);
// Fill topmost (least sig) cell with high value
vh = _mm_or_si128(vh, vhilsw);
// For each character in the reference text:
size_t j;
for(j = 0; j < iter; j++) {
// Load cells from E, calculated previously
ve = _mm_load_si128(pvELeft);
vhd = _mm_load_si128(pvHLeft);
assert_all_lt(ve, vhi);
pvELeft += ROWSTRIDE_2COL;
// Store cells in F, calculated previously
vf = _mm_subs_epu8(vf, pvScore[1]); // veto some ref gap extensions
_mm_store_si128(pvFRight, vf);
pvFRight += ROWSTRIDE_2COL;
// Factor in query profile (matches and mismatches)
vh = _mm_subs_epu8(vh, pvScore[0]);
// Update H, factoring in E and F
vh = _mm_max_epu8(vh, vf);
// Update vE value
vhdtmp = vhd;
vhd = _mm_subs_epu8(vhd, rdgapo);
vhd = _mm_subs_epu8(vhd, pvScore[1]); // veto some read gap opens
ve = _mm_subs_epu8(ve, rdgape);
ve = _mm_max_epu8(ve, vhd);
vh = _mm_max_epu8(vh, ve);
// Save the new vH values
_mm_store_si128(pvHRight, vh);
pvHRight += ROWSTRIDE_2COL;
vtmp = vh;
assert_all_lt(ve, vhi);
// Load the next h value
vh = vhdtmp;
pvHLeft += ROWSTRIDE_2COL;
// Save E values
_mm_store_si128(pvERight, ve);
pvERight += ROWSTRIDE_2COL;
// Update vf value
vtmp = _mm_subs_epu8(vtmp, rfgapo);
vf = _mm_subs_epu8(vf, rfgape);
assert_all_lt(vf, vhi);
vf = _mm_max_epu8(vf, vtmp);
pvScore += 2; // move on to next query profile / gap veto
}
// pvHStore, pvELoad, pvEStore have all rolled over to the next column
pvFRight -= colstride; // reset to start of column
vtmp = _mm_load_si128(pvFRight);
pvHRight -= colstride; // reset to start of column
vh = _mm_load_si128(pvHRight);
pvScore = d.profbuf_.ptr() + off + 1; // reset veto vector
// vf from last row gets shifted down by one to overlay the first row
// rfgape has already been subtracted from it.
vf = _mm_slli_si128(vf, NBYTES_PER_WORD);
vf = _mm_subs_epu8(vf, *pvScore); // veto some ref gap extensions
vf = _mm_max_epu8(vtmp, vf);
vtmp = _mm_subs_epu8(vf, vtmp);
vtmp = _mm_cmpeq_epi8(vtmp, vzero);
int cmp = _mm_movemask_epi8(vtmp);
// If any element of vtmp is greater than H - gap-open...
j = 0;
while(cmp != 0xffff) {
// Store this vf
_mm_store_si128(pvFRight, vf);
pvFRight += ROWSTRIDE_2COL;
// Update vh w/r/t new vf
vh = _mm_max_epu8(vh, vf);
// Save vH values
_mm_store_si128(pvHRight, vh);
pvHRight += ROWSTRIDE_2COL;
pvScore += 2;
assert_lt(j, iter);
if(++j == iter) {
pvFRight -= colstride;
vtmp = _mm_load_si128(pvFRight); // load next vf ASAP
pvHRight -= colstride;
vh = _mm_load_si128(pvHRight); // load next vh ASAP
pvScore = d.profbuf_.ptr() + off + 1;
j = 0;
vf = _mm_slli_si128(vf, NBYTES_PER_WORD);
} else {
vtmp = _mm_load_si128(pvFRight); // load next vf ASAP
vh = _mm_load_si128(pvHRight); // load next vh ASAP
}
// Update F with another gap extension
vf = _mm_subs_epu8(vf, rfgape);
vf = _mm_subs_epu8(vf, *pvScore); // veto some ref gap extensions
vf = _mm_max_epu8(vtmp, vf);
vtmp = _mm_subs_epu8(vf, vtmp);
vtmp = _mm_cmpeq_epi8(vtmp, vzero);
cmp = _mm_movemask_epi8(vtmp);
nfixup++;
}
// Check in the last row for the maximum so far
__m128i *vtmp = vbuf_r + 2 /* H */ + (d.lastIter_ * ROWSTRIDE_2COL);
// Note: we may not want to extract from the final row
TCScore lr = ((TCScore*)(vtmp))[d.lastWord_];
found = true;
if(lr > lrmax) {
lrmax = lr;
}
// Now we'd like to know whether the bottommost element of the right
// column is a candidate we might backtrace from. First question is:
// did it exceed the minimum score threshold?
TAlScore score = (TAlScore)(lr - 0xff);
if(lr == MIN_U8) {
score = MIN_I64;
}
if(!debug && score >= minsc_) {
DpBtCandidate cand(dpRows() - 1, i - rfi_, score);
btdiag_.add(i - rfi_, cand);
}
// Save some elements to checkpoints
if(checkpoint) {
__m128i *pvE = vbuf_r + 0;
__m128i *pvF = vbuf_r + 1;
__m128i *pvH = vbuf_r + 2;
size_t coli = i - rfi_;
if(coli < cper_.locol_) cper_.locol_ = coli;
if(coli > cper_.hicol_) cper_.hicol_ = coli;
if(cperTri_) {
size_t rc_mod = coli & cper_.lomask_;
assert_lt(rc_mod, cper_.per_);
int64_t row = -rc_mod-1;
int64_t row_mod = row;
int64_t row_div = 0;
size_t idx = coli >> cper_.perpow2_;
size_t idxrow = idx * cper_.nrow_;
assert_eq(4, ROWSTRIDE_2COL);
bool done = false;
while(true) {
row += (cper_.per_ - 2);
row_mod += (cper_.per_ - 2);
for(size_t j = 0; j < 2; j++) {
row++;
row_mod++;
if(row >= 0 && (size_t)row < cper_.nrow_) {
// Update row divided by iter_ and mod iter_
while(row_mod >= (int64_t)iter) {
row_mod -= (int64_t)iter;
row_div++;
}
size_t delt = idxrow + row;
size_t vecoff = (row_mod << 6) + row_div;
assert_lt(row_div, 16);
int16_t h_sc = ((uint8_t*)pvH)[vecoff];
int16_t e_sc = ((uint8_t*)pvE)[vecoff];
int16_t f_sc = ((uint8_t*)pvF)[vecoff];
if(h_sc == 0) h_sc = MIN_I16;
else h_sc -= 0xff;
if(e_sc == 0) e_sc = MIN_I16;
else e_sc -= 0xff;
if(f_sc == 0) f_sc = MIN_I16;
else f_sc -= 0xff;
assert_leq(h_sc, cper_.perf_);
assert_leq(e_sc, cper_.perf_);
assert_leq(f_sc, cper_.perf_);
CpQuad *qdiags = ((j == 0) ? cper_.qdiag1s_.ptr() : cper_.qdiag2s_.ptr());
qdiags[delt].sc[0] = h_sc;
qdiags[delt].sc[1] = e_sc;
qdiags[delt].sc[2] = f_sc;
} // if(row >= 0 && row < nrow_)
else if(row >= 0 && (size_t)row >= cper_.nrow_) {
done = true;
break;
}
} // end of loop over anti-diags
if(done) {
break;
}
idx++;
idxrow += cper_.nrow_;
}
} else {
// If this is the first column, take this opportunity to
// pre-calculate the coordinates of the elements we're going to
// checkpoint.
if(coli == 0) {
size_t cpi = cper_.per_-1;
size_t cpimod = cper_.per_-1;
size_t cpidiv = 0;
cper_.commitMap_.clear();
while(cpi < cper_.nrow_) {
while(cpimod >= iter) {
cpimod -= iter;
cpidiv++;
}
size_t vecoff = (cpimod << 6) + cpidiv;
cper_.commitMap_.push_back(vecoff);
cpi += cper_.per_;
cpimod += cper_.per_;
}
}
// Save all the rows
size_t rowoff = 0;
size_t sz = cper_.commitMap_.size();
for(size_t i = 0; i < sz; i++, rowoff += cper_.ncol_) {
size_t vecoff = cper_.commitMap_[i];
int16_t h_sc = ((uint8_t*)pvH)[vecoff];
//int16_t e_sc = ((uint8_t*)pvE)[vecoff];
int16_t f_sc = ((uint8_t*)pvF)[vecoff];
if(h_sc == 0) h_sc = MIN_I16;
else h_sc -= 0xff;
//if(e_sc == 0) e_sc = MIN_I16;
//else e_sc -= 0xff;
if(f_sc == 0) f_sc = MIN_I16;
else f_sc -= 0xff;
assert_leq(h_sc, cper_.perf_);
//assert_leq(e_sc, cper_.perf_);
assert_leq(f_sc, cper_.perf_);
CpQuad& dst = cper_.qrows_[rowoff + coli];
dst.sc[0] = h_sc;
//dst.sc[1] = e_sc;
dst.sc[2] = f_sc;
}
// Is this a column we'd like to checkpoint?
if((coli & cper_.lomask_) == cper_.lomask_) {
// Save the column using memcpys
assert_gt(coli, 0);
size_t wordspercol = cper_.niter_ * ROWSTRIDE_2COL;
size_t coloff = (coli >> cper_.perpow2_) * wordspercol;
__m128i *dst = cper_.qcols_.ptr() + coloff;
memcpy(dst, vbuf_r, sizeof(__m128i) * wordspercol);
}
}
if(cper_.debug_) {
// Save the column using memcpys
size_t wordspercol = cper_.niter_ * ROWSTRIDE_2COL;
size_t coloff = coli * wordspercol;
__m128i *dst = cper_.qcolsD_.ptr() + coloff;
memcpy(dst, vbuf_r, sizeof(__m128i) * wordspercol);
}
}
}
// Update metrics
if(!debug) {
size_t ninner = (rff_ - rfi_) * iter;
met.col += (rff_ - rfi_); // DP columns
met.cell += (ninner * NWORDS_PER_REG); // DP cells
met.inner += ninner; // DP inner loop iters
met.fixup += nfixup; // DP fixup loop iters
}
flag = 0;
// Did we find a solution?
TAlScore score = MIN_I64;
if(!found) {
flag = -1; // no
if(!debug) met.dpfail++;
return MIN_I64;
} else {
score = (TAlScore)(lrmax - 0xff);
if(score < minsc_) {
flag = -1; // no
if(!debug) met.dpfail++;
return score;
}
}
// Could we have saturated?
if(lrmax == MIN_U8) {
flag = -2; // yes
if(!debug) met.dpsat++;
return MIN_I64;
}
// Now take all the backtrace candidates in the btdaig_ structure and
// dump them into the btncand_ array. They'll be sorted later.
if(!debug) {
btdiag_.dump(btncand_);
assert(!btncand_.empty());
}
// Return largest score
if(!debug) met.dpsucc++;
return score;
}
/**
* Solve the current alignment problem using SSE instructions that operate on 16
* unsigned 8-bit values packed into a single 128-bit register.
*/
TAlScore SwAligner::alignNucleotidesEnd2EndSseU8(int& flag, bool debug) {
assert_leq(rdf_, rd_->length());
assert_leq(rdf_, qu_->length());
assert_lt(rfi_, rff_);
assert_lt(rdi_, rdf_);
assert_eq(rd_->length(), qu_->length());
assert_geq(sc_->gapbar, 1);
assert(repOk());
#ifndef NDEBUG
for(size_t i = (size_t)rfi_; i < (size_t)rff_; i++) {
assert_range(0, 16, (int)rf_[i]);
}
#endif
SSEData& d = fw_ ? sseU8fw_ : sseU8rc_;
SSEMetrics& met = extend_ ? sseU8ExtendMet_ : sseU8MateMet_;
if(!debug) met.dp++;
buildQueryProfileEnd2EndSseU8(fw_);
assert(!d.profbuf_.empty());
assert_eq(0, d.maxBonus_);
size_t iter =
(dpRows() + (NWORDS_PER_REG-1)) / NWORDS_PER_REG; // iter = segLen
int dup;
// Many thanks to Michael Farrar for releasing his striped Smith-Waterman
// implementation:
//
// http://sites.google.com/site/farrarmichael/smith-waterman
//
// Much of the implmentation below is adapted from Michael's code.
// Set all elts to reference gap open penalty
__m128i rfgapo = _mm_setzero_si128();
__m128i rfgape = _mm_setzero_si128();
__m128i rdgapo = _mm_setzero_si128();
__m128i rdgape = _mm_setzero_si128();
__m128i vlo = _mm_setzero_si128();
__m128i vhi = _mm_setzero_si128();
__m128i ve = _mm_setzero_si128();
__m128i vf = _mm_setzero_si128();
__m128i vh = _mm_setzero_si128();
#if 0
__m128i vhd = _mm_setzero_si128();
__m128i vhdtmp = _mm_setzero_si128();
#endif
__m128i vtmp = _mm_setzero_si128();
__m128i vzero = _mm_setzero_si128();
__m128i vhilsw = _mm_setzero_si128();
assert_gt(sc_->refGapOpen(), 0);
assert_leq(sc_->refGapOpen(), MAX_U8);
dup = (sc_->refGapOpen() << 8) | (sc_->refGapOpen() & 0x00ff);
rfgapo = _mm_insert_epi16(rfgapo, dup, 0);
rfgapo = _mm_shufflelo_epi16(rfgapo, 0);
rfgapo = _mm_shuffle_epi32(rfgapo, 0);
// Set all elts to reference gap extension penalty
assert_gt(sc_->refGapExtend(), 0);
assert_leq(sc_->refGapExtend(), MAX_U8);
assert_leq(sc_->refGapExtend(), sc_->refGapOpen());
dup = (sc_->refGapExtend() << 8) | (sc_->refGapExtend() & 0x00ff);
rfgape = _mm_insert_epi16(rfgape, dup, 0);
rfgape = _mm_shufflelo_epi16(rfgape, 0);
rfgape = _mm_shuffle_epi32(rfgape, 0);
// Set all elts to read gap open penalty
assert_gt(sc_->readGapOpen(), 0);
assert_leq(sc_->readGapOpen(), MAX_U8);
dup = (sc_->readGapOpen() << 8) | (sc_->readGapOpen() & 0x00ff);
rdgapo = _mm_insert_epi16(rdgapo, dup, 0);
rdgapo = _mm_shufflelo_epi16(rdgapo, 0);
rdgapo = _mm_shuffle_epi32(rdgapo, 0);
// Set all elts to read gap extension penalty
assert_gt(sc_->readGapExtend(), 0);
assert_leq(sc_->readGapExtend(), MAX_U8);
assert_leq(sc_->readGapExtend(), sc_->readGapOpen());
dup = (sc_->readGapExtend() << 8) | (sc_->readGapExtend() & 0x00ff);
rdgape = _mm_insert_epi16(rdgape, dup, 0);
rdgape = _mm_shufflelo_epi16(rdgape, 0);
rdgape = _mm_shuffle_epi32(rdgape, 0);
vhi = _mm_cmpeq_epi16(vhi, vhi); // all elts = 0xffff
vlo = _mm_xor_si128(vlo, vlo); // all elts = 0
// vhilsw: topmost (least sig) word set to 0x7fff, all other words=0
vhilsw = _mm_shuffle_epi32(vhi, 0);
vhilsw = _mm_srli_si128(vhilsw, NBYTES_PER_REG - NBYTES_PER_WORD);
// Points to a long vector of __m128i where each element is a block of
// contiguous cells in the E, F or H matrix. If the index % 3 == 0, then
// the block of cells is from the E matrix. If index % 3 == 1, they're
// from the F matrix. If index % 3 == 2, then they're from the H matrix.
// Blocks of cells are organized in the same interleaved manner as they are
// calculated by the Farrar algorithm.
const __m128i *pvScore; // points into the query profile
d.mat_.init(dpRows(), rff_ - rfi_, NWORDS_PER_REG);
const size_t colstride = d.mat_.colstride();
//const size_t rowstride = d.mat_.rowstride();
assert_eq(ROWSTRIDE, colstride / iter);
// Initialize the H and E vectors in the first matrix column
__m128i *pvHTmp = d.mat_.tmpvec(0, 0);
__m128i *pvETmp = d.mat_.evec(0, 0);
// Maximum score in final row
bool found = false;
TCScore lrmax = MIN_U8;
for(size_t i = 0; i < iter; i++) {
_mm_store_si128(pvETmp, vlo);
_mm_store_si128(pvHTmp, vlo); // start high in end-to-end mode
pvETmp += ROWSTRIDE;
pvHTmp += ROWSTRIDE;
}
// These are swapped just before the innermost loop
__m128i *pvHStore = d.mat_.hvec(0, 0);
__m128i *pvHLoad = d.mat_.tmpvec(0, 0);
__m128i *pvELoad = d.mat_.evec(0, 0);
__m128i *pvEStore = d.mat_.evecUnsafe(0, 1);
__m128i *pvFStore = d.mat_.fvec(0, 0);
__m128i *pvFTmp = NULL;
assert_gt(sc_->gapbar, 0);
size_t nfixup = 0;
// Fill in the table as usual but instead of using the same gap-penalty
// vector for each iteration of the inner loop, load words out of a
// pre-calculated gap vector parallel to the query profile. The pre-
// calculated gap vectors enforce the gap barrier constraint by making it
// infinitely costly to introduce a gap in barrier rows.
//
// AND use a separate loop to fill in the first row of the table, enforcing
// the st_ constraints in the process. This is awkward because it
// separates the processing of the first row from the others and might make
// it difficult to use the first-row results in the next row, but it might
// be the simplest and least disruptive way to deal with the st_ constraint.
colstop_ = rff_ - 1;
lastsolcol_ = 0;
for(size_t i = (size_t)rfi_; i < (size_t)rff_; i++) {
assert(pvFStore == d.mat_.fvec(0, i - rfi_));
assert(pvHStore == d.mat_.hvec(0, i - rfi_));
// Fetch the appropriate query profile. Note that elements of rf_ must
// be numbers, not masks.
const int refc = (int)rf_[i];
size_t off = (size_t)firsts5[refc] * iter * 2;
pvScore = d.profbuf_.ptr() + off; // even elts = query profile, odd = gap barrier
// Set all cells to low value
vf = _mm_xor_si128(vf, vf);
// Load H vector from the final row of the previous column
vh = _mm_load_si128(pvHLoad + colstride - ROWSTRIDE);
// Shift 2 bytes down so that topmost (least sig) cell gets 0
vh = _mm_slli_si128(vh, NBYTES_PER_WORD);
// Fill topmost (least sig) cell with high value
vh = _mm_or_si128(vh, vhilsw);
// For each character in the reference text:
size_t j;
for(j = 0; j < iter; j++) {
// Load cells from E, calculated previously
ve = _mm_load_si128(pvELoad);
#if 0
vhd = _mm_load_si128(pvHLoad);
#endif
assert_all_lt(ve, vhi);
pvELoad += ROWSTRIDE;
// Store cells in F, calculated previously
vf = _mm_subs_epu8(vf, pvScore[1]); // veto some ref gap extensions
_mm_store_si128(pvFStore, vf);
pvFStore += ROWSTRIDE;
// Factor in query profile (matches and mismatches)
vh = _mm_subs_epu8(vh, pvScore[0]);
// Update H, factoring in E and F
vh = _mm_max_epu8(vh, ve);
vh = _mm_max_epu8(vh, vf);
// Save the new vH values
_mm_store_si128(pvHStore, vh);
pvHStore += ROWSTRIDE;
// Update vE value
vtmp = vh;
#if 0
vhdtmp = vhd;
vhd = _mm_subs_epu8(vhd, rdgapo);
vhd = _mm_subs_epu8(vhd, pvScore[1]); // veto some read gap opens
ve = _mm_subs_epu8(ve, rdgape);
ve = _mm_max_epu8(ve, vhd);
#else
vh = _mm_subs_epu8(vh, rdgapo);
vh = _mm_subs_epu8(vh, pvScore[1]); // veto some read gap opens
ve = _mm_subs_epu8(ve, rdgape);
ve = _mm_max_epu8(ve, vh);
#endif
assert_all_lt(ve, vhi);
// Load the next h value
#if 0
vh = vhdtmp;