Add a new indel caller.

This is a combination of PR samtools#1679 from 2022 coupled with more recent
changes to replace BAQ with edlib for evaluating the alignment of
reads against candidate alleles.  Due to the complexity of rebasing
many commits with many conflicts, these have been squashed together
for simplicity.  Please consult the original commit log in the PR for
more detailed history.

PR 1679 is a major restructuring of the indel caller that produces
proper consensus sequences.  TODO: paste collected commit messages
here.

The newer edlib changes are designed to dramtically increase the speed
of indel calling on long read data by use of a modern algorithm.  Note
edlib's alignments are very crude, with just +1 for all differences
including substitutions and each indel base.  However if we've
computed our candidate allele consensus sequences correctly then the
aligning against the real allele should be a minimal score regardless
of scoring regime, so the cost of an HMM or even affine weights do not
help much.

The lack of quality values is compensated for by detection of minimum
quality within STRs.

Author: jkbonfield

LICENSE

@@ -772,3 +772,28 @@ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+-----------------------------------------------------------------------------
+
+License for edlib.[ch]
+
+The MIT License (MIT)
+
+Copyright (c) 2014 Martin Šošić
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+the Software, and to permit persons to whom the Software is furnished to do so,
+subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Makefile

@@ -40,9 +40,10 @@ OBJS     = main.o vcfindex.o tabix.o \
            vcfcall.o mcall.o vcmp.o gvcf.o reheader.o convert.o vcfconvert.o tsv2vcf.o \
            vcfcnv.o vcfhead.o HMM.o consensus.o ploidy.o bin.o hclust.o version.o \
            regidx.o smpl_ilist.o csq.o vcfbuf.o \
-           mpileup.o bam2bcf.o bam2bcf_indel.o bam2bcf_iaux.o read_consensus.o bam_sample.o \
+           mpileup.o bam2bcf.o bam2bcf_indel.o bam2bcf_iaux.o bam2bcf_edlib.o \
+	   read_consensus.o bam_sample.o \
            vcfsort.o cols.o extsort.o dist.o abuf.o \
-           ccall.o em.o prob1.o kmin.o str_finder.o gff.o
+           ccall.o em.o prob1.o kmin.o str_finder.o gff.o edlib.o
 PLUGIN_OBJS = vcfplugin.o
 
 prefix      = /usr/local
@@ -234,6 +235,7 @@ vcfbuf_h = vcfbuf.h $(htslib_vcf_h)
 abuf_h = abuf.h $(htslib_vcf_h)
 dbuf_h = dbuf.h $(htslib_vcf_h)
 bam2bcf_h = bam2bcf.h $(htslib_hts_h) $(htslib_vcf_h)
+edlib.h = edlib.h
 bam_sample_h = bam_sample.h $(htslib_sam_h)
 cigar_state_h = cigar_state.h $(htslib_hts_h) $(htslib_sam_h)
 read_consensus_h = read_consensus.h $(htslib_hts_h) $(htslib_sam_h)
@@ -285,6 +287,7 @@ mpileup.o: mpileup.c $(htslib_sam_h) $(htslib_faidx_h) $(htslib_kstring_h) $(hts
 bam2bcf.o: bam2bcf.c $(htslib_hts_h) $(htslib_sam_h) $(htslib_kstring_h) $(htslib_kfunc_h) $(bam2bcf_h) mw.h
 bam2bcf_indel.o: bam2bcf_indel.c $(htslib_hts_h) $(htslib_sam_h) $(htslib_khash_str2int_h) $(bam2bcf_h) $(htslib_ksort_h) $(str_finder_h)
 bam2bcf_iaux.o: bam2bcf_iaux.c $(htslib_hts_h) $(htslib_sam_h) $(htslib_khash_str2int_h) $(bcftools_h) $(bam2bcf_h) $(htslib_ksort_h) $(read_consensus_h) $(cigar_state_h)
+bam2bcf_edlib.o: bam2bcf_edlib.c $(htslib_hts_h) $(htslib_sam_h) $(htslib_khash_str2int_h) $(bcftools_h) $(bam2bcf_h) $(htslib_ksort_h) $(read_consensus_h) $(cigar_state_h) $(edlib.h)
 read_consensus.o: read_consensus.c $(read_consensus_h) $(cigar_state_h) $(bcftools_h) kheap.h
 bam_sample.o: bam_sample.c $(htslib_hts_h) $(htslib_kstring_h) $(htslib_khash_str2int_h) $(khash_str2str_h) $(bam_sample_h) $(bcftools_h)
 version.o: version.h version.c

bam2bcf.c

@@ -268,6 +268,15 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
         bca->bases = (uint16_t*)realloc(bca->bases, 2 * bca->max_bases);
     }
 
+    // Detect if indel occurs anywhere in this sample
+    int indel_in_sample = 0;
+    if (bca->edlib) {
+        for (i = n = 0; i < _n; ++i) {
+            const bam_pileup1_t *p = pl + i;
+            if (p->indel) indel_in_sample = 1;
+        }
+    }
+
     // fill the bases array
     double nqual_over_60 = bca->nqual / 60.0;
     int ADR_ref_missed[4] = {0};
@@ -298,7 +307,19 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
             b = p->aux>>16&0x3f;        // indel type
             seqQ = q = (p->aux & 0xff); // mp2 + builtin indel-bias
 
-            if ( !bca->indels_v20 )
+            if (bca->edlib) {
+                if (indel_in_sample) {
+                    seqQ = q = (p->aux & 0xff); // mp2 + builtin indel-bias
+                } else {
+                    // An indel in another sample, but not this.  So just use
+                    // basic sequence confidences.
+                    q = bam_get_qual(p->b)[p->qpos];
+                    if (q > bca->max_baseQ) q = bca->max_baseQ;
+                    seqQ = 99;
+                }
+            }
+
+            if ( !bca->indels_v20 && !bca->edlib )
             {
                 /*
                     This heuristics was introduced by e4e161068 and claims to fix #1446. However, we obtain
@@ -341,6 +362,25 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
                 }
                 continue;
             }
+
+            // FIXME: CHECK if this is still needed with edlib mode
+            // It's a slight variant on the one above guarded by --indels-2.0
+            if (bca->edlib) {
+                if (indel_in_sample && p->indel == 0 && (q < _n/2 || _n > 20)) {
+                    // high quality indel calls without p->indel set aren't
+                    // particularly indicative of being a good REF match either,
+                    // at least not in low coverage.  So require solid coverage
+                    // before we start utilising such quals.
+                    if (b != 0)
+                        b = 5;
+                    q = (int)bam_get_qual(p->b)[p->qpos];
+                    seqQ = (3*seqQ + 2*q)/8;
+                }
+                if (_n > 20 && seqQ > 40) seqQ = 40;
+            }
+
+            // Note baseQ changes some output fields such as I16, but has no
+            // significant affect on "call".
             baseQ  = p->aux>>8&0xff;
         }
         else
@@ -478,9 +518,19 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
             for (i=0; i<4; i++) r->ADF[i] += lroundf((float)dp_ambig * r->ADF[i]/dp);
     }
 
+    // Else consider downgrading bca->bases[] scores by AD vs AD_ref_missed
+    // ratios.  This is detrimental on Illumina, but beneficial on PacBio CCS.
+    // It's possibly related to the homopolyer error likelihoods or overall
+    // Indel accuracy.  Maybe tie this in to the -h option?
+
     r->ori_depth = ori_depth;
     // glfgen
     errmod_cal(bca->e, n, 5, bca->bases, r->p); // calculate PL of each genotype
+
+    // TODO: account for the number of unassigned reads.  If depth is 50,
+    // but AD is 5,7 then it may look like a variant but it probably
+    // should be low quality.
+
     return n;
 }
 

bam2bcf.h

@@ -122,14 +122,15 @@ typedef struct __bcf_callaux_t {
     // for internal uses
     int max_bases;
     int indel_types[4];     // indel lengths
-    int indel_win_size, indels_v20;
+    int indel_win_size, indels_v20, edlib;
     int maxins, indelreg;
     int read_len;
     char *inscns;
     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 del_bias;    // (-.9 < x < .9) error profile; >0 => more del, <0 => more ins
     int32_t *ref_nm, *alt_nm;   // pointers to bcf_call_t.{ref_nm,alt_nm}
     unsigned int nnm[2];        // number of nm observations
     float nm[2];                // cumulative count of mismatches in ref and alt reads
@@ -193,11 +194,35 @@ extern "C" {
                      const bcf_callaux_t *bca, const char *ref);
     int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref);
     int bcf_iaux_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref);
+    int bcf_edlib_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos,
+                           bcf_callaux_t *bca, const char *ref, int ref_len);
+
     void bcf_callaux_clean(bcf_callaux_t *bca, bcf_call_t *call);
 
     int bcf_cgp_l_run(const char *ref, int pos);
     int est_indelreg(int pos, const char *ref, int l, char *ins4);
 
+/* ----------------------------------------------------------------------
+ * Shared between bam2bcf_indel.c and bam2bcf_edlib.c
+ */
+
+// Take a reference position tpos and convert to a query position (returned).
+// This uses the CIGAR string plus alignment c->pos to do the mapping.
+//
+// *_tpos is returned as tpos if query overlaps tpos, but for deletions
+// it'll be either the start (is_left) or end (!is_left) ref position.
+int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos);
+
+// Identify spft-clip length, position in seq, and clipped seq len
+void get_pos(const bcf_callaux_t *bca, bam_pileup1_t *p,
+             int *sc_len_r, int *slen_r, int *epos_r, int *end);
+
+// Compute the consensus for this sample 's', minus indels which
+// get added later.
+char *bcf_cgp_calc_cons(int n, int *n_plp, bam_pileup1_t **plp,
+                        int pos, int *types, int n_types,
+                        int max_ins, int s);
+
 #ifdef __cplusplus
 }
 #endif