- for (i = 0; i < b->n_alleles; ++i)
- p[i] += (int)pi[(i+1)*(i+2)/2-1];
- }
- for (i = 0; i < b->n_alleles; ++i) p[i] = p[i]<<4 | i;
- for (i = 1; i < b->n_alleles; ++i) // insertion sort
- for (j = i; j > 0 && p[j] < p[j-1]; --j)
- tmp = p[j], p[j] = p[j-1], p[j-1] = tmp;
- for (i = b->n_alleles - 1; i >= 0; --i)
- if ((p[i]&0xf) == 0) break;
- return i;
+ for (i=0; i<n; i++) lk[i] += pi[i];
+ }
+
+ double norm=lk[0];
+ for (i=1; i<n; i++) if (lk[i]<norm) norm=lk[i];
+ #if DBG
+ for (i=0,j=0; i<b->n_alleles; i++)
+ {
+ int k; for (k=0; k<=i; k++) printf("%.0f\t", lk[j++]);
+ printf("\n");
+ }
+ #endif
+ for (i=0; i<n; i++) lk[i] = pow(10,-0.1*(lk[i]-norm));
+
+ // Find out the most likely alleles. In contrast to the original version,
+ // ALT alleles may not be printed when they are more likely than REF but
+ // significantly less likely than the most likely ALT. The only criterion
+ // is the LK ratio now. To obtain behaviour similar to the original one,
+ // use the pref variant below.
+ double pmax=0; //,pref=0;
+ n = ma->is_indel ? b->n_alleles : b->n_alleles-1;
+ for (i=0; i<n; i++)
+ {
+ double pr=0;
+ int k=i*(i+1)/2;
+ for (j=0; j<=i; j++) { pr+=lk[k]; k++; }
+ for (j=i+1; j<b->n_alleles; j++) { k=j*(j+1)/2+i; pr+=lk[k]; }
+ #if DBG
+ printf("%d\t%e\n", i,pr);
+ #endif
+ if (pmax<pr) pmax=pr;
+ // if (i==0) pref=pr;
+ // if (pr<pref && pr/pmax < 1e-4) break;
+ if (pr/pmax < 1e-4) break; // Assuming the alleles are sorted by the lk
+ }
+ return i-1;