1 #include "IncludeDefine.h"
2 #include "Parameters.h"
3 #include "ReadAlign.h"
4 #include "ErrorWarning.h"
5
assignAlignToWindow(uint a1,uint aLength,uint aStr,uint aNrep,uint aFrag,uint aRstart,bool aAnchor,uint sjA)6 void ReadAlign::assignAlignToWindow(uint a1, uint aLength, uint aStr, uint aNrep, uint aFrag, uint aRstart, bool aAnchor, uint sjA) {
7
8 uint iW=winBin[aStr][a1>>P.winBinNbits];
9
10 if (iW==uintWinBinMax || (!aAnchor && aLength < WALrec[iW]) ) return; //alignment does not belong to any window, or it's shorter than rec-length
11
12 //check if this alignment overlaps with any other alignment in the window, record the longest of the two
13 {//do not check for overlap if this is an sj-align
14 uint iA;
15 for (iA=0; iA<nWA[iW]; iA++) {
16 if (aFrag==WA[iW][iA][WA_iFrag] && WA[iW][iA][WA_sjA]==sjA \
17 && a1+WA[iW][iA][WA_rStart]==WA[iW][iA][WA_gStart]+aRstart \
18 && ( (aRstart>=WA[iW][iA][WA_rStart] && aRstart<WA[iW][iA][WA_rStart]+WA[iW][iA][WA_Length]) \
19 || (aRstart+aLength>=WA[iW][iA][WA_rStart] && aRstart+aLength<WA[iW][iA][WA_rStart]+WA[iW][iA][WA_Length]) ) ) {//this piece overlaps with iA
20 break;
21 };
22 };
23 if (iA<nWA[iW]) {//found overlap
24 if (aLength>WA[iW][iA][WA_Length]) {//replace
25
26 uint iA0;//iA0 is where the align has to be inserted
27 for (iA0=0;iA0<nWA[iW];iA0++)
28 {//find the insertion point TODO binary search
29 if (iA0!=iA && aRstart<WA[iW][iA0][WA_rStart])
30 {//do not compare with the piece to be removed
31 break;
32 };
33 };
34
35 if (iA0>iA)
36 {//true insertion place since iA will be removed
37 --iA0;
38 };
39
40 if (iA0<iA) {//shift aligns down
41 for (uint iA1=iA;iA1>iA0;iA1--) {//shift aligns to free up insertion point
42 for (uint ii=0;ii<WA_SIZE;ii++) {
43 WA[iW][iA1][ii]=WA[iW][iA1-1][ii];
44 };
45 };
46 } else if (iA0>iA) {//shift aligns up
47 for (uint iA1=iA;iA1<iA0;iA1++) {//shift aligns to free up insertion point
48 for (uint ii=0;ii<WA_SIZE;ii++) {
49 WA[iW][iA1][ii]=WA[iW][iA1+1][ii];
50 };
51 };
52 };
53
54
55 WA[iW][iA0][WA_rStart]=aRstart;
56 WA[iW][iA0][WA_Length]=aLength;
57 WA[iW][iA0][WA_gStart]=a1;
58 WA[iW][iA0][WA_Nrep]=aNrep;
59 WA[iW][iA0][WA_Anchor]=int(aAnchor);//=0 if not, =1 if yes
60 WA[iW][iA0][WA_iFrag]=aFrag;
61 WA[iW][iA0][WA_sjA]=sjA;
62
63 };
64 return; //do not record new align
65 };
66 };
67
68 if (nWA[iW]==P.seedPerWindowNmax) {//too many aligns per window, re-calcualte min-length, remove the shortest one,
69
70 WALrec[iW]=Lread+1;
71 for (uint iA=0; iA<nWA[iW]; iA++) {//find the new min-length
72 if (WA[iW][iA][WA_Anchor]!=1) WALrec[iW]=min(WALrec[iW],WA[iW][iA][WA_Length]); //protect the anchors - they are not counted for min-length
73 };
74
75
76 if (WALrec[iW]==Lread+1) {//this could happen if there are too many anchors
77 mapMarker=MARKER_TOO_MANY_ANCHORS_PER_WINDOW;
78 nW=0;
79 return;
80 };
81
82
83 if (!aAnchor && aLength < WALrec[iW]) return; //alignment is shorter than min-length, do not record - unless it's an anchor
84
85 uint iA1=0;
86 for (uint iA=0; iA<nWA[iW]; iA++) {//remove the shortest aligns
87 if ( WA[iW][iA][WA_Anchor]==1 || WA[iW][iA][WA_Length] > WALrec[iW] ) {//re-record the anchors and long aligns
88 for (uint ii=0; ii<WA_SIZE; ii++) WA[iW][iA1][ii]=WA[iW][iA][ii]; //re-record the iA-th alignment into iA1-th place
89 iA1++;
90 };
91 };
92
93 nWA[iW]=iA1;
94
95 if (!aAnchor && aLength <= WALrec[iW]) {//current align was removed, zero out its nWAP
96 nWAP[iW]=0;
97 };
98
99 };
100
101 if ( aAnchor || aLength > WALrec[iW] ) {
102 if (nWA[iW]>=P.seedPerWindowNmax) {
103 exitWithError("BUG: iA>=P.seedPerWindowNmax in stitchPieces, exiting",std::cerr, P.inOut->logMain, EXIT_CODE_BUG, P);
104 };
105
106 uint iA;
107 for (iA=0;iA<nWA[iW];iA++) {//find the insertion point in case aligns are not sorted by aRstart
108 //TODO binary search
109 if (aRstart<WA[iW][iA][WA_rStart]) break;
110 };
111 for (uint iA1=nWA[iW];iA1>iA;iA1--) {//shift aligns for free up insertion point
112 for (uint ii=0;ii<WA_SIZE;ii++) {
113 WA[iW][iA1][ii]=WA[iW][iA1-1][ii];
114 };
115 };
116
117 // now iW is the window to which this align belongs, record it
118 WA[iW][iA][WA_rStart]=aRstart;
119 WA[iW][iA][WA_Length]=aLength;
120 WA[iW][iA][WA_gStart]=a1;
121 WA[iW][iA][WA_Nrep]=aNrep;
122 WA[iW][iA][WA_Anchor]=int(aAnchor);//=0 if not, =1 if yes
123 WA[iW][iA][WA_iFrag]=aFrag;
124 WA[iW][iA][WA_sjA]=sjA;
125
126 nWA[iW]++;
127 nWAP[iW]++;
128 if (aAnchor && WlastAnchor[iW]<iA) WlastAnchor[iW]=iA; //record the index of the last anchor
129 };
130 };
131