X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_md.c;h=d2c84dcccd084c3399bf94ac6085af8902eb8532;hb=6cda28b3a4c6afcf8c739a478b6bcf1e514f96e9;hp=eee767e08fb290ab1ea8baae778201e797c923c3;hpb=51f014165fb63b2b9b86db80ea4681f9351453a1;p=samtools.git

diff --git a/bam_md.c b/bam_md.c
index eee767e..d2c84dc 100644
--- a/bam_md.c
+++ b/bam_md.c
@@ -110,14 +110,14 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal)
 	bam_fillmd1_core(b, ref, is_equal, 0);
 }
 
-int bam_cap_mapQ(bam1_t *b, char *ref)
+int bam_cap_mapQ(bam1_t *b, char *ref, int thres)
 {
 	uint8_t *seq = bam1_seq(b), *qual = bam1_qual(b);
 	uint32_t *cigar = bam1_cigar(b);
 	bam1_core_t *c = &b->core;
 	int i, x, y, mm, q, len, clip_l, clip_q;
 	double t;
-
+	if (thres < 0) thres = 40; // set the default
 	mm = q = len = clip_l = clip_q = 0;
 	for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) {
 		int j, l = cigar[i]>>4, op = cigar[i]&0xf;
@@ -154,10 +154,10 @@ int bam_cap_mapQ(bam1_t *b, char *ref)
 	for (i = 0, t = 1; i < mm; ++i)
 		t *= (double)len / (i+1);
 	t = q - 4.343 * log(t) + clip_q / 5.;
-	if (t > 40) return -1;
+	if (t > thres) return -1;
 	if (t < 0) t = 0;
-	t = sqrt((40 - t) / 40) * 40;
-	fprintf(stderr, "%s %lf %d\n", bam1_qname(b), t, q);
+	t = sqrt((thres - t) / thres) * thres;
+//	fprintf(stderr, "%s %lf %d\n", bam1_qname(b), t, q);
 	return (int)(t + .499);
 }
 
@@ -320,7 +320,7 @@ int bam_fillmd(int argc, char *argv[])
 			}
 			if (is_realn) bam_realn(b, ref);
 			if (is_capQ) {
-				int q = bam_cap_mapQ(b, ref);
+				int q = bam_cap_mapQ(b, ref, 40);
 				if (b->core.qual > q) b->core.qual = q;
 			}
 			if (ref) bam_fillmd1_core(b, ref, is_equal, max_nm);