typedef struct {
float esum[4], fsum[4];
uint32_t c[4];
- uint32_t mapQ_max;
+ uint32_t rms_mapQ;
} glf_call_aux_t;
char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
int i, j, k, w[8], c, n;
glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
float p[16], min_p = 1e30;
+ uint64_t rms;
g->ref_base = ref_base;
if (_n == 0) return g;
// generate esum and fsum
b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
for (k = 0; k != 8; ++k) w[k] = 0;
- b->mapQ_max = 0;
+ rms = 0;
for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
uint32_t info = bm->aux->info[j];
if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
if (w[k] < 0xff) ++w[k];
++b->c[k&3];
}
- if (b->mapQ_max < (info&0x7f)) b->mapQ_max = info&0x7f;
+ rms += (int)(info&0x7f) * (info&0x7f);
}
+ b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
// rescale ->c[]
for (j = c = 0; j != 4; ++j) c += b->c[j];
if (c > 255) {
}
// convert necessary information to glf1_t
- g->ref_base = ref_base; g->max_mapQ = b->mapQ_max;
+ g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
g->depth = n > 16777215? 16777215 : n;
for (j = 0; j != 4; ++j)
for (k = j; k < 4; ++k)
#define MINUS_CONST 0x10000000
-bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref)
+bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
+ int _n_types, int *_types)
{
int i, j, n_types, *types, left, right;
bam_maqindel_ret_t *ret = 0;
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pl + i;
- if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
+ // if there is no proposed indel, check if there is an indel from the alignment
+ if (_n_types == 0) {
+ for (i = 0; i < n; ++i) {
+ const bam_pileup1_t *p = pl + i;
+ if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
+ }
+ if (i == n) return 0; // no indel
}
- if (i == n) return 0; // no indel
{ // calculate how many types of indels are available (set n_types and types)
int m;
uint32_t *aux;
- aux = (uint32_t*)calloc(n+1, 4);
+ aux = (uint32_t*)calloc(n + _n_types + 1, 4);
m = 0;
aux[m++] = MINUS_CONST; // zero indel is always a type
for (i = 0; i < n; ++i) {
if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
aux[m++] = MINUS_CONST + p->indel;
}
+ if (_n_types) // then also add this to aux[]
+ for (i = 0; i < _n_types; ++i)
+ if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
ks_introsort(uint32_t, m, aux);
- n_types = 1;
- for (i = 1; i < m; ++i)
+ // squeeze out identical types
+ for (i = 1, n_types = 1; i < m; ++i)
if (aux[i] != aux[i-1]) ++n_types;
types = (int*)calloc(n_types, sizeof(int));
j = 0;