0);
for(int r=y+1; r<=lim; r++, gpos++){
gref[gpos]=ref[r];
}
for(int g=0; ggref.length){
System.err.println("gref buffer overflow: "+lim+" > "+gref.length);
return null;
}
for(int i=greflimit, r=refEndLoc+1; i "+j);}
return j;
}
j+=(c==GAPC ? GAPLEN : 1);
// if(c!=GAPC){j++;}
// else{j+=GAPLEN;}
}
System.err.println(grefRefOrigin);
System.err.println(point);
System.err.println(new String(gref));
throw new RuntimeException("Out of bounds.");
}
private final int translateToGappedCoordinate(int point, byte[] gref){
if(verbose){System.err.println("translateToGappedCoordinate("+point+"), gro="+grefRefOrigin+", grl="+greflimit);}
if(point<=grefRefOrigin){return point-grefRefOrigin;}
for(int i=0, j=grefRefOrigin; i "+i);}
return i;
}
j+=(c==GAPC ? GAPLEN : 1);
// if(c!=GAPC){j++;}
// else{j+=GAPLEN;}
}
System.err.println(grefRefOrigin);
System.err.println(point);
System.err.println(new String(gref));
throw new RuntimeException("Out of bounds.");
}
public final int[] fillAndScoreLimited(byte[] read, int chrom, int start, int stop, int thresh, int minScore, int[] gaps){
return fillAndScoreLimited(read, Data.getChromosome(chrom).array, start-thresh, stop+thresh, minScore, gaps);
}
@Deprecated
public final int[] fillAndScoreQ(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, byte[] baseScores){
int a=Tools.max(0, refStartLoc);
int b=Tools.min(ref.length-1, refEndLoc);
assert(b>=a);
if(b-a>=maxColumns){
System.err.println("Warning: Max alignment columns exceeded; restricting range. "+(b-a+1)+" > "+maxColumns);
b=Tools.min(ref.length-1, a+maxColumns-1);
}
int[] max=fillQ(read, ref, baseScores, a, b);
// int[] score=score(read, ref, a, b, max[0], max[1], max[2]);
// return score;
return null;
}
@Deprecated
public final int[] fillAndScoreQ(byte[] read, SiteScore ss, int thresh, byte[] baseScores){
return fillAndScoreQ(read, ss.chrom, ss.start, ss.stop, thresh, baseScores);
}
@Deprecated
public final int[] fillAndScoreQ(byte[] read, int chrom, int start, int stop, int thresh, byte[] baseScores){
return fillAndScoreQ(read, Data.getChromosome(chrom).array, start-thresh, stop+thresh, baseScores);
}
// public final int scoreNoIndels(byte[] read, SiteScore ss){
//
// ChromosomeArray cha=Data.getChromosome(ss.chrom);
// final int refStart=ss.start;
//
// int score=0;
// int mode=MODE_START;
// int timeInMode=0;
// if(refStart<0 || refStart+read.length>cha.maxIndex+1){return -99999;} //TODO: Partial match
//
// for(int i=0; icha.maxIndex+1){
// int dif=(cha.maxIndex+1-refStop);
// readStop-=dif;
// score+=POINTSoff_NOREF*dif;
// }
//
//// if(refStart<0 || refStart+read.length>cha.maxIndex+1){return -99999;} //No longer needed.
//
// for(int i=readStart; i0){//match
if(loc==lastValue){//contiguous match
score+=POINTS_MATCH2;
}else if(loc==lastLoc || lastLoc<0){//match after a sub, or first match
score+=POINTS_MATCH;
}else if(loc=0);
score+=POINTS_MATCH;
score+=POINTS_DEL;
int dif=lastLoc-loc+1;
if(dif>MINGAP){
int rem=dif%GAPLEN;
int div=(dif-GAPBUFFER2)/GAPLEN;
score+=(div*POINTS_GAP);
assert(rem+GAPBUFFER2LIMIT_FOR_COST_4); //and probably LIMIT_FOR_COST_5
// assert(false) : div;
}
if(dif>LIMIT_FOR_COST_5){
score+=((dif-LIMIT_FOR_COST_5+MASK5)/TIMESLIP)*POINTS_DEL5;
dif=LIMIT_FOR_COST_5;
}
if(dif>LIMIT_FOR_COST_4){
score+=(dif-LIMIT_FOR_COST_4)*POINTS_DEL4;
dif=LIMIT_FOR_COST_4;
}
if(dif>LIMIT_FOR_COST_3){
score+=(dif-LIMIT_FOR_COST_3)*POINTS_DEL3;
dif=LIMIT_FOR_COST_3;
}
if(dif>1){
score+=(dif-1)*POINTS_DEL2;
}
timeInMode=1;
}else if(loc>lastLoc){//insertion
assert(lastLoc>=0);
score+=POINTS_MATCH;
score+=POINTS_INS;
int dif=Tools.min(loc-lastLoc+1, 5);
assert(dif>0);
if(dif>LIMIT_FOR_COST_4){
score+=(dif-LIMIT_FOR_COST_4)*POINTS_INS4;
dif=LIMIT_FOR_COST_4;
}
if(dif>LIMIT_FOR_COST_3){
score+=(dif-LIMIT_FOR_COST_3)*POINTS_INS3;
dif=LIMIT_FOR_COST_3;
}
if(dif>1){
score+=(dif-1)*POINTS_INS2;
}
timeInMode=1;
}else{
assert(false);
}
lastLoc=loc;
}else{//substitution
if(lastValue<0 && timeInMode>0){//contiguous
if(timeInMode0){//match
if(loc==lastValue){//contiguous match
score+=(POINTS_MATCH2+baseScores[i]);
}else if(loc==lastLoc || lastLoc<0){//match after a sub, or first match
score+=(POINTS_MATCH+baseScores[i]);
}else if(loc=0);
score+=(POINTS_MATCH+baseScores[i]);
score+=POINTS_DEL;
int dif=lastLoc-loc+1;
if(dif>MINGAP){
int rem=dif%GAPLEN;
int div=(dif-GAPBUFFER2)/GAPLEN;
score+=(div*POINTS_GAP);
assert(rem+GAPBUFFER2LIMIT_FOR_COST_4); //and probably LIMIT_FOR_COST_5
// assert(false) : div;
}
if(dif>LIMIT_FOR_COST_5){
score+=((dif-LIMIT_FOR_COST_5+MASK5)/TIMESLIP)*POINTS_DEL5;
dif=LIMIT_FOR_COST_5;
}
if(dif>LIMIT_FOR_COST_4){
score+=(dif-LIMIT_FOR_COST_4)*POINTS_DEL4;
dif=LIMIT_FOR_COST_4;
}
if(dif>LIMIT_FOR_COST_3){
score+=(dif-LIMIT_FOR_COST_3)*POINTS_DEL3;
dif=LIMIT_FOR_COST_3;
}
if(dif>1){
score+=(dif-1)*POINTS_DEL2;
}
timeInMode=1;
}else if(loc>lastLoc){//insertion
assert(lastLoc>=0);
score+=(POINTS_MATCH+baseScores[i]);
score+=POINTS_INS;
int dif=Tools.min(loc-lastLoc+1, 5);
assert(dif>0);
if(dif>LIMIT_FOR_COST_4){
score+=(dif-LIMIT_FOR_COST_4)*POINTS_INS4;
dif=LIMIT_FOR_COST_4;
}
if(dif>LIMIT_FOR_COST_3){
score+=(dif-LIMIT_FOR_COST_3)*POINTS_INS3;
dif=LIMIT_FOR_COST_3;
}
if(dif>1){
score+=(dif-1)*POINTS_INS2;
}
timeInMode=1;
}else{
assert(false);
}
lastLoc=loc;
}else{//substitution
if(lastValue<0 && timeInMode>0){//contiguous
if(timeInModeref.length){
int dif=(refStop-ref.length);
readStop-=dif;
score+=POINTS_NOREF*dif;
}
// if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
for(int i=readStart; iref.length){
int dif=(refStop-ref.length);
readStop-=dif;
score+=POINTS_NOREF*dif;
}
// if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
for(int i=readStart; iref.length){
int dif=(refStop-ref.length);
System.err.println("dif="+dif+", ref.length="+ref.length+", refStop="+refStop);
readStop-=dif;
score+=POINTS_NOREF*dif;
}
assert(refStart+readStop<=ref.length) : "readStart="+readStart+", readStop="+readStop+
", refStart="+refStart+", refStop="+refStop+", ref.length="+ref.length+", read.length="+read.length;
assert(matchReturn!=null);
assert(matchReturn.length==1);
if(matchReturn[0]==null || matchReturn[0].length!=read.length){
assert(matchReturn[0]==null || matchReturn[0].lengthref.length){return -99999;} //No longer needed.
for(int i=readStart; iref.length){
int dif=(refStop-ref.length);
System.err.println("dif="+dif+", ref.length="+ref.length+", refStop="+refStop);
readStop-=dif;
score+=POINTS_NOREF*dif;
}
assert(refStart+readStop<=ref.length) : "readStart="+readStart+", readStop="+readStop+
", refStart="+refStart+", refStop="+refStop+", ref.length="+ref.length+", read.length="+read.length;
assert(matchReturn!=null);
assert(matchReturn.length==1);
if(matchReturn[0]==null || matchReturn[0].length!=read.length){
assert(matchReturn[0]==null || matchReturn[0].lengthref.length){return -99999;} //No longer needed.
for(int i=readStart; i=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
for(int i=0; i=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
for(int i=0; i>SCOREOFFSET;
String s=" "+num;
if(s.length()>width){s=num>0 ? maxString : minString;}
int spaces=width-s.length();
assert(spaces>=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
for(int i=0; i=0);
for(int i=0; iLIMIT_FOR_COST_5){
score+=((len-LIMIT_FOR_COST_5+MASK5)/TIMESLIP)*POINTS_DEL5;
len=LIMIT_FOR_COST_5;
}
if(len>LIMIT_FOR_COST_4){
score+=(len-LIMIT_FOR_COST_4)*POINTS_DEL4;
len=LIMIT_FOR_COST_4;
}
if(len>LIMIT_FOR_COST_3){
score+=(len-LIMIT_FOR_COST_3)*POINTS_DEL3;
len=LIMIT_FOR_COST_3;
}
if(len>1){
score+=(len-1)*POINTS_DEL2;
}
return score;
}
private static int calcDelScoreOffset(int len){
if(len<=0){return 0;}
int score=POINTSoff_DEL;
if(len>LIMIT_FOR_COST_5){
score+=((len-LIMIT_FOR_COST_5+MASK5)/TIMESLIP)*POINTSoff_DEL5;
len=LIMIT_FOR_COST_5;
}
if(len>LIMIT_FOR_COST_4){
score+=(len-LIMIT_FOR_COST_4)*POINTSoff_DEL4;
len=LIMIT_FOR_COST_4;
}
if(len>LIMIT_FOR_COST_3){
score+=(len-LIMIT_FOR_COST_3)*POINTSoff_DEL3;
len=LIMIT_FOR_COST_3;
}
if(len>1){
score+=(len-1)*POINTSoff_DEL2;
}
return score;
}
public static int calcInsScore(int len){
if(len<=0){return 0;}
int score=POINTS_INS;
if(len>LIMIT_FOR_COST_4){
score+=(len-LIMIT_FOR_COST_4)*POINTS_INS4;
len=LIMIT_FOR_COST_4;
}
if(len>LIMIT_FOR_COST_3){
score+=(len-LIMIT_FOR_COST_3)*POINTS_INS3;
len=LIMIT_FOR_COST_3;
}
if(len>1){
score+=(len-1)*POINTS_INS2;
}
return score;
}
private static int calcInsScoreOffset(int len){
if(len<=0){return 0;}
int score=POINTSoff_INS;
if(len>LIMIT_FOR_COST_4){
score+=(len-LIMIT_FOR_COST_4)*POINTSoff_INS4;
len=LIMIT_FOR_COST_4;
}
if(len>LIMIT_FOR_COST_3){
score+=(len-LIMIT_FOR_COST_3)*POINTSoff_INS3;
len=LIMIT_FOR_COST_3;
}
if(len>1){
score+=(len-1)*POINTSoff_INS2;
}
return score;
}
public final int maxRows;
public final int maxColumns;
private final int[][][] packed;
private final byte[] grefbuffer;
private int greflimit=-1;
private int greflimit2=-1;
private int grefRefOrigin=-1;
public static final int GAPBUFFER=Shared.GAPBUFFER;
public static final int GAPBUFFER2=Shared.GAPBUFFER2;
public static final int GAPLEN=Shared.GAPLEN;
public static final int MINGAP=Shared.MINGAP;
public static final int GAPCOST=Shared.GAPCOST;
public static final byte GAPC=Shared.GAPC;
private static final int GREFLIMIT2_CUSHION=128; //Tools.max(GAPBUFFER2, GAPLEN);
/**DO NOT MODIFY*/
public final byte[] getGrefbuffer(){
return grefbuffer;
}
public final int[] vertLimit;
public final int[] horizLimit;
CharSequence showVertLimit(){
StringBuilder sb=new StringBuilder();
for(int i=0; i<=rows; i++){sb.append(vertLimit[i]>>SCOREOFFSET).append(",");}
return sb;
}
CharSequence showHorizLimit(){
StringBuilder sb=new StringBuilder();
for(int i=0; i<=columns; i++){sb.append(horizLimit[i]>>SCOREOFFSET).append(",");}
return sb;
}
// public static final int MODEBITS=2;
public static final int TIMEBITS=12;
public static final int SCOREBITS=32-TIMEBITS;
public static final int MAX_TIME=((1<