#include "bam.h"
#include "kstring.h"
#include "bam2bcf.h"
-#include "bcf.h"
+#include "bcftools/bcf.h"
+
+#define MAX_END_DIST 30
extern void ks_introsort_uint32_t(size_t n, uint32_t a[]);
#define CALL_DEFTHETA 0.85f
struct __bcf_callaux_t {
- int max_info, capQ;
+ int max_info, capQ, min_baseQ;
double *fk;
uint32_t *info;
};
-bcf_callaux_t *bcf_call_init(double theta)
+bcf_callaux_t *bcf_call_init(double theta, int min_baseQ)
{
bcf_callaux_t *bca;
int n;
if (theta <= 0.) theta = CALL_DEFTHETA;
bca = calloc(1, sizeof(bcf_callaux_t));
bca->capQ = 60;
+ bca->min_baseQ = min_baseQ;
bca->fk = calloc(CALL_MAX, sizeof(double));
bca->fk[0] = 1.;
for (n = 1; n < CALL_MAX; ++n)
}
typedef struct {
- float esum[4], fsum[4];
- uint32_t c[4];
+ float esum[5], fsum[5];
+ uint32_t c[5];
int w[8];
} auxaux_t;
*/
int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf_callaux_t *bca, bcf_callret1_t *r)
{
- int i, j, k, c, n;
+ int i, j, k, c, n, ref4;
float *p = r->p;
auxaux_t aux;
memset(r, 0, sizeof(bcf_callret1_t));
+ ref4 = bam_nt16_nt4_table[ref_base];
if (_n == 0) return -1;
// enlarge the aux array if necessary
for (i = n = 0; i < _n; ++i) {
const bam_pileup1_t *p = pl + i;
uint32_t q, x = 0, qq;
+ int min_dist;
if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue; // skip unmapped reads and deleted bases
q = (uint32_t)bam1_qual(p->b)[p->qpos]; // base quality
+ if (q < bca->min_baseQ) continue;
x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
if (p->b->core.qual < q) q = p->b->core.qual; // cap the overall quality at mapping quality
x |= q << 24;
qq = bam1_seqi(bam1_seq(p->b), p->qpos); // base
q = bam_nt16_nt4_table[qq? qq : ref_base]; // q is the 2-bit base
if (q < 4) x |= 1 << 21 | q << 16;
-
+ k = (ref4 < 4 && q == ref4)? 0 : 1;
+ k = k<<1 | bam1_strand(p->b);
+ ++r->d[k];
bca->info[n++] = x;
+ // calculate min_dist
+ min_dist = p->b->core.l_qseq - 1 - p->qpos;
+ if (min_dist > p->qpos) min_dist = p->qpos;
+ if (min_dist > MAX_END_DIST) min_dist = MAX_END_DIST;
+ k >>= 1;
+ r->dsum[k] += min_dist;
+ r->d2sum[k] += min_dist * min_dist;
}
ks_introsort_uint32_t(n, bca->info);
r->depth = n;
tmp = (int)(info&0xff) < bca->capQ? (int)(info&0xff) : bca->capQ;
r->sum_Q2 += tmp * tmp;
}
- memcpy(r->esum, aux.esum, 4 * sizeof(float));
+ memcpy(r->esum, aux.esum, 5 * sizeof(float));
// rescale ->c[]
for (j = c = 0; j != 4; ++j) c += aux.c[j];
if (c > 255) {
for (j = c = 0; j != 4; ++j) c += aux.c[j];
}
// generate likelihood
- for (j = 0; j != 4; ++j) {
+ for (j = 0; j != 5; ++j) {
float tmp;
// homozygous
- for (k = 0, tmp = 0.0; k != 4; ++k)
+ for (k = 0, tmp = 0.0; k != 5; ++k)
if (j != k) tmp += aux.esum[k];
- p[j<<2|j] = tmp; // anything that is not j
+ p[j*5+j] = tmp; // anything that is not j
// heterozygous
- for (k = j + 1; k < 4; ++k) {
- for (i = 0, tmp = 0.0; i != 4; ++i)
+ for (k = j + 1; k < 5; ++k) {
+ for (i = 0, tmp = 0.0; i != 5; ++i)
if (i != j && i != k) tmp += aux.esum[i];
- p[j<<2|k] = p[k<<2|j] = 3.01 * (aux.c[j] + aux.c[k]) + tmp;
+ p[j*5+k] = p[k*5+j] = 3.01 * (aux.c[j] + aux.c[k]) + tmp;
}
}
- return 0;
+ return r->depth;
}
-/*
- 1) Find the top 2 bases (from esum[4]).
-
- 2) If the reference base is among the top 2, consider the locus is
- potentially biallelic and set call->a[2] as -1; otherwise, the
- locus is potentially triallelic. If the reference is ambiguous,
- take the weakest call as the pseudo-reference.
- */
int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call)
{
int ref4, i, j;
- int64_t sum[4], tmp;
+ int64_t sum[5], tmp;
call->ori_ref = ref4 = bam_nt16_nt4_table[ref_base];
if (ref4 > 4) ref4 = 4;
{ // calculate esum
- double esum[4];
+ double esum[5];
memset(esum, 0, sizeof(double) * 4);
for (i = 0; i < n; ++i) {
for (j = 0; j < 4; ++j)
for (j = i; j > 0 && sum[j] < sum[j-1]; --j)
tmp = sum[j], sum[j] = sum[j-1], sum[j-1] = tmp;
// set the reference allele and alternative allele(s)
- for (i = 0; i < 4; ++i) call->a[i] = -1;
+ for (i = 0; i < 5; ++i) call->a[i] = -1;
call->unseen = -1;
- if (ref4 < 4) {
- call->a[0] = ref4;
- for (i = 3, j = 1; i >= 0; --i) {
- if ((sum[i]&3) != ref4) {
- if (sum[i]>>2 != 0) call->a[j++] = sum[i]&3;
- else break;
- }
- }
- if (j < 4 && i >= 0) call->unseen = j, call->a[j++] = sum[i]&3;
- call->n_alleles = j;
- } else {
- for (i = 3, j = 0; i >= 0; --i)
+ call->a[0] = ref4;
+ for (i = 3, j = 1; i >= 0; --i) {
+ if ((sum[i]&3) != ref4) {
if (sum[i]>>2 != 0) call->a[j++] = sum[i]&3;
else break;
- if (j < 4 && i >= 0) call->unseen = j, call->a[j++] = sum[i]&3;
- call->n_alleles = j;
+ }
}
+ if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
+ call->unseen = j, call->a[j++] = sum[i]&3;
+ call->n_alleles = j;
// set the PL array
if (call->n < n) {
call->n = n;
- call->PL = realloc(call->PL, 6 * n);
+ call->PL = realloc(call->PL, 15 * n);
}
{
- int x, g[6], z;
+ int x, g[15], z;
double sum_min = 0.;
x = call->n_alleles * (call->n_alleles + 1) / 2;
// get the possible genotypes
for (i = z = 0; i < call->n_alleles; ++i)
for (j = i; j < call->n_alleles; ++j)
- g[z++] = call->a[i]<<2 | call->a[j];
+ g[z++] = call->a[i] * 5 + call->a[j];
for (i = 0; i < n; ++i) {
uint8_t *PL = call->PL + x * i;
const bcf_callret1_t *r = calls + i;
}
call->shift = (int)(sum_min + .499);
}
+ memset(call->d, 0, 4 * sizeof(int));
+ call->davg[0] = call->davg[1] = call->dstd[0] = call->dstd[1] = 0.;
for (i = call->depth = 0, tmp = 0; i < n; ++i) {
call->depth += calls[i].depth;
+ for (j = 0; j < 4; ++j) call->d[j] += calls[i].d[j];
+ for (j = 0; j < 2; ++j) {
+ call->davg[j] += calls[i].dsum[j];
+ call->dstd[j] += calls[i].d2sum[j];
+ }
tmp += calls[i].sum_Q2;
}
+ for (j = 0; j < 2; ++j) {
+ int x = call->d[j*2] + call->d[j*2+1];
+ if (x == 0) {
+ call->davg[j] = call->dstd[j] = 0.;
+ } else {
+ call->davg[j] /= x;
+ call->dstd[j] = sqrt(call->dstd[j] / x - call->davg[j] * call->davg[j]);
+ }
+ }
call->rmsQ = (int)(sqrt((double)tmp / call->depth) + .499);
return 0;
}
int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b)
{
kstring_t s;
- int i, beg;
+ int i;
b->tid = tid; b->pos = pos; b->qual = 0;
s.s = b->str; s.m = b->m_str; s.l = 0;
kputc('\0', &s);
kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s);
- beg = bc->ori_ref > 3? 0 : 1;
- for (i = beg; i < 4; ++i) {
+ for (i = 1; i < 5; ++i) {
if (bc->a[i] < 0) break;
- if (i > beg) kputc(',', &s);
+ if (i > 1) kputc(',', &s);
// kputc(bc->unseen == i && i != 3? 'X' : "ACGT"[bc->a[i]], &s);
kputc(bc->unseen == i? 'X' : "ACGT"[bc->a[i]], &s);
}
kputc('\0', &s);
kputc('\0', &s);
- kputs("MQ=", &s); kputw(bc->rmsQ, &s); kputs(";DP=", &s); kputw(bc->depth, &s); kputc('\0', &s);
+ // INFO
+ kputs("MQ=", &s); kputw(bc->rmsQ, &s); // kputs(";DP=", &s); kputw(bc->depth, &s);
+ kputs(";DP4=", &s);
+ for (i = 0; i < 4; ++i) {
+ if (i) kputc(',', &s);
+ kputw(bc->d[i], &s);
+ }
+ if (bc->d[2] + bc->d[3] > 1)
+ ksprintf(&s, ";MED=%.1lf,%.1lf,%.1lf,%.1lf", bc->davg[0], bc->dstd[0], bc->davg[1], bc->dstd[1]);
+ kputc('\0', &s);
+ // FMT
kputs("PL", &s); kputc('\0', &s);
b->m_str = s.m; b->str = s.s; b->l_str = s.l;
bcf_sync(bc->n, b);
projects / samtools.git / blobdiff