Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/sketch/Comparison.java @ 68:5028fdace37b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
rev   line source
jpayne@68 1 package sketch;
jpayne@68 2
jpayne@68 3 import java.util.Comparator;
jpayne@68 4 import java.util.Locale;
jpayne@68 5
jpayne@68 6 import aligner.Aligner;
jpayne@68 7 import prok.GeneCaller;
jpayne@68 8 import shared.Tools;
jpayne@68 9 import tax.TaxNode;
jpayne@68 10
jpayne@68 11 public final class Comparison extends SketchObject implements Comparable<Comparison> {
jpayne@68 12
jpayne@68 13 /*--------------------------------------------------------------*/
jpayne@68 14 /*---------------- Constructors ----------------*/
jpayne@68 15 /*--------------------------------------------------------------*/
jpayne@68 16
jpayne@68 17 // public Comparison(CompareBuffer buffer){
jpayne@68 18 // this(buffer, null, null);
jpayne@68 19 // }
jpayne@68 20
jpayne@68 21 // public Comparison(Sketch a_, Sketch b_){
jpayne@68 22 // this(null, a_, b_);
jpayne@68 23 // }
jpayne@68 24
jpayne@68 25 public Comparison(CompareBuffer buffer, Sketch a_, Sketch b_){
jpayne@68 26
jpayne@68 27 a=a_;
jpayne@68 28 b=b_;
jpayne@68 29
jpayne@68 30 if(buffer!=null){setFrom(buffer);}
jpayne@68 31
jpayne@68 32 if(b!=null){
jpayne@68 33 taxName=b.taxName();
jpayne@68 34 taxID=b.taxID;
jpayne@68 35 }
jpayne@68 36
jpayne@68 37 // System.err.println(this);
jpayne@68 38 // System.err.println(b.present);
jpayne@68 39 }
jpayne@68 40
jpayne@68 41 /*--------------------------------------------------------------*/
jpayne@68 42 /*---------------- Mutators ----------------*/
jpayne@68 43 /*--------------------------------------------------------------*/
jpayne@68 44
jpayne@68 45 public void setFrom(CompareBuffer buffer){
jpayne@68 46 hits=buffer.hits();
jpayne@68 47 multiHits=buffer.multiHits();
jpayne@68 48 unique2=buffer.unique2();
jpayne@68 49 unique3=buffer.unique3();
jpayne@68 50 noHits=buffer.noHits();
jpayne@68 51
jpayne@68 52 contamHits=buffer.contamHits();
jpayne@68 53 contam2Hits=buffer.contam2Hits();
jpayne@68 54 multiContamHits=buffer.multiContamHits();
jpayne@68 55
jpayne@68 56 refDivisor=buffer.refDivisor();
jpayne@68 57 queryDivisor=buffer.queryDivisor();
jpayne@68 58
jpayne@68 59 refSize=buffer.refSize();
jpayne@68 60 querySize=buffer.querySize();
jpayne@68 61
jpayne@68 62 depth=buffer.depth();
jpayne@68 63 depth2=buffer.depth2();
jpayne@68 64 float x=buffer.avgRefHits();
jpayne@68 65 if(x>0){avgRefHits=x;}
jpayne@68 66 // volume=volume0();
jpayne@68 67 score=score0();
jpayne@68 68
jpayne@68 69 hits1=buffer.hits1();
jpayne@68 70 qSeen1=buffer.qSeen1();
jpayne@68 71 rSeen1=buffer.rSeen1();
jpayne@68 72 }
jpayne@68 73
jpayne@68 74 public void recompare(CompareBuffer buffer, int[][] taxHits, int contamLevel){
jpayne@68 75
jpayne@68 76 // for(int[] row : taxHits){
jpayne@68 77 // if(row!=null){
jpayne@68 78 // System.err.println(Arrays.toString(row));
jpayne@68 79 // }
jpayne@68 80 // }//Tested; correctly indicates most rows have octopus but some have other things.
jpayne@68 81
jpayne@68 82 assert(a.merged());
jpayne@68 83 // int oldContam2=contam2Hits;
jpayne@68 84 int x=a.countMatches(b, buffer, a.compareBitSet(), false, taxHits, contamLevel);
jpayne@68 85 assert(x==hits);
jpayne@68 86 setFrom(buffer);
jpayne@68 87 // contam2Hits=oldContam2;
jpayne@68 88 }
jpayne@68 89
jpayne@68 90 /*--------------------------------------------------------------*/
jpayne@68 91 /*---------------- Methods ----------------*/
jpayne@68 92 /*--------------------------------------------------------------*/
jpayne@68 93
jpayne@68 94 @Override
jpayne@68 95 public boolean equals(Object b){
jpayne@68 96 if(b==null || b.getClass()!=this.getClass()){return false;}
jpayne@68 97 return scoreComparator.compare(this, (Comparison)b)==0;
jpayne@68 98 }
jpayne@68 99
jpayne@68 100 // //WKID
jpayne@68 101 // public float wkid(){return idMinDivisor();}
jpayne@68 102 // public float idMinDivisor(){
jpayne@68 103 // return hits/(float)minDivisor();
jpayne@68 104 // }
jpayne@68 105
jpayne@68 106 // public float k1Fraction(){
jpayne@68 107 // return a.k1Fraction();
jpayne@68 108 // }
jpayne@68 109
jpayne@68 110 // //KID
jpayne@68 111 // public float kid(){return idMaxDivisor();}
jpayne@68 112 // public float idMaxDivisor(){
jpayne@68 113 // return hits/(float)maxDivisor();
jpayne@68 114 // }
jpayne@68 115
jpayne@68 116 public float idQueryDivisor(){
jpayne@68 117 return hits/(float)(Tools.max(1, refDivisor));
jpayne@68 118 }
jpayne@68 119
jpayne@68 120 public float idRefDivisor(){
jpayne@68 121 return hits/(float)(Tools.max(1, refDivisor));
jpayne@68 122 }
jpayne@68 123
jpayne@68 124 public float completeness(){
jpayne@68 125 float complt=Tools.min(1, (queryDivisor-contamHits)/(float)Tools.max(1, refDivisor));
jpayne@68 126 return complt;
jpayne@68 127 // float c2=hits/(float)refDivisor;
jpayne@68 128 // assert(queryDivisor-contamHits>=hits);
jpayne@68 129 // assert(c1>=c2);
jpayne@68 130 // System.err.println(hits+", "+contamHits+", "+refDivisor+", "+queryDivisor+", "+c1+", "+c2);
jpayne@68 131 // return Tools.max(c1, c2);
jpayne@68 132 // float kid=idMaxDivisor(), wkid=idMinDivisor();
jpayne@68 133 // return kid==0 ? 0 : kid/wkid;
jpayne@68 134 }
jpayne@68 135
jpayne@68 136 public float contamFraction(){
jpayne@68 137 return Tools.min(1, contamHits/(float)Tools.max(1, queryDivisor));
jpayne@68 138 }
jpayne@68 139
jpayne@68 140 public float contam2Fraction(){
jpayne@68 141 return Tools.min(1, contam2Hits/(float)Tools.max(1, queryDivisor));
jpayne@68 142 }
jpayne@68 143
jpayne@68 144 public float uContamFraction() {
jpayne@68 145 int uContamHits=contamHits-multiContamHits;
jpayne@68 146 return Tools.min(1, uContamHits/(float)Tools.max(1, queryDivisor));
jpayne@68 147 }
jpayne@68 148
jpayne@68 149
jpayne@68 150
jpayne@68 151 /*--------------------------------------------------------------*/
jpayne@68 152
jpayne@68 153 final float wkid(){
jpayne@68 154 final int div=minDivisor();
jpayne@68 155 return hits/(float)div;
jpayne@68 156 }
jpayne@68 157 final float kid(){
jpayne@68 158 final int div=maxDivisor();
jpayne@68 159 return hits/(float)div;
jpayne@68 160 }
jpayne@68 161 final float aniOld(){
jpayne@68 162 float wkid=wkid();
jpayne@68 163 // final float ani=wkidToAni(wkid, k1Fraction());
jpayne@68 164 final float ani=wkidToAni(wkid);
jpayne@68 165 return ani;
jpayne@68 166 }
jpayne@68 167 final float ani(){
jpayne@68 168 if(hits<1){return 0;}
jpayne@68 169 final float ani;
jpayne@68 170 if(k2>0 && useToValue2){
jpayne@68 171 float ani1=ani1();
jpayne@68 172 float ani2=ani2();
jpayne@68 173 // ani=0.5f*(ani1+ani2);
jpayne@68 174 ani=0.5f*(Tools.max(0.9f*ani2, ani1)+Tools.max(0.8f*ani1, ani2));
jpayne@68 175 // return (ani1*qSeen1+ani2*qSeen2())/queryDivisor;
jpayne@68 176
jpayne@68 177 // System.err.println("ani="+ani+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)aniDual()+"\n"
jpayne@68 178 //// +"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
jpayne@68 179 // + "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
jpayne@68 180 // + "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
jpayne@68 181 }else{
jpayne@68 182 ani=aniOld();
jpayne@68 183 }
jpayne@68 184 return ani;
jpayne@68 185 }
jpayne@68 186
jpayne@68 187 final float wkid1(){
jpayne@68 188 final int div=minDivisor1();
jpayne@68 189 return hits1()/(float)div;
jpayne@68 190 }
jpayne@68 191 final float kid1(){
jpayne@68 192 final int div=maxDivisor1();
jpayne@68 193 return hits1()/(float)div;
jpayne@68 194 }
jpayne@68 195 final float ani1(){
jpayne@68 196 float wkid=wkid1();
jpayne@68 197 final float ani=wkidToAniExact(wkid, k);
jpayne@68 198 return ani;
jpayne@68 199 }
jpayne@68 200
jpayne@68 201 final float wkid2(){
jpayne@68 202 final int div=minDivisor2();
jpayne@68 203 return hits2()/(float)div;
jpayne@68 204 }
jpayne@68 205 final float kid2(){
jpayne@68 206 final int div=maxDivisor2();
jpayne@68 207 return hits2()/(float)div;
jpayne@68 208 }
jpayne@68 209 final float ani2(){
jpayne@68 210 assert(k2>0);
jpayne@68 211 float wkid=wkid2();
jpayne@68 212 final float ani=wkidToAniExact(wkid, k2);
jpayne@68 213 return ani;
jpayne@68 214 }
jpayne@68 215
jpayne@68 216 final float aniDual(){
jpayne@68 217 assert(k2>0);
jpayne@68 218 float wkid1=wkid1();
jpayne@68 219 float wkid2=wkid2();
jpayne@68 220 float ratio=(wkid1/wkid2);
jpayne@68 221 float exp=1f/(k-k2);//TODO - make this initialized
jpayne@68 222 double ani=Math.pow(ratio, exp);
jpayne@68 223 double gf=wkid2/Math.pow(ani, k2);
jpayne@68 224
jpayne@68 225 // System.err.println("ani="+ani()+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)ani+"\n"
jpayne@68 226 // +"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
jpayne@68 227 // + "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
jpayne@68 228 // + "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
jpayne@68 229
jpayne@68 230 return (float)ani;
jpayne@68 231 }
jpayne@68 232
jpayne@68 233 /*--------------------------------------------------------------*/
jpayne@68 234
jpayne@68 235 int hits1(){return hits1;}
jpayne@68 236 int qSeen1(){return qSeen1;}
jpayne@68 237 int rSeen1(){return rSeen1;}
jpayne@68 238 int minDivisor1(){return Tools.max(1, Tools.min(qSeen1, rSeen1));}
jpayne@68 239 int maxDivisor1(){return Tools.max(1, qSeen1, rSeen1);}
jpayne@68 240
jpayne@68 241 int hits2(){return hits-hits1;}
jpayne@68 242 int qSeen2(){return queryDivisor-qSeen1;}
jpayne@68 243 int rSeen2(){return refDivisor-rSeen1;}
jpayne@68 244 int minDivisor2(){return Tools.max(1, Tools.min(qSeen2(), rSeen2()));}
jpayne@68 245 int maxDivisor2(){return Tools.max(1, qSeen2(), rSeen2());}
jpayne@68 246
jpayne@68 247 /*--------------------------------------------------------------*/
jpayne@68 248
jpayne@68 249 // public float aniOld(){
jpayne@68 250 // if(hits<1){return 0;}
jpayne@68 251 //
jpayne@68 252 //// double wkid=aniFromWKID ? idMinDivisor() : idMaxDivisor();
jpayne@68 253 // double wkid=idMinDivisor();
jpayne@68 254 // return wkidToAni(wkid, k1Fraction());
jpayne@68 255 //
jpayne@68 256 //// final float rID=hits/(float)(refDivisor);
jpayne@68 257 //// final float qID=hits/(float)(queryDivisor-contamHits);
jpayne@68 258 //// final float wkid2=Tools.max(qID, rID);
jpayne@68 259 //// final float ani=wkidToAni(wkid2);
jpayne@68 260 ////
jpayne@68 261 ////// System.err.println("rid: "+wkidToAni(rID)+", qid: "+wkidToAni(qID)+", qid2: "+wkidToAni(hits/(float)(queryDivisor)));
jpayne@68 262 ////
jpayne@68 263 //// return ani;
jpayne@68 264 // }
jpayne@68 265
jpayne@68 266 int minDivisor(){return Tools.max(1, Tools.min(refDivisor, queryDivisor));}
jpayne@68 267 int maxDivisor(){return Tools.max(1, refDivisor, queryDivisor);}
jpayne@68 268
jpayne@68 269 private float score0_old(){
jpayne@68 270 long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
jpayne@68 271 float wkid=wkid();
jpayne@68 272 float kid=kid();
jpayne@68 273 float complt=completeness();
jpayne@68 274 float contam=contamFraction();
jpayne@68 275 float refHits=Tools.max(avgRefHits, 1f);
jpayne@68 276 float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
jpayne@68 277 return (float)((Math.log(hits+2)*.25f+0.5f)*(refHitMult*0.2*Math.log(hits+2+(1.2*uHits()+0.25*unique2()+0.1*unique3()))
jpayne@68 278 *Math.sqrt(40*(20000+hits+uHits())*(wkid*kid*Math.pow(est*complt, 0.2)*(1-contam*0.1)))+0.1));
jpayne@68 279 }
jpayne@68 280
jpayne@68 281 private float score0(){
jpayne@68 282 final long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
jpayne@68 283 final float wkid=wkid();
jpayne@68 284 final float kid=kid();
jpayne@68 285 final float ani=ani();
jpayne@68 286 final float complt=completeness();
jpayne@68 287 final float contam=contamFraction();
jpayne@68 288 final float refHits=Tools.max(avgRefHits, 1f);
jpayne@68 289 final float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
jpayne@68 290 final float uhits=uHits();
jpayne@68 291 final float uhits2=unique2();
jpayne@68 292 final float uhits3=unique3();
jpayne@68 293 final float contamMult=(1-contam*0.95f);
jpayne@68 294 final float estMult=(float)(Math.pow(est, 0.2)*Math.sqrt(complt));
jpayne@68 295 final float aniMult=(float)(ani*Math.sqrt(wkid*kid));
jpayne@68 296
jpayne@68 297 final float hitsSum=1+hits+uhits+0.5f*uhits2+0.25f*uhits3;
jpayne@68 298
jpayne@68 299 final float score=(float)(Math.log(Tools.max(1.2f, hits-1))*hitsSum*refHitMult*contamMult*aniMult*estMult);//+(Math.log(hitsSum+1)*wkid*complt));
jpayne@68 300 return (float)(8*Math.sqrt(score));
jpayne@68 301 }
jpayne@68 302
jpayne@68 303 private long range(){//TODO Make sure these are calculated correctly; it seems like one divisor might be 1 higher than necessary.
jpayne@68 304 long maxA=a.keys[Tools.max(0, queryDivisor-1)];
jpayne@68 305 long maxB=b.keys[Tools.max(0, refDivisor-1)];
jpayne@68 306 // assert(false) : Tools.max(0, queryDivisor-1)+", "+Tools.max(0, refDivisor-1)+
jpayne@68 307 // ", "+a.array[Tools.max(0, queryDivisor-1)]+", "+b.array[Tools.max(0, refDivisor-1)]+", "+Tools.max(maxA, maxB);//+"\n\n"+Arrays.toString(a.array)+"\n\n"+Arrays.toString(b.array);
jpayne@68 308 return Tools.min(maxA, maxB);
jpayne@68 309 }
jpayne@68 310
jpayne@68 311 private static double eValue(int hits, int minDiv, int maxDiv, long range){
jpayne@68 312 if(hits>=range || maxDiv>=range){return 1.0;}
jpayne@68 313 double probHit=maxDiv/(double)range;//Saturation of range
jpayne@68 314 // double probNoHit=1-probHit;
jpayne@68 315 double eValue=Math.pow(probHit, hits); //This is a simplification, assuming hits are very improbable.
jpayne@68 316 //Note that this does not take into account minDiv, the number of attempts... but it should.
jpayne@68 317 // System.err.println("hits="+hits+", minDiv="+minDiv+", maxDiv="+maxDiv+", range="+range+", eValue="+eValue);
jpayne@68 318 return eValue;
jpayne@68 319 }
jpayne@68 320
jpayne@68 321 public double eValue(){
jpayne@68 322 double eValue=eValue1()*eValue2();
jpayne@68 323 // System.err.println("eValue="+eValue);
jpayne@68 324 return eValue;
jpayne@68 325 }
jpayne@68 326
jpayne@68 327 public double eValue1(){
jpayne@68 328 long range0=range();
jpayne@68 329 int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k;
jpayne@68 330 double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
jpayne@68 331 int hits=hits1;
jpayne@68 332 int minDiv=Tools.min(qSeen1, rSeen1);
jpayne@68 333 int maxDiv=Tools.max(qSeen1, rSeen1);
jpayne@68 334 long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
jpayne@68 335 return eValue(hits, minDiv, maxDiv, range);
jpayne@68 336 }
jpayne@68 337
jpayne@68 338 public double eValue2(){
jpayne@68 339 if(k2<1){return 1.0;}
jpayne@68 340
jpayne@68 341 long range0=range();
jpayne@68 342 int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k2;
jpayne@68 343 double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
jpayne@68 344 int hits=hits2();
jpayne@68 345 int minDiv=Tools.min(qSeen2(), rSeen2());
jpayne@68 346 int maxDiv=Tools.max(qSeen2(), rSeen2());
jpayne@68 347 long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
jpayne@68 348 // assert(false) : missingBits+", "+quantizer+", "+range0+", "+range+", "+eValue(hits, minDiv, maxDiv, range);
jpayne@68 349 return eValue(hits, minDiv, maxDiv, range);
jpayne@68 350 }
jpayne@68 351
jpayne@68 352 public String scoreS(){
jpayne@68 353 float x=score;
jpayne@68 354 return format3(x);
jpayne@68 355 }
jpayne@68 356
jpayne@68 357 public double depth(boolean observedToActual){
jpayne@68 358 return observedToActual ? Tools.observedToActualCoverage(depth) : depth;
jpayne@68 359 }
jpayne@68 360
jpayne@68 361 public double depth2(boolean observedToActual){
jpayne@68 362 return observedToActual ? Tools.observedToActualCoverage(depth2) : depth2;
jpayne@68 363 }
jpayne@68 364
jpayne@68 365 public String depthS(boolean observedToActual){
jpayne@68 366 float x=depth;
jpayne@68 367 if(observedToActual){x=(float)Tools.observedToActualCoverage(x);}
jpayne@68 368 return format3(x);
jpayne@68 369 }
jpayne@68 370
jpayne@68 371 public float avgRefHits(){
jpayne@68 372 return avgRefHits;
jpayne@68 373 }
jpayne@68 374
jpayne@68 375 public String avgRefHitsS(){
jpayne@68 376 return format2(avgRefHits);
jpayne@68 377 }
jpayne@68 378
jpayne@68 379 public String depth2S(boolean observedToActual){
jpayne@68 380 float x=depth2;
jpayne@68 381 if(observedToActual){
jpayne@68 382 x=(float)(Tools.observedToActualCoverage(depth)*(depth2/depth));
jpayne@68 383 }
jpayne@68 384 return format3(x);
jpayne@68 385 }
jpayne@68 386
jpayne@68 387 public String volumeS(){
jpayne@68 388 double x=volume()*0.001;
jpayne@68 389 return format3(x);
jpayne@68 390 }
jpayne@68 391
jpayne@68 392 static String format3(double x){
jpayne@68 393 if(x>=999.95){
jpayne@68 394 return(""+(long)Math.round(x));
jpayne@68 395 }else if(x>=99.995){
jpayne@68 396 return String.format(Locale.ROOT, "%.1f", x);
jpayne@68 397 }else if(x>=9.9995){
jpayne@68 398 return String.format(Locale.ROOT, "%.2f", x);
jpayne@68 399 }
jpayne@68 400 return String.format(Locale.ROOT, "%.3f", x);
jpayne@68 401 }
jpayne@68 402
jpayne@68 403 static String format2(double x){
jpayne@68 404 if(x>=999.95){
jpayne@68 405 return(""+(long)Math.round(x));
jpayne@68 406 }else if(x>=99.995){
jpayne@68 407 return String.format(Locale.ROOT, "%.1f", x);
jpayne@68 408 }else if(x>=9.9995){
jpayne@68 409 return String.format(Locale.ROOT, "%.2f", x);
jpayne@68 410 }
jpayne@68 411 return String.format(Locale.ROOT, "%.2f", x);
jpayne@68 412 }
jpayne@68 413
jpayne@68 414 float volume(){
jpayne@68 415 return Tools.max(1f, depth)*hits;
jpayne@68 416 }
jpayne@68 417
jpayne@68 418 @Override
jpayne@68 419 public String toString(){
jpayne@68 420 return "hits="+hits+", refDivisor="+refDivisor+", queryDivisor="+queryDivisor+", refSize="+refSize+", querySize="+querySize+
jpayne@68 421 ", contamHits="+contamHits+", contam2Hits="+contam2Hits+", multiContamHits="+multiContamHits+", depth="+depth+", depth2="+depth2+", volume="+volume()+
jpayne@68 422 ", hits="+hits+", multiHits="+multiHits+", unique2="+unique2+", unique3="+unique3+", noHits="+noHits+", taxID="+taxID+", taxName="+taxName;
jpayne@68 423 }
jpayne@68 424
jpayne@68 425 /*--------------------------------------------------------------*/
jpayne@68 426 /*---------------- Getters ----------------*/
jpayne@68 427 /*--------------------------------------------------------------*/
jpayne@68 428
jpayne@68 429 // public boolean passesFilter(TaxFilter white, TaxFilter black){
jpayne@68 430 // if(white==null && black==null){return true;}
jpayne@68 431 // int id=taxID();
jpayne@68 432 // String s=name();
jpayne@68 433 // return passesFilter(white, id, s) && passesFilter(black, id, s);
jpayne@68 434 // }
jpayne@68 435 //
jpayne@68 436 // private boolean passesFilter(TaxFilter filter, int id, String s){
jpayne@68 437 // if(filter==null){return true;}
jpayne@68 438 // if(id>0 && !filter.passesFilter(id)){return false;}
jpayne@68 439 // if(s!=null && !filter.passesFilterByNameOnly(s)){return false;}
jpayne@68 440 // return true;
jpayne@68 441 // }
jpayne@68 442
jpayne@68 443 public String name(){return taxName!=null ? taxName : name0()!=null ? name0() : fname()!=null ? fname() : ""+taxID();}
jpayne@68 444 public String taxName(){return taxName;}
jpayne@68 445 String name0(){return b.name0();}
jpayne@68 446 String fname(){return b.fname();}
jpayne@68 447
jpayne@68 448 // public int taxID(){return b.taxID<minFakeID ? b.taxID : 0;}
jpayne@68 449 public int taxID(){return (taxID<minFakeID && taxID>=0) ? taxID : -1;}
jpayne@68 450 long imgID(){return (b.imgID>0 ? b.imgID : -1);}
jpayne@68 451
jpayne@68 452 long genomeSizeBases(){return b.genomeSizeBases;}
jpayne@68 453 long genomeSizeKmers(){return b.genomeSizeKmers;}
jpayne@68 454 long genomeSequences(){return b.genomeSequences;}
jpayne@68 455 long genomeSizeEstimate(){return b.genomeSizeEstimate();}
jpayne@68 456 float gc(){return b.gc();}
jpayne@68 457 boolean hasGC(){return b.baseCounts!=null;}
jpayne@68 458
jpayne@68 459 public boolean hasSSU() {
jpayne@68 460 return (a.r16S()!=null && b.r16S()!=null) || (a.r18S()!=null && b.r18S()!=null);
jpayne@68 461 }
jpayne@68 462 public boolean hasSSUIdentity() {
jpayne@68 463 return ssuIdentity>=0;
jpayne@68 464 }
jpayne@68 465 public boolean needsAlignment() {
jpayne@68 466 return hasSSU() && !hasSSUIdentity();
jpayne@68 467 }
jpayne@68 468 public boolean hasQueryTaxID() {
jpayne@68 469 return a.taxID>0 && a.taxID<minFakeID;
jpayne@68 470 }
jpayne@68 471
jpayne@68 472
jpayne@68 473 public int uHits() {return hits-multiHits;}
jpayne@68 474
jpayne@68 475 /** Common ancestor TaxID, if both Sketches have a TaxID */
jpayne@68 476 public int commonAncestor() {
jpayne@68 477 if(a.taxID<1 || b.taxID<1){return -1;}
jpayne@68 478 assert(taxtree!=null);
jpayne@68 479 int id=taxtree.commonAncestor(a.taxID, b.taxID);
jpayne@68 480 return id;
jpayne@68 481 }
jpayne@68 482
jpayne@68 483 /** Common ancestor node tax level, if both Sketches have a TaxID */
jpayne@68 484 public String commonAncestorLevel() {
jpayne@68 485 int id=commonAncestor();
jpayne@68 486 if(id<1){return ".";}
jpayne@68 487 TaxNode tn=taxtree.getNode(id);
jpayne@68 488 while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
jpayne@68 489 String s=tn.levelStringExtended(false);
jpayne@68 490 return s;
jpayne@68 491 }
jpayne@68 492
jpayne@68 493 /** Common ancestor node tax level, if both Sketches have a TaxID */
jpayne@68 494 public int commonAncestorLevelInt() {
jpayne@68 495 int id=commonAncestor();
jpayne@68 496 if(id<1){return 0;}
jpayne@68 497 TaxNode tn=taxtree.getNode(id);
jpayne@68 498 while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
jpayne@68 499 return tn.levelExtended;
jpayne@68 500 }
jpayne@68 501
jpayne@68 502 /*--------------------------------------------------------------*/
jpayne@68 503 /*---------------- Comparators ----------------*/
jpayne@68 504 /*--------------------------------------------------------------*/
jpayne@68 505
jpayne@68 506
jpayne@68 507
jpayne@68 508 static class ScoreComparator implements Comparator<Comparison>{
jpayne@68 509
jpayne@68 510 @Override
jpayne@68 511 public int compare(Comparison a, Comparison b) {
jpayne@68 512 {
jpayne@68 513 float pa=a.score, pb=b.score;
jpayne@68 514 if(pa>pb){
jpayne@68 515 return 1;
jpayne@68 516 }else if (pa<pb){
jpayne@68 517 return -1;
jpayne@68 518 }
jpayne@68 519 }
jpayne@68 520
jpayne@68 521 int x=a.hits-b.hits;
jpayne@68 522 if(x!=0){return x;}
jpayne@68 523 x=b.minDivisor()-a.minDivisor();
jpayne@68 524 if(x!=0){return x;}
jpayne@68 525 x=b.maxDivisor()-a.maxDivisor();
jpayne@68 526 if(x!=0){return x;}
jpayne@68 527 x=b.refDivisor-a.refDivisor;
jpayne@68 528 if(x!=0){return x;}
jpayne@68 529 x=a.taxID()-b.taxID();
jpayne@68 530 if(x!=0){return x;}
jpayne@68 531 if(a.name0()!=null && b.name0()!=null){
jpayne@68 532 return a.name0().compareTo(b.name0());
jpayne@68 533 }
jpayne@68 534 if(a.taxName()!=null && b.taxName()!=null){
jpayne@68 535 return a.taxName().compareTo(b.taxName());
jpayne@68 536 }
jpayne@68 537 return 0;
jpayne@68 538 }
jpayne@68 539
jpayne@68 540 @Override
jpayne@68 541 public String toString(){return "sortByScore";}
jpayne@68 542
jpayne@68 543 }
jpayne@68 544
jpayne@68 545 static class DepthComparator implements Comparator<Comparison>{
jpayne@68 546
jpayne@68 547 @Override
jpayne@68 548 public int compare(Comparison a, Comparison b) {
jpayne@68 549 final float da=Tools.max(0.1f, a.depth-0.5f), db=Tools.max(0.1f, b.depth-0.5f);
jpayne@68 550 final float sa, sb;
jpayne@68 551 if(sqrt){
jpayne@68 552 sa=da*(float)Math.sqrt(a.score);
jpayne@68 553 sb=db*(float)Math.sqrt(b.score);
jpayne@68 554 }else{
jpayne@68 555 sa=da*a.score;
jpayne@68 556 sb=db*b.score;
jpayne@68 557 }
jpayne@68 558 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 559 }
jpayne@68 560
jpayne@68 561 @Override
jpayne@68 562 public String toString(){return "sortByDepth";}
jpayne@68 563
jpayne@68 564 }
jpayne@68 565
jpayne@68 566 static class Depth2Comparator implements Comparator<Comparison>{
jpayne@68 567
jpayne@68 568 @Override
jpayne@68 569 public int compare(Comparison a, Comparison b) {
jpayne@68 570 final float da=Tools.max(0.1f, a.depth2-0.8f), db=Tools.max(0.1f, b.depth2-0.8f);
jpayne@68 571 final float sa, sb;
jpayne@68 572 if(sqrt){
jpayne@68 573 sa=da*(float)Math.sqrt(a.score);
jpayne@68 574 sb=db*(float)Math.sqrt(b.score);
jpayne@68 575 }else{
jpayne@68 576 sa=da*a.score;
jpayne@68 577 sb=db*b.score;
jpayne@68 578 }
jpayne@68 579 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 580 }
jpayne@68 581
jpayne@68 582 @Override
jpayne@68 583 public String toString(){return "sortByDepth2";}
jpayne@68 584
jpayne@68 585 }
jpayne@68 586
jpayne@68 587 static class VolumeComparator implements Comparator<Comparison>{
jpayne@68 588
jpayne@68 589 @Override
jpayne@68 590 public int compare(Comparison a, Comparison b) {
jpayne@68 591 final float da=a.volume(), db=b.volume();
jpayne@68 592 final float sa, sb;
jpayne@68 593 if(sqrt){
jpayne@68 594 sa=da*(float)Math.sqrt(a.score);
jpayne@68 595 sb=db*(float)Math.sqrt(b.score);
jpayne@68 596 }else{
jpayne@68 597 sa=da*a.score;
jpayne@68 598 sb=db*b.score;
jpayne@68 599 }
jpayne@68 600 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 601 }
jpayne@68 602
jpayne@68 603 @Override
jpayne@68 604 public String toString(){return "sortByVolume";}
jpayne@68 605
jpayne@68 606 }
jpayne@68 607
jpayne@68 608 static class KIDComparator implements Comparator<Comparison>{
jpayne@68 609
jpayne@68 610 @Override
jpayne@68 611 public int compare(Comparison a, Comparison b) {
jpayne@68 612 final float sa=a.kid(), sb=b.kid();
jpayne@68 613 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 614 }
jpayne@68 615
jpayne@68 616 @Override
jpayne@68 617 public String toString(){return "sortByKID";}
jpayne@68 618
jpayne@68 619 }
jpayne@68 620
jpayne@68 621 static class WKIDComparator implements Comparator<Comparison>{
jpayne@68 622
jpayne@68 623 @Override
jpayne@68 624 public int compare(Comparison a, Comparison b) {
jpayne@68 625 final float sa=a.wkid(), sb=b.wkid();
jpayne@68 626 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 627 }
jpayne@68 628
jpayne@68 629 @Override
jpayne@68 630 public String toString(){return "sortByWKID";}
jpayne@68 631
jpayne@68 632 }
jpayne@68 633
jpayne@68 634 static class SSUComparator implements Comparator<Comparison>{
jpayne@68 635
jpayne@68 636 @Override
jpayne@68 637 public int compare(Comparison a, Comparison b) {
jpayne@68 638 // if((a.has16S() || a.has18S()) && !a.hasSSUIdentity()){
jpayne@68 639 // synchronized(a){a.ssuIdentity();}
jpayne@68 640 // assert(a.hasSSUIdentity());
jpayne@68 641 // }
jpayne@68 642 // if((b.has16S() || b.has18S()) && !b.hasSSUIdentity()){
jpayne@68 643 // synchronized(b){b.ssuIdentity();}
jpayne@68 644 // assert(b.hasSSUIdentity());
jpayne@68 645 // }
jpayne@68 646
jpayne@68 647 if(a.hasSSUIdentity() && b.hasSSUIdentity()){
jpayne@68 648 final float sa=a.ssuIdentity(), sb=b.ssuIdentity();
jpayne@68 649 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 650 }else if(a.hasSSUIdentity()){
jpayne@68 651 return 1;
jpayne@68 652 }else if(b.hasSSUIdentity()){
jpayne@68 653 return -1;
jpayne@68 654 }else{
jpayne@68 655 return scoreComparator.compare(a, b);
jpayne@68 656 }
jpayne@68 657
jpayne@68 658 }
jpayne@68 659
jpayne@68 660 @Override
jpayne@68 661 public String toString(){return "sortBySSU";}
jpayne@68 662
jpayne@68 663 }
jpayne@68 664
jpayne@68 665 static class HitsComparator implements Comparator<Comparison>{
jpayne@68 666
jpayne@68 667 @Override
jpayne@68 668 public int compare(Comparison a, Comparison b) {
jpayne@68 669 final float sa=a.hits(), sb=b.hits();
jpayne@68 670 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68 671 }
jpayne@68 672
jpayne@68 673 @Override
jpayne@68 674 public String toString(){return "sortByHits";}
jpayne@68 675
jpayne@68 676 }
jpayne@68 677
jpayne@68 678 @Override
jpayne@68 679 public int compareTo(Comparison b) {
jpayne@68 680 assert(false) : "Please use comparators instead.";
jpayne@68 681 return scoreComparator.compare(this, b);
jpayne@68 682 }
jpayne@68 683
jpayne@68 684 @Override
jpayne@68 685 public int hashCode() {
jpayne@68 686 assert(false) : "TODO";
jpayne@68 687 return super.hashCode();
jpayne@68 688 }
jpayne@68 689
jpayne@68 690 public static final ScoreComparator scoreComparator=new ScoreComparator();
jpayne@68 691 public static final DepthComparator depthComparator=new DepthComparator();
jpayne@68 692 public static final Depth2Comparator depth2Comparator=new Depth2Comparator();
jpayne@68 693 public static final VolumeComparator volumeComparator=new VolumeComparator();
jpayne@68 694 public static final KIDComparator KIDComparator=new KIDComparator();
jpayne@68 695 public static final SSUComparator SSUComparator=new SSUComparator();
jpayne@68 696 public static final WKIDComparator WKIDComparator=new WKIDComparator();
jpayne@68 697 public static final HitsComparator HitsComparator=new HitsComparator();
jpayne@68 698 private static final boolean sqrt=false;
jpayne@68 699 private static final double aaBitValue=0.86438561897747246957406388589788; //(log(20)/log(2))/5;
jpayne@68 700
jpayne@68 701 /*--------------------------------------------------------------*/
jpayne@68 702
jpayne@68 703 int hits(){return hits;}
jpayne@68 704 int multiHits(){return multiHits;}
jpayne@68 705 int noHits(){return noHits;}
jpayne@68 706 int unique2(){return unique2;}
jpayne@68 707 int unique3(){return unique3;}
jpayne@68 708
jpayne@68 709 float depth(){return depth;}
jpayne@68 710 float depth2(){return depth2;}
jpayne@68 711 float score(){return score;}
jpayne@68 712
jpayne@68 713 int contamHits(){return contamHits;}
jpayne@68 714 int contam2Hits(){return contam2Hits;}
jpayne@68 715 int multiContamHits(){return multiContamHits;}
jpayne@68 716
jpayne@68 717 int queryDivisor(){return queryDivisor;}
jpayne@68 718 int refDivisor(){return refDivisor;}
jpayne@68 719
jpayne@68 720 int querySize(){return querySize;}
jpayne@68 721 int refSize(){return refSize;}
jpayne@68 722
jpayne@68 723 int ssuType(){return has18S() ? 18 : has16S() ? 16 : 0;}
jpayne@68 724 int ssuLen(){return has18S() ? b.r18SLen() : has16S() ? b.r16SLen() : 0;}
jpayne@68 725 boolean has16S(){return a.r16S()!=null && b.r16S()!=null;}
jpayne@68 726 boolean has18S(){return a.r18S()!=null && b.r18S()!=null;}
jpayne@68 727
jpayne@68 728 float ssuIdentity(){
jpayne@68 729 if(ssuIdentity>0){return ssuIdentity;}
jpayne@68 730 if(has18S()){ssuIdentity=calcIdentity(a.r18S(), b.r18S());}
jpayne@68 731 else if(has16S()){ssuIdentity=calcIdentity(a.r16S(), b.r16S());}
jpayne@68 732 return ssuIdentity;
jpayne@68 733 }
jpayne@68 734
jpayne@68 735 private static float calcIdentity(byte[] ssuA, byte[] ssuB){
jpayne@68 736 //assert(false);
jpayne@68 737 if(ssuA.length>ssuB.length){
jpayne@68 738 byte[] c=ssuA;
jpayne@68 739 ssuA=ssuB;
jpayne@68 740 ssuB=c;
jpayne@68 741 }
jpayne@68 742 if(useSSA){
jpayne@68 743 Aligner ssa=(useSSA3 ? GeneCaller.getSSA3() : GeneCaller.getSSA());
jpayne@68 744 int[] max=ssa.fillUnlimited(ssuA, ssuB, 0, ssuB.length-1, 0);
jpayne@68 745 if(max==null){return 0;}
jpayne@68 746
jpayne@68 747 final int rows=max[0];
jpayne@68 748 final int maxCol=max[1];
jpayne@68 749 final int maxState=max[2];
jpayne@68 750
jpayne@68 751 //returns {score, bestRefStart, bestRefStop}
jpayne@68 752 //padded: {score, bestRefStart, bestRefStop, padLeft, padRight};
jpayne@68 753 int[] score=ssa.score(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
jpayne@68 754 int rstart=score[1];
jpayne@68 755 int rstop=score[2];
jpayne@68 756
jpayne@68 757 // byte[] match=ssa.traceback(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
jpayne@68 758 // float id=Read.identity(match);
jpayne@68 759 float id=ssa.tracebackIdentity(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState, null);
jpayne@68 760 return id;
jpayne@68 761 }else{
jpayne@68 762 GlocalAligner ga=new GlocalAligner();
jpayne@68 763 return ga.alignForward(ssuA, ssuB);
jpayne@68 764 }
jpayne@68 765 }
jpayne@68 766
jpayne@68 767 /*--------------------------------------------------------------*/
jpayne@68 768 /*---------------- Fields ----------------*/
jpayne@68 769 /*--------------------------------------------------------------*/
jpayne@68 770
jpayne@68 771 public Sketch a, b;
jpayne@68 772
jpayne@68 773 String taxName;
jpayne@68 774 int taxID;
jpayne@68 775
jpayne@68 776 private int hits;
jpayne@68 777 private int multiHits;
jpayne@68 778 private int unique2;
jpayne@68 779 private int unique3;
jpayne@68 780 private int noHits;
jpayne@68 781
jpayne@68 782 private float depth;
jpayne@68 783 private float depth2;
jpayne@68 784 private float avgRefHits;
jpayne@68 785 private float score;
jpayne@68 786
jpayne@68 787 private int contamHits;
jpayne@68 788 private int contam2Hits;
jpayne@68 789 private int multiContamHits;
jpayne@68 790
jpayne@68 791 private int refDivisor;
jpayne@68 792 private int queryDivisor;
jpayne@68 793
jpayne@68 794 private int refSize;
jpayne@68 795 private int querySize;
jpayne@68 796
jpayne@68 797 private int hits1;
jpayne@68 798 private int qSeen1;
jpayne@68 799 private int rSeen1;
jpayne@68 800
jpayne@68 801 private float ssuIdentity=-1;
jpayne@68 802
jpayne@68 803 }