X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam2bcf.c;h=1f7a0c042b2dad49abe5bb9d1668d90ddc96edf5;hb=f27c5d351576e00cece3708adcb386da61fe7b7f;hp=a78e7d282c8ff2280a819ca32570a107a076d32a;hpb=67184476fcd4a7aba449556a2333586465fd6102;p=samtools.git

diff --git a/bam2bcf.c b/bam2bcf.c
index a78e7d2..1f7a0c0 100644
--- a/bam2bcf.c
+++ b/bam2bcf.c
@@ -3,7 +3,9 @@
 #include "bam.h"
 #include "kstring.h"
 #include "bam2bcf.h"
-#include "bcf.h"
+#include "bcftools/bcf.h"
+
+#define END_DIST_THRES 11
 
 extern	void ks_introsort_uint32_t(size_t n, uint32_t a[]);
 
@@ -12,18 +14,19 @@ 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)
@@ -40,8 +43,8 @@ void bcf_call_destroy(bcf_callaux_t *bca)
 }
 
 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;
 
@@ -70,11 +73,12 @@ typedef struct {
  */
 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
@@ -87,16 +91,25 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf
 	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;
+		k = (k&2) | (min_dist <= END_DIST_THRES);
+		++r->ed[k];
 	}
 	ks_introsort_uint32_t(n, bca->info);
 	r->depth = n;
@@ -116,7 +129,7 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf
 		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) {
@@ -124,38 +137,30 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf
 		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)
@@ -169,38 +174,31 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/,
 		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;
@@ -217,8 +215,11 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/,
 		}
 		call->shift = (int)(sum_min + .499);
 	}
+	memset(call->d, 0, 4 * sizeof(int));
+	memset(call->ed, 0, 4 * sizeof(int));
 	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], call->ed[j] += calls[i].ed[j];
 		tmp += calls[i].sum_Q2;
 	}
 	call->rmsQ = (int)(sqrt((double)tmp / call->depth) + .499);
@@ -228,21 +229,33 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/,
 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);
+	}
+	kputs(";ED4=", &s);
+	for (i = 0; i < 4; ++i) {
+		if (i) kputc(',', &s);
+		kputw(bc->ed[i], &s);
+	}
+	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);