* make tview more friendly

[samtools.git] / bam_maqcns.c

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 9a25380f36b77411ecc34336b5a92fe2d2f734d4..cad63d772268db289467c6a4fd625e50b7ac3889 100644 (file)
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -310,6 +310,7 @@ bam_maqindel_opt_t *bam_maqindel_opt_init()
        bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
        mi->q_indel = 40;
        mi->r_indel = 0.00015;
+       mi->r_snp = 0.001;
        //
        mi->mm_penalty = 3;
        mi->indel_err = 4;
@@ -406,7 +407,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
        }
        { // the core part
                char *ref2, *rs, *inscns = 0;
-               int k, l, *score, *pscore, max_ins = types[n_types-1];
+               int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1];
+               qr_snp = (int)(-4.343 * log(mi->r_snp) + .499);
                if (max_ins > 0) { // get the consensus of inserted sequences
                        int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
                        // count occurrences
@@ -488,7 +490,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c
                                                if (op == BAM_CMATCH) {
                                                        int k;
                                                        for (k = 0; k < len; ++k)
-                                                               if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4) ps += bam1_qual(p->b)[y+k];
+                                                               if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4)
+                                                                       ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp;
                                                        x += len; y += len;
                                                } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
                                                        if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;