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 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 if(loc==-1){//substitution
if(lastValue<0 && timeInMode>0){//contiguous
if(timeInMode1) : minContig;
int contig=0;
int maxContig=0;
int score=0;
int lastLoc=-3; //Last true location
int lastValue=-1;
int timeInMode=0;
for(int i=0; i0){//match
if(loc==lastValue){//contiguous match
contig++;
score+=(POINTS_MATCH2+baseScores[i]);
}else if(loc==lastLoc || lastLoc<0){//match after a sub, or first match
maxContig=Tools.max(maxContig, contig);
contig=1;
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
maxContig=Tools.max(maxContig, contig);
contig=0;
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 if(loc==-1){//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; i=ref.length) ? (byte)'N' : ref[j];
if(c=='N' || r=='N'){match[i]='N';}
else if(c==r){match[i]='m';}
else{match[i]='S';}
}
return match;
}
@Override
public final int scoreNoIndels(byte[] read, byte[] ref, byte[] baseScores, final int refStart){
int score=0;
int mode=-1;
int timeInMode=0;
//This block handles cases where the read runs outside the reference
//Of course, padding the reference with 'N' would be better, but...
int readStart=0;
int readStop=read.length;
final int refStop=refStart+read.length;
if(refStart<0){
readStart=0-refStart;
score+=POINTS_NOREF*readStart;
}
if(refStop>ref.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){return -99999;}
if(refStart<0){
readStart=0-refStart;
score+=POINTS_NOREF*readStart;
}
if(refStop>ref.length){
int dif=(refStop-ref.length);
System.err.println(new String(read)+"\ndif="+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){return -99999;}
if(refStart<0){
readStart=0-refStart;
score+=POINTS_NOREF*readStart;
}
if(refStop>ref.length){
int dif=(refStop-ref.length);
System.err.println(new String(read)+"\ndif="+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; iMINGAP){
int rem=len%GAPLEN;
int div=(len-GAPBUFFER2)/GAPLEN;
score+=(div*POINTS_GAP);
assert(rem+GAPBUFFER2LIMIT_FOR_COST_4); //and probably LIMIT_FOR_COST_5
// assert(false) : div;
}
if(len>LIMIT_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;
}
@Override
public 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;
}
private final int[][][] packed;
/** Banded packed matrix */
private final int[][][] bpacked;
/** Start locations of banded matrix */
private final int[][][] startmatrix;
/** Score for column zero */
private final int[] col0score;
private final byte[] grefbuffer;
private int greflimit=-1;
private int greflimit2=-1;
private int grefRefOrigin=-1;
/**DO NOT MODIFY*/
public final byte[] getGrefbuffer(){
return grefbuffer;
}
public final int[] vertLimit;
public final int[] horizLimit;
public CharSequence showVertLimit(){
StringBuilder sb=new StringBuilder();
for(int i=0; i<=rows; i++){sb.append(vertLimit[i]>>SCOREOFFSET).append(",");}
return sb;
}
public CharSequence showHorizLimit(){
StringBuilder sb=new StringBuilder();
for(int i=0; i<=columns; i++){sb.append(horizLimit[i]>>SCOREOFFSET).append(",");}
return sb;
}
public static float minIdToMinRatio(double minid){
if(minid>1){minid=minid/100;}
assert(minid>0 && minid<=1) : "Min identity should be between 0 and 1. Values above 1 will be assumed to be percent and divided by 100.";
double matchdif=POINTS_MATCH-POINTS_MATCH2;
double match=POINTS_MATCH2;
double sub=-POINTS_MATCH2+0.5*(matchdif+POINTS_SUB)+0.5*POINTS_SUB2;
double del=0.1*(matchdif+POINTS_DEL)+0.2*POINTS_DEL2+0.4*POINTS_DEL3+0.3*POINTS_DEL4;
double ins=-POINTS_MATCH2+0.4*(matchdif+POINTS_INS)+0.3*(POINTS_INS2)+0.3*(POINTS_INS3);
double badAvg=.7*sub+.2*del+.1*ins;
double badFraction=1-minid;
double minratio=(match+badFraction*badAvg)/match;
assert(minratio<=1);
minratio=Tools.max(0.05, minratio);
return (float)minratio;
}
public static final int TIMEBITS=11;
public static final int SCOREBITS=32-TIMEBITS;
public static final int MAX_TIME=((1<