jpayne@68
|
1 package aligner;
|
jpayne@68
|
2
|
jpayne@68
|
3 import dna.AminoAcid;
|
jpayne@68
|
4 import shared.KillSwitch;
|
jpayne@68
|
5 import shared.Tools;
|
jpayne@68
|
6
|
jpayne@68
|
7 /**
|
jpayne@68
|
8 * Based on MSA9PBA, but reduced to a single matrix. */
|
jpayne@68
|
9 public final class SingleStateAlignerFlat implements Aligner {
|
jpayne@68
|
10
|
jpayne@68
|
11
|
jpayne@68
|
12 public SingleStateAlignerFlat(){}
|
jpayne@68
|
13
|
jpayne@68
|
14 private void prefillTopRow(){
|
jpayne@68
|
15 final int[] header=packed[0];
|
jpayne@68
|
16 final int qlen=rows;
|
jpayne@68
|
17 for(int i=0; i<=columns; i++){
|
jpayne@68
|
18 int x=columns-i+1;
|
jpayne@68
|
19 int qbases=qlen-x;
|
jpayne@68
|
20
|
jpayne@68
|
21 //Minimal points to prefer a leftmost alignment
|
jpayne@68
|
22 header[i]=qbases<=0 ? 0 : ((-qbases)<<STARTOFFSET);
|
jpayne@68
|
23
|
jpayne@68
|
24 //Forces consumption of query, but does not allow for insertions...
|
jpayne@68
|
25 // header[i]=qbases<=0 ? 0 : calcDelScoreOffset(qbases);
|
jpayne@68
|
26 }
|
jpayne@68
|
27 }
|
jpayne@68
|
28
|
jpayne@68
|
29 private void prefillLeftColumnStartingAt(int i){
|
jpayne@68
|
30 packed[0][0]=MODE_MATCH;
|
jpayne@68
|
31 i=Tools.max(1, i);
|
jpayne@68
|
32 for(int score=MODE_INS+(POINTSoff_INS*i); i<=maxRows; i++){//Fill column 0 with insertions
|
jpayne@68
|
33 score+=POINTSoff_INS;
|
jpayne@68
|
34 packed[i][0]=score;
|
jpayne@68
|
35 }
|
jpayne@68
|
36 }
|
jpayne@68
|
37
|
jpayne@68
|
38 // private void clearInner(){//NotmaxRows needed
|
jpayne@68
|
39 // for(int i=1; i<=maxRows; i++){
|
jpayne@68
|
40 // for(int j=1; j<=maxColumns; j++){
|
jpayne@68
|
41 // packed[i][j]=0;
|
jpayne@68
|
42 // }
|
jpayne@68
|
43 // }
|
jpayne@68
|
44 // }
|
jpayne@68
|
45 //
|
jpayne@68
|
46 // private void clearAll(){//NotmaxRows needed
|
jpayne@68
|
47 // for(int i=0; i<=maxRows; i++){
|
jpayne@68
|
48 // for(int j=0; j<=maxColumns; j++){
|
jpayne@68
|
49 // packed[i][j]=0;
|
jpayne@68
|
50 // }
|
jpayne@68
|
51 // }
|
jpayne@68
|
52 // }
|
jpayne@68
|
53
|
jpayne@68
|
54 private void initialize(int rows_, int columns_){
|
jpayne@68
|
55 rows=rows_;
|
jpayne@68
|
56 columns=columns_;
|
jpayne@68
|
57 if(rows<=maxRows && columns<=maxColumns){
|
jpayne@68
|
58 prefillTopRow();
|
jpayne@68
|
59 // prefillLeftColumn();
|
jpayne@68
|
60 return;
|
jpayne@68
|
61 }
|
jpayne@68
|
62
|
jpayne@68
|
63 final int maxRows0=maxRows;
|
jpayne@68
|
64 final int maxColumns0=maxColumns;
|
jpayne@68
|
65 final int[][] packed0=packed;
|
jpayne@68
|
66
|
jpayne@68
|
67 //Monotonic increase
|
jpayne@68
|
68 maxRows=Tools.max(maxRows, rows+10);
|
jpayne@68
|
69 maxColumns=Tools.max(maxColumns, columns+10);
|
jpayne@68
|
70
|
jpayne@68
|
71 if(packed==null || maxColumns>maxColumns0){//Make a new matrix
|
jpayne@68
|
72 packed=KillSwitch.allocInt2D(maxRows+1, maxColumns+1);
|
jpayne@68
|
73 prefillLeftColumnStartingAt(1);
|
jpayne@68
|
74 }else{//Copy old rows
|
jpayne@68
|
75 assert(maxRows0>0 && maxColumns0>0);
|
jpayne@68
|
76 assert(maxRows>maxRows0 && maxColumns<=maxColumns0);
|
jpayne@68
|
77 packed=KillSwitch.allocInt2D(maxRows+1);
|
jpayne@68
|
78 for(int i=0; i<packed.length; i++){
|
jpayne@68
|
79 if(i<packed0.length){
|
jpayne@68
|
80 packed[i]=packed0[i];
|
jpayne@68
|
81 }else{
|
jpayne@68
|
82 packed[i]=KillSwitch.allocInt1D(maxColumns+1);
|
jpayne@68
|
83 }
|
jpayne@68
|
84 }
|
jpayne@68
|
85 //Fill column 0 with insertions
|
jpayne@68
|
86 prefillLeftColumnStartingAt(maxRows0);
|
jpayne@68
|
87 }
|
jpayne@68
|
88 prefillTopRow();
|
jpayne@68
|
89 }
|
jpayne@68
|
90
|
jpayne@68
|
91 /** return new int[] {rows, maxCol, maxState, maxScore, maxStart};
|
jpayne@68
|
92 * Will not fill areas that cannot match minScore */
|
jpayne@68
|
93 @Override
|
jpayne@68
|
94 public final int[] fillLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
|
jpayne@68
|
95 return fillUnlimited(read, ref, refStartLoc, refEndLoc, minScore);
|
jpayne@68
|
96 }
|
jpayne@68
|
97
|
jpayne@68
|
98 @Override
|
jpayne@68
|
99 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc){
|
jpayne@68
|
100 return fillUnlimited(read, ref, refStartLoc, refEndLoc, -999999);
|
jpayne@68
|
101 }
|
jpayne@68
|
102
|
jpayne@68
|
103 /** return new int[] {rows, maxCol, maxState, maxScore, maxStart};
|
jpayne@68
|
104 * Min score is optional */
|
jpayne@68
|
105 @Override
|
jpayne@68
|
106 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
|
jpayne@68
|
107 initialize(read.length, refEndLoc-refStartLoc+1);
|
jpayne@68
|
108
|
jpayne@68
|
109 //temporary, for finding a bug
|
jpayne@68
|
110 if(rows>maxRows || columns>maxColumns){
|
jpayne@68
|
111 throw new RuntimeException("rows="+rows+", maxRows="+maxRows+", cols="+columns+", maxCols="+maxColumns+"\n"+new String(read)+"\n");
|
jpayne@68
|
112 }
|
jpayne@68
|
113
|
jpayne@68
|
114 assert(rows<=maxRows) : "Check that values are in-bounds before calling this function: "+rows+", "+maxRows;
|
jpayne@68
|
115 assert(columns<=maxColumns) : "Check that values are in-bounds before calling this function: "+columns+", "+maxColumns;
|
jpayne@68
|
116
|
jpayne@68
|
117 assert(refStartLoc>=0) : "Check that values are in-bounds before calling this function: "+refStartLoc;
|
jpayne@68
|
118 assert(refEndLoc<ref.length) : "Check that values are in-bounds before calling this function: "+refEndLoc+", "+ref.length;
|
jpayne@68
|
119
|
jpayne@68
|
120 final int refOffset=refStartLoc-1;
|
jpayne@68
|
121 for(int row=1; row<=rows; row++){
|
jpayne@68
|
122
|
jpayne@68
|
123 final byte qBase=read[row-1];
|
jpayne@68
|
124 for(int col=1; col<=columns; col++){
|
jpayne@68
|
125
|
jpayne@68
|
126 final byte rBase=ref[refOffset+col];
|
jpayne@68
|
127
|
jpayne@68
|
128 final boolean match=(qBase==rBase);
|
jpayne@68
|
129 final boolean defined=(qBase!='N' && rBase!='N');
|
jpayne@68
|
130
|
jpayne@68
|
131 final int valueFromDiag=packed[row-1][col-1];
|
jpayne@68
|
132 final int valueFromDel=packed[row][col-1];
|
jpayne@68
|
133 final int valueFromIns=packed[row-1][col];
|
jpayne@68
|
134 // final int stateFromDiag=valueFromDiag&MODEMASK;
|
jpayne@68
|
135 final int scoreFromDiag=valueFromDiag&HIGHMASK;
|
jpayne@68
|
136 // final int stateFromDel=valueFromDel&MODEMASK;
|
jpayne@68
|
137 final int scoreFromDel=valueFromDel&HIGHMASK;
|
jpayne@68
|
138 // final int stateFromIns=valueFromIns&MODEMASK;
|
jpayne@68
|
139 final int scoreFromIns=valueFromIns&HIGHMASK;
|
jpayne@68
|
140
|
jpayne@68
|
141 // final boolean prevMatch=stateFromDiag==MODE_MATCH;
|
jpayne@68
|
142 // final boolean prevSub=stateFromDiag==MODE_SUB;
|
jpayne@68
|
143 // final boolean prevDel=stateFromDel==MODE_DEL;
|
jpayne@68
|
144 // final boolean prevIns=stateFromIns==MODE_INS;
|
jpayne@68
|
145
|
jpayne@68
|
146 //Old conditional code, replaced by faster arrays
|
jpayne@68
|
147 // final int diagScoreM=MODE_MATCH|(prevMatch ? POINTSoff_MATCH2 : POINTSoff_MATCH);
|
jpayne@68
|
148 // final int diagScoreS=MODE_SUB|(prevSub ? POINTSoff_SUB2 : POINTSoff_SUB);
|
jpayne@68
|
149 // final int delScore=scoreFromDel+(prevDel ? POINTSoff_DEL2 : POINTSoff_DEL)|MODE_DEL;
|
jpayne@68
|
150 // final int insScore=scoreFromIns+(prevIns ? POINTSoff_INS2 : POINTSoff_INS)|MODE_INS;
|
jpayne@68
|
151
|
jpayne@68
|
152 final int diagScoreM=POINTSoff_MATCH|MODE_MATCH;
|
jpayne@68
|
153 final int diagScoreS=POINTSoff_SUB|MODE_SUB;
|
jpayne@68
|
154 final int delScore=scoreFromDel+POINTSoff_DEL|MODE_DEL;
|
jpayne@68
|
155 final int insScore=scoreFromIns+POINTSoff_INS|MODE_INS;
|
jpayne@68
|
156
|
jpayne@68
|
157 int diagScore=(match ? diagScoreM : diagScoreS);
|
jpayne@68
|
158 diagScore=scoreFromDiag+(defined ? diagScore : POINTSoff_NOREF_MODE_SUB);
|
jpayne@68
|
159
|
jpayne@68
|
160 int score=diagScore>=delScore ? diagScore : delScore;
|
jpayne@68
|
161 score=score>=insScore ? score : insScore;
|
jpayne@68
|
162
|
jpayne@68
|
163 packed[row][col]=score;
|
jpayne@68
|
164 }
|
jpayne@68
|
165 // iterationsUnlimited+=columns;
|
jpayne@68
|
166 }
|
jpayne@68
|
167
|
jpayne@68
|
168
|
jpayne@68
|
169 int maxCol=-1;
|
jpayne@68
|
170 int maxState=-1;
|
jpayne@68
|
171 int maxStart=-1;
|
jpayne@68
|
172 int maxScore=Integer.MIN_VALUE;
|
jpayne@68
|
173
|
jpayne@68
|
174 for(int col=1; col<=columns; col++){
|
jpayne@68
|
175 int x=packed[rows][col];
|
jpayne@68
|
176 if((x&SCOREMASK)>maxScore){
|
jpayne@68
|
177 maxScore=x;
|
jpayne@68
|
178 maxCol=col;
|
jpayne@68
|
179 maxState=x&MODEMASK;
|
jpayne@68
|
180 maxStart=x&STARTMASK;
|
jpayne@68
|
181 }
|
jpayne@68
|
182 }
|
jpayne@68
|
183 maxScore>>=SCOREOFFSET;
|
jpayne@68
|
184
|
jpayne@68
|
185 // System.err.println("Returning "+rows+", "+maxCol+", "+maxState+", "+maxScore+"; minScore="+minScore);
|
jpayne@68
|
186 return maxScore<minScore ? null : new int[] {rows, maxCol, maxState, maxScore, maxStart};
|
jpayne@68
|
187 }
|
jpayne@68
|
188
|
jpayne@68
|
189 /** Generates the match string */
|
jpayne@68
|
190 @Override
|
jpayne@68
|
191 public final byte[] traceback(byte[] query, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state){
|
jpayne@68
|
192 // assert(false);
|
jpayne@68
|
193 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
|
jpayne@68
|
194 assert(row==rows);
|
jpayne@68
|
195
|
jpayne@68
|
196 byte[] out=new byte[row+col-1]; //TODO if an out of bound crash occurs, try removing the "-1".
|
jpayne@68
|
197 int outPos=0;
|
jpayne@68
|
198
|
jpayne@68
|
199 // assert(state==(packed[row][col]&MODEMASK));
|
jpayne@68
|
200
|
jpayne@68
|
201 while(row>0 && col>0){
|
jpayne@68
|
202 byte q=query[row-1];
|
jpayne@68
|
203 byte r=ref[refStartLoc+col-1];
|
jpayne@68
|
204 boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
|
jpayne@68
|
205 state=(packed[row][col]&MODEMASK);
|
jpayne@68
|
206 if(state==MODE_MATCH){
|
jpayne@68
|
207 col--;
|
jpayne@68
|
208 row--;
|
jpayne@68
|
209 out[outPos]=defined ? (byte)'m' : (byte)'N';
|
jpayne@68
|
210 }else if(state==MODE_SUB){
|
jpayne@68
|
211 col--;
|
jpayne@68
|
212 row--;
|
jpayne@68
|
213 out[outPos]=defined ? (byte)'S' : (byte)'N';
|
jpayne@68
|
214 }else if(state==MODE_N){
|
jpayne@68
|
215 col--;
|
jpayne@68
|
216 row--;
|
jpayne@68
|
217 out[outPos]='N';
|
jpayne@68
|
218 }else if(state==MODE_DEL){
|
jpayne@68
|
219 col--;
|
jpayne@68
|
220 out[outPos]='D';
|
jpayne@68
|
221 }else if(state==MODE_INS){
|
jpayne@68
|
222 row--;
|
jpayne@68
|
223 // out[outPos]='I';
|
jpayne@68
|
224 if(col>=0 && col<columns){
|
jpayne@68
|
225 out[outPos]='I';
|
jpayne@68
|
226 }else{
|
jpayne@68
|
227 out[outPos]='C';
|
jpayne@68
|
228 col--;
|
jpayne@68
|
229 }
|
jpayne@68
|
230 }else{
|
jpayne@68
|
231 assert(false) : state;
|
jpayne@68
|
232 }
|
jpayne@68
|
233 outPos++;
|
jpayne@68
|
234 }
|
jpayne@68
|
235
|
jpayne@68
|
236 assert(row==0 || col==0);
|
jpayne@68
|
237 if(col!=row){//Not sure what this is doing
|
jpayne@68
|
238 while(row>0){
|
jpayne@68
|
239 out[outPos]='C';
|
jpayne@68
|
240 outPos++;
|
jpayne@68
|
241 row--;
|
jpayne@68
|
242 col--;
|
jpayne@68
|
243 }
|
jpayne@68
|
244 if(col>0){
|
jpayne@68
|
245 //do nothing
|
jpayne@68
|
246 }
|
jpayne@68
|
247 }
|
jpayne@68
|
248
|
jpayne@68
|
249 //Shrink and reverse the string
|
jpayne@68
|
250 byte[] out2=new byte[outPos];
|
jpayne@68
|
251 for(int i=0; i<outPos; i++){
|
jpayne@68
|
252 out2[i]=out[outPos-i-1];
|
jpayne@68
|
253 }
|
jpayne@68
|
254 out=null;
|
jpayne@68
|
255
|
jpayne@68
|
256 return out2;
|
jpayne@68
|
257 }
|
jpayne@68
|
258
|
jpayne@68
|
259 @Override
|
jpayne@68
|
260 /** Generates identity;
|
jpayne@68
|
261 * fills 'extra' with {match, sub, del, ins, N, clip} if present */
|
jpayne@68
|
262 public float tracebackIdentity(byte[] query, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state, int[] extra){
|
jpayne@68
|
263
|
jpayne@68
|
264 // assert(false);
|
jpayne@68
|
265 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
|
jpayne@68
|
266 assert(row==rows);
|
jpayne@68
|
267
|
jpayne@68
|
268 // assert(state==(packed[row][col]&MODEMASK));
|
jpayne@68
|
269 int match=0, sub=0, del=0, ins=0, noref=0, clip=0;
|
jpayne@68
|
270
|
jpayne@68
|
271 while(row>0 && col>0){
|
jpayne@68
|
272 byte q=query[row-1];
|
jpayne@68
|
273 byte r=ref[refStartLoc+col-1];
|
jpayne@68
|
274 boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
|
jpayne@68
|
275 state=(packed[row][col]&MODEMASK);
|
jpayne@68
|
276 if(state==MODE_MATCH){
|
jpayne@68
|
277 col--;
|
jpayne@68
|
278 row--;
|
jpayne@68
|
279 match+=(defined ? 1 : 0);
|
jpayne@68
|
280 noref+=(defined ? 0 : 1);
|
jpayne@68
|
281 }else if(state==MODE_SUB){
|
jpayne@68
|
282 col--;
|
jpayne@68
|
283 row--;
|
jpayne@68
|
284 sub+=(defined ? 1 : 0);
|
jpayne@68
|
285 noref+=(defined ? 0 : 1);
|
jpayne@68
|
286 }else if(state==MODE_N){//Can't currently happen
|
jpayne@68
|
287 col--;
|
jpayne@68
|
288 row--;
|
jpayne@68
|
289 noref++;
|
jpayne@68
|
290 }else if(state==MODE_DEL){
|
jpayne@68
|
291 col--;
|
jpayne@68
|
292 del++;
|
jpayne@68
|
293 }else if(state==MODE_INS){
|
jpayne@68
|
294 row--;
|
jpayne@68
|
295 boolean edge=(col<=1 || col>=columns);
|
jpayne@68
|
296 ins+=(edge ? 0 : 1);
|
jpayne@68
|
297 clip+=(edge ? 1 : 0);
|
jpayne@68
|
298 }else{
|
jpayne@68
|
299 assert(false) : state;
|
jpayne@68
|
300 }
|
jpayne@68
|
301 }
|
jpayne@68
|
302
|
jpayne@68
|
303 assert(row==0 || col==0);
|
jpayne@68
|
304 if(col!=row){//Not sure what this is doing
|
jpayne@68
|
305 while(row>0){
|
jpayne@68
|
306 clip++;
|
jpayne@68
|
307 row--;
|
jpayne@68
|
308 col--;
|
jpayne@68
|
309 }
|
jpayne@68
|
310 if(col>0){
|
jpayne@68
|
311 //do nothing
|
jpayne@68
|
312 }
|
jpayne@68
|
313 }
|
jpayne@68
|
314
|
jpayne@68
|
315 if(extra!=null){
|
jpayne@68
|
316 assert(extra.length==5);
|
jpayne@68
|
317 extra[0]=match;
|
jpayne@68
|
318 extra[1]=sub;
|
jpayne@68
|
319 extra[2]=del;
|
jpayne@68
|
320 extra[3]=ins;
|
jpayne@68
|
321 extra[4]=noref;
|
jpayne@68
|
322 extra[5]=clip;
|
jpayne@68
|
323 }
|
jpayne@68
|
324
|
jpayne@68
|
325 float len=match+sub+ins+del+noref*0.1f;
|
jpayne@68
|
326 float id=match/Tools.max(1.0f, len);
|
jpayne@68
|
327 return id;
|
jpayne@68
|
328 }
|
jpayne@68
|
329
|
jpayne@68
|
330 /** @return {score, bestRefStart, bestRefStop} */
|
jpayne@68
|
331 @Override
|
jpayne@68
|
332 public final int[] score(final byte[] read, final byte[] ref, final int refStartLoc, final int refEndLoc,
|
jpayne@68
|
333 final int maxRow, final int maxCol, final int maxState/*, final int maxScore, final int maxStart*/){
|
jpayne@68
|
334
|
jpayne@68
|
335 int row=maxRow;
|
jpayne@68
|
336 int col=maxCol;
|
jpayne@68
|
337 int state=maxState;
|
jpayne@68
|
338
|
jpayne@68
|
339 assert(maxState>=0 && maxState<packed.length) :
|
jpayne@68
|
340 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
|
jpayne@68
|
341 assert(maxRow>=0 && maxRow<packed.length) :
|
jpayne@68
|
342 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
|
jpayne@68
|
343 assert(maxCol>=0 && maxCol<packed[0].length) :
|
jpayne@68
|
344 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
|
jpayne@68
|
345
|
jpayne@68
|
346 int score=packed[maxRow][maxCol]&SCOREMASK; //Or zero, if it is to be recalculated
|
jpayne@68
|
347
|
jpayne@68
|
348 if(row<rows){
|
jpayne@68
|
349 int difR=rows-row;
|
jpayne@68
|
350 int difC=columns-col;
|
jpayne@68
|
351
|
jpayne@68
|
352 while(difR>difC){
|
jpayne@68
|
353 score+=POINTSoff_NOREF;
|
jpayne@68
|
354 difR--;
|
jpayne@68
|
355 }
|
jpayne@68
|
356
|
jpayne@68
|
357 row+=difR;
|
jpayne@68
|
358 col+=difR;
|
jpayne@68
|
359
|
jpayne@68
|
360 }
|
jpayne@68
|
361
|
jpayne@68
|
362 assert(refStartLoc<=refEndLoc);
|
jpayne@68
|
363 assert(row==rows);
|
jpayne@68
|
364
|
jpayne@68
|
365
|
jpayne@68
|
366 final int bestRefStop=refStartLoc+col-1;
|
jpayne@68
|
367
|
jpayne@68
|
368 while(row>0 && col>0){
|
jpayne@68
|
369 // final byte c=read[row-1];
|
jpayne@68
|
370 // final byte r=ref[refStartLoc+col-1];
|
jpayne@68
|
371 // final boolean defined=(AminoAcid.isFullyDefined(c) && AminoAcid.isFullyDefined(r));
|
jpayne@68
|
372 state=(packed[row][col]&MODEMASK);
|
jpayne@68
|
373 if(state==MODE_MATCH){
|
jpayne@68
|
374 col--;
|
jpayne@68
|
375 row--;
|
jpayne@68
|
376 }else if(state==MODE_SUB){
|
jpayne@68
|
377 col--;
|
jpayne@68
|
378 row--;
|
jpayne@68
|
379 }else if(state==MODE_N){
|
jpayne@68
|
380 col--;
|
jpayne@68
|
381 row--;
|
jpayne@68
|
382 }else if(state==MODE_DEL){
|
jpayne@68
|
383 col--;
|
jpayne@68
|
384 }else if(state==MODE_INS){
|
jpayne@68
|
385 row--;
|
jpayne@68
|
386 }else{
|
jpayne@68
|
387 assert(false) : state;
|
jpayne@68
|
388 }
|
jpayne@68
|
389 }
|
jpayne@68
|
390 // assert(false) : row+", "+col;
|
jpayne@68
|
391 if(row>col){
|
jpayne@68
|
392 col-=row;
|
jpayne@68
|
393 }
|
jpayne@68
|
394
|
jpayne@68
|
395 final int bestRefStart=refStartLoc+col;
|
jpayne@68
|
396
|
jpayne@68
|
397 score>>=SCOREOFFSET;
|
jpayne@68
|
398 // System.err.println("t2\t"+score+", "+maxScore+", "+maxStart+", "+bestRefStart);
|
jpayne@68
|
399 int[] rvec;
|
jpayne@68
|
400 if(bestRefStart<refStartLoc || bestRefStop>refEndLoc){ //Suggest extra padding in cases of overflow
|
jpayne@68
|
401 int padLeft=Tools.max(0, refStartLoc-bestRefStart);
|
jpayne@68
|
402 int padRight=Tools.max(0, bestRefStop-refEndLoc);
|
jpayne@68
|
403 rvec=new int[] {score, bestRefStart, bestRefStop, padLeft, padRight};
|
jpayne@68
|
404 }else{
|
jpayne@68
|
405 rvec=new int[] {score, bestRefStart, bestRefStop};
|
jpayne@68
|
406 }
|
jpayne@68
|
407 return rvec;
|
jpayne@68
|
408 }
|
jpayne@68
|
409
|
jpayne@68
|
410
|
jpayne@68
|
411 /** Will not fill areas that cannot match minScore.
|
jpayne@68
|
412 * @return {score, bestRefStart, bestRefStop} */
|
jpayne@68
|
413 @Override
|
jpayne@68
|
414 public final int[] fillAndScoreLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
|
jpayne@68
|
415 int a=Tools.max(0, refStartLoc);
|
jpayne@68
|
416 int b=Tools.min(ref.length-1, refEndLoc);
|
jpayne@68
|
417 assert(b>=a);
|
jpayne@68
|
418
|
jpayne@68
|
419 if(b-a>=maxColumns){
|
jpayne@68
|
420 System.err.println("Warning: Max alignment columns exceeded; restricting range. "+(b-a+1)+" > "+maxColumns);
|
jpayne@68
|
421 assert(false) : refStartLoc+", "+refEndLoc;
|
jpayne@68
|
422 b=Tools.min(ref.length-1, a+maxColumns-1);
|
jpayne@68
|
423 }
|
jpayne@68
|
424 int[] max=fillLimited(read, ref, a, b, minScore);
|
jpayne@68
|
425 // return max==null ? null : new int[] {max[3], 0, max[1]};
|
jpayne@68
|
426
|
jpayne@68
|
427 int[] score=(max==null ? null : score(read, ref, a, b, max[0], max[1], max[2]/*, max[3], max[4]*/));
|
jpayne@68
|
428
|
jpayne@68
|
429 return score;
|
jpayne@68
|
430 }
|
jpayne@68
|
431
|
jpayne@68
|
432 public static final String toString(byte[] ref, int startLoc, int stopLoc){
|
jpayne@68
|
433 StringBuilder sb=new StringBuilder(stopLoc-startLoc+1);
|
jpayne@68
|
434 for(int i=startLoc; i<=stopLoc; i++){sb.append((char)ref[i]);}
|
jpayne@68
|
435 return sb.toString();
|
jpayne@68
|
436 }
|
jpayne@68
|
437
|
jpayne@68
|
438 // public int maxScoreByIdentity(int len, float identity){
|
jpayne@68
|
439 // assert(identity>=0 && identity<=1);
|
jpayne@68
|
440 // return (int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
|
jpayne@68
|
441 // }
|
jpayne@68
|
442
|
jpayne@68
|
443 @Override
|
jpayne@68
|
444 public int minScoreByIdentity(int len, float identity){
|
jpayne@68
|
445 assert(identity>=0 && identity<=1);
|
jpayne@68
|
446
|
jpayne@68
|
447 int a=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
|
jpayne@68
|
448 int b=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_INS));
|
jpayne@68
|
449 int c=(int)(len*(1*POINTS_MATCH+((1/(Tools.max(identity, 0.000001f)))-1)*POINTS_DEL));
|
jpayne@68
|
450 return Tools.min(a, b, c);
|
jpayne@68
|
451 }
|
jpayne@68
|
452
|
jpayne@68
|
453 private static int calcDelScoreOffset(int len){
|
jpayne@68
|
454 if(len<=0){return 0;}
|
jpayne@68
|
455 int score=POINTSoff_DEL*len;
|
jpayne@68
|
456 return score;
|
jpayne@68
|
457 }
|
jpayne@68
|
458 //
|
jpayne@68
|
459 // public static int calcInsScore(int len){
|
jpayne@68
|
460 // if(len<=0){return 0;}
|
jpayne@68
|
461 // int score=POINTS_INS;
|
jpayne@68
|
462 //
|
jpayne@68
|
463 // if(len>1){
|
jpayne@68
|
464 // score+=(len-1)*POINTS_INS2;
|
jpayne@68
|
465 // }
|
jpayne@68
|
466 // return score;
|
jpayne@68
|
467 // }
|
jpayne@68
|
468 //
|
jpayne@68
|
469 // private static int calcInsScoreOffset(int len){
|
jpayne@68
|
470 // if(len<=0){return 0;}
|
jpayne@68
|
471 // int score=POINTSoff_INS;
|
jpayne@68
|
472 //
|
jpayne@68
|
473 // if(len>1){
|
jpayne@68
|
474 // score+=(len-1)*POINTSoff_INS2;
|
jpayne@68
|
475 // }
|
jpayne@68
|
476 // return score;
|
jpayne@68
|
477 // }
|
jpayne@68
|
478
|
jpayne@68
|
479 @Override
|
jpayne@68
|
480 public int rows(){return rows;}
|
jpayne@68
|
481 @Override
|
jpayne@68
|
482 public int columns(){return columns;}
|
jpayne@68
|
483
|
jpayne@68
|
484
|
jpayne@68
|
485 private int maxRows;
|
jpayne@68
|
486 private int maxColumns;
|
jpayne@68
|
487
|
jpayne@68
|
488 private int[][] packed;
|
jpayne@68
|
489
|
jpayne@68
|
490 public static final int MODEBITS=3;
|
jpayne@68
|
491 public static final int STARTBITS=9;
|
jpayne@68
|
492 public static final int LOWBITS=MODEBITS+STARTBITS;
|
jpayne@68
|
493 public static final int SCOREBITS=32-STARTBITS;
|
jpayne@68
|
494 public static final int MAX_START=((1<<STARTBITS)-1);
|
jpayne@68
|
495 public static final int MAX_SCORE=((1<<(SCOREBITS-1))-1)-2000;
|
jpayne@68
|
496 public static final int MIN_SCORE=0-MAX_SCORE; //Keeps it 1 point above "BAD".
|
jpayne@68
|
497
|
jpayne@68
|
498 public static final int STARTOFFSET=MODEBITS;
|
jpayne@68
|
499 public static final int SCOREOFFSET=LOWBITS;
|
jpayne@68
|
500
|
jpayne@68
|
501 public static final int MODEMASK=~((-1)<<MODEBITS);
|
jpayne@68
|
502 public static final int STARTMASK=(~((-1)<<STARTBITS))<<STARTOFFSET;
|
jpayne@68
|
503 public static final int SCOREMASK=(~((-1)<<SCOREBITS))<<SCOREOFFSET;
|
jpayne@68
|
504 public static final int HIGHMASK=SCOREMASK|STARTMASK;
|
jpayne@68
|
505
|
jpayne@68
|
506 //For some reason changing MODE_DEL from 1 to 0 breaks everything
|
jpayne@68
|
507 private static final byte MODE_DEL=1;
|
jpayne@68
|
508 private static final byte MODE_INS=2;
|
jpayne@68
|
509 private static final byte MODE_SUB=3;
|
jpayne@68
|
510 private static final byte MODE_MATCH=4;
|
jpayne@68
|
511 private static final byte MODE_N=5;
|
jpayne@68
|
512
|
jpayne@68
|
513 public static final int POINTS_NOREF=-15;
|
jpayne@68
|
514 public static final int POINTS_MATCH=100;
|
jpayne@68
|
515 public static final int POINTS_SUB=-50;
|
jpayne@68
|
516 public static final int POINTS_INS=-121;
|
jpayne@68
|
517 public static final int POINTS_DEL=-111;
|
jpayne@68
|
518
|
jpayne@68
|
519 public static final int BAD=MIN_SCORE-1;
|
jpayne@68
|
520
|
jpayne@68
|
521
|
jpayne@68
|
522 public static final int POINTSoff_NOREF=(POINTS_NOREF<<SCOREOFFSET);
|
jpayne@68
|
523 public static final int POINTSoff_NOREF_MODE_SUB=POINTSoff_NOREF|MODE_SUB;
|
jpayne@68
|
524 public static final int POINTSoff_MATCH=(POINTS_MATCH<<SCOREOFFSET);
|
jpayne@68
|
525 public static final int POINTSoff_SUB=(POINTS_SUB<<SCOREOFFSET);
|
jpayne@68
|
526 public static final int POINTSoff_INS=(POINTS_INS<<SCOREOFFSET);
|
jpayne@68
|
527 public static final int POINTSoff_DEL=(POINTS_DEL<<SCOREOFFSET);
|
jpayne@68
|
528 public static final int BADoff=(BAD<<SCOREOFFSET);
|
jpayne@68
|
529
|
jpayne@68
|
530 private int rows;
|
jpayne@68
|
531 private int columns;
|
jpayne@68
|
532
|
jpayne@68
|
533 // public long iterationsLimited=0;
|
jpayne@68
|
534 // public long iterationsUnlimited=0;
|
jpayne@68
|
535
|
jpayne@68
|
536 public boolean verbose=false;
|
jpayne@68
|
537 public boolean verbose2=false;
|
jpayne@68
|
538
|
jpayne@68
|
539 }
|