Make --indel-bias more sensitive to indels

by switching off an existing heuristics which reduces indelQ
and, in the extreme case, has the power to discard an indel
entirely when the reference alignment has low score.

This is a problem for long reads, so newly `--indel-bias 1.01`
is added to 'ont' and 'pacbio-ccs' presets.

Tentative fix to #1459

Author: pd3

bam2bcf.h

@@ -1,7 +1,8 @@
 /*  bam2bcf.h -- variant calling.
+                mplp.indel_bias = 1.01;
 
     Copyright (C) 2010-2012 Broad Institute.
-    Copyright (C) 2012-2021 Genome Research Ltd.
+    Copyright (C) 2012-2022 Genome Research Ltd.
 
     Author: Heng Li <[email protected]>
 
@@ -99,7 +100,7 @@ typedef struct __bcf_callaux_t {
     uint16_t *bases;        // 5bit: unused, 6:quality, 1:is_rev, 4:2-bit base or indel allele (index to bcf_callaux_t.indel_types)
     errmod_t *e;
     void *rghash;
-    float indel_bias;  // adjusts indel score threshold; lower => call more.
+    float indel_bias_inverted;  // adjusts indel score threshold, 1/--indel-bias, so lower => call more.
 } bcf_callaux_t;
 
 // per-sample values

bam2bcf_indel.c

@@ -1,7 +1,7 @@
 /*  bam2bcf_indel.c -- indel caller.
 
     Copyright (C) 2010, 2011 Broad Institute.
-    Copyright (C) 2012-2014,2016-2017, 2021 Genome Research Ltd.
+    Copyright (C) 2012-2014,2016-2017,2021-2022 Genome Research Ltd.
 
     Author: Heng Li <[email protected]>
 
@@ -540,7 +540,7 @@ static int bcf_cgp_align_score(bam_pileup1_t *p, bcf_callaux_t *bca,
     }
 
     // used for adjusting indelQ below
-    l = (int)(100. * sc / (qend - qbeg) + .499) * bca->indel_bias;
+    l = (int)((100. * sc / (qend - qbeg) + .499) * bca->indel_bias_inverted);
     *score = sc<<8 | MIN(255, l);
 
     rep_ele *reps, *elt, *tmp;
@@ -623,8 +623,14 @@ static int bcf_cgp_compute_indelQ(int n, int *n_plp, bam_pileup1_t **plp,
                 seqQ = est_seqQ(bca, types[sc[0]&0x3f], l_run);
             }
             tmp = sc[0]>>6 & 0xff;
+
+            // Don't know how this indelQ reduction threshold of 111 was derived,
+            // but it does not function well for longer reads that span multiple
+            // events.
+            //
             // reduce indelQ
-            indelQ = tmp > 111? 0 : (int)((1. - tmp/111.) * indelQ + .499);
+            if ( bca->indel_bias_inverted >= 1 )
+                indelQ = tmp > 111? 0 : (int)((1. - tmp/111.) * indelQ + .499);
 
             // Doesn't really help accuracy, but permits -h to take
             // affect still.
@@ -633,7 +639,7 @@ static int bcf_cgp_compute_indelQ(int n, int *n_plp, bam_pileup1_t **plp,
             if (seqQ > 255) seqQ = 255;
             p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ; // use 22 bits in total
             sumq[sc[0]&0x3f] += indelQ < seqQ? indelQ : seqQ;
-            //              fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d indelQ=%d seqQ=%d\n", pos, s, i, bam1_qname(p->b), types[sc[0]&0x3f], indelQ, seqQ);
+            // fprintf(stderr, "  read=%d:%d name=%s call=%d indelQ=%d seqQ=%d\n", s, i, bam_get_qname(p->b), types[sc[0]&0x3f], indelQ, seqQ);
         }
     }
     // determine bca->indel_types[] and bca->inscns
@@ -922,7 +928,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos,
                 fprintf(stderr, "pos=%d type=%d read=%d:%d name=%s "
                         "qbeg=%d tbeg=%d score=%d\n",
                         pos, types[t], s, i, bam_get_qname(p->b),
-                        qbeg, tbeg, sc);
+                        qbeg, tbeg, score[K*n_types + t]);
 #endif
             }
         }

doc/bcftools.txt

@@ -2150,8 +2150,8 @@ INFO/DPR    .. Deprecated in favor of INFO/AD; Number of high-quality bases for
 
         1.12           -Q13 -h100 -m1
         illumina       [ default values ]
-        ont	           -B -Q5 --max-BQ 30 -I
-        pacbio-ccs     -D -Q5 --max-BQ 50 -F0.1 -o25 -e1 -M99999
+        ont	           -B -Q5 --max-BQ 30 --indel-bias 1.01 -I
+        pacbio-ccs     -D -Q5 --max-BQ 50 --indel-bias 1.01 -F0.1 -o25 -e1 -M99999
 
 *--ar, --ambig-reads* 'drop'|'incAD'|'incAD0'::
     What to do with ambiguous indel reads that do not span an entire

mpileup.c

@@ -842,7 +842,7 @@ static int mpileup(mplp_conf_t *conf)
                               conf->delta_baseQ);
     conf->bcr = (bcf_callret1_t*) calloc(nsmpl, sizeof(bcf_callret1_t));
     conf->bca->openQ = conf->openQ, conf->bca->extQ = conf->extQ, conf->bca->tandemQ = conf->tandemQ;
-    conf->bca->indel_bias = conf->indel_bias;
+    conf->bca->indel_bias_inverted = 1 / conf->indel_bias;
     conf->bca->min_frac = conf->min_frac;
     conf->bca->min_support = conf->min_support;
     conf->bca->per_sample_flt = conf->flag & MPLP_PER_SAMPLE;
@@ -1180,8 +1180,8 @@ static void print_usage(FILE *fp, const mplp_conf_t *mplp)
         "Configuration profiles activated with -X, --config:\n"
         "    1.12:        -Q13 -h100 -m1 -F0.002\n"
         "    illumina:    [ default values ]\n"
-        "    ont:         -B -Q5 --max-BQ 30 -I [also try eg |bcftools call -P0.01]\n"
-        "    pacbio-ccs:  -D -Q5 --max-BQ 50 -F0.1 -o25 -e1 --delta-BQ 10 -M99999\n"
+        "    ont:         -B -Q5 --max-BQ 30 --indel-bias 1.01 -I [also try eg |bcftools call -P0.01]\n"
+        "    pacbio-ccs:  -D -Q5 --max-BQ 50 --indel-bias 1.01 -F0.1 -o25 -e1 --delta-BQ 10 -M99999\n"
         "\n"
         "Notes: Assuming diploid individuals.\n"
         "\n"
@@ -1400,11 +1400,9 @@ int main_mpileup(int argc, char *argv[])
             break;
         case 'e': mplp.extQ = atoi(optarg); break;
         case 'h': mplp.tandemQ = atoi(optarg); break;
-        case 10: // --indel-bias (inverted so higher => more indels called)
-            if (atof(optarg) < 1e-2)
-                mplp.indel_bias = 1/1e2;
-            else
-                mplp.indel_bias = 1/atof(optarg);
+        case 10: // --indel-bias
+            mplp.indel_bias =  atof(optarg);
+            if ( mplp.indel_bias < 1e-2 ) mplp.indel_bias = 1e-2;
             break;
         case  15: {
                 char *tmp;
@@ -1441,13 +1439,15 @@ int main_mpileup(int argc, char *argv[])
                 mplp.extQ = 1;
                 mplp.flag |= MPLP_REALN_PARTIAL;
                 mplp.max_read_len = 99999;
+                mplp.indel_bias = 1.01;
             } else if (strcasecmp(optarg, "ont") == 0) {
                 fprintf(stderr, "For ONT it may be beneficial to also run bcftools call with "
                         "a higher -P, eg -P0.01 or -P 0.1\n");
                 mplp.min_baseQ = 5;
                 mplp.max_baseQ = 30;
                 mplp.flag &= ~MPLP_REALN;
                 mplp.flag |= MPLP_NO_INDEL;
+                mplp.indel_bias = 1.01;
             } else if (strcasecmp(optarg, "1.12") == 0) {
                 // 1.12 and earlier
                 mplp.min_frac = 0.002;