Add new --ar, --ambig-reads option

which allows to retain the original behavior (the default) or
increment AD counters proportionally (--ar incAD) or just the
reference counter (--ar incAD0).

Note that the reference counter is not incremented explicitly by the
incAD0 mode, for now we are assuming that it is always the case that
the competing indel allele is the reference allele. When this is not
the case, the correct counter will be incremented, only the 'incAD0'
is somewhat misleading in such (rare) cases.

Author: pd3

bam2bcf.c

@@ -205,6 +205,8 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
 
     // fill the bases array
     double nqual_over_60 = bca->nqual / 60.0;
+    int ADR_ref_missed[4] = {0};
+    int ADF_ref_missed[4] = {0};
     for (i = n = 0; i < _n; ++i) {
         const bam_pileup1_t *p = pl + i;
         int q, b, mapQ, baseQ, is_diff, min_dist, seqQ;
@@ -218,12 +220,12 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
             seqQ = q = (p->aux & 0xff); // mp2 + builtin indel-bias
             if (q < bca->min_baseQ)
             {
-                if (!p->indel && r->ADF)
+                if (!p->indel && b < 4)
                 {
-                    if ( bam_is_rev(p->b) )
-                        r->ADR[b]++;
+                    if (bam_is_rev(p->b))
+                        ADR_ref_missed[b]++;
                     else
-                        r->ADF[b]++;
+                        ADF_ref_missed[b]++;
                 }
                 continue;
             }
@@ -331,6 +333,30 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t
             bca->alt_scl[sc_len]++;
         }
     }
+
+    // Compensate for AD not being counted on low quality REF indel matches.
+    if ( r->ADF && bca->ambig_reads==B2B_INC_AD0 )
+    {
+        for (i=0; i<4; i++) // verify: are the counters ever non-zero for i!=0?
+        {
+            r->ADR[i] += ADR_ref_missed[i];
+            r->ADF[i] += ADF_ref_missed[i];
+        }
+    }
+    else if ( r->ADF && bca->ambig_reads==B2B_INC_AD )
+    {
+        int dp = 0, dp_ambig = 0;
+        for (i=0; i<4; i++) dp += r->ADR[i];
+        for (i=0; i<4; i++) dp_ambig += ADR_ref_missed[i];
+        if ( dp )
+            for (i=0; i<4; i++) r->ADR[i] += lroundf((float)dp_ambig * r->ADR[i]/dp);
+        dp = 0, dp_ambig = 0;
+        for (i=0; i<4; i++) dp += r->ADF[i];
+        for (i=0; i<4; i++) dp_ambig += ADF_ref_missed[i];
+        if ( dp )
+            for (i=0; i<4; i++) r->ADF[i] += lroundf((float)dp_ambig * r->ADF[i]/dp);
+    }
+
     r->ori_depth = ori_depth;
     // glfgen
     errmod_cal(bca->e, n, 5, bca->bases, r->p); // calculate PL of each genotype

bam2bcf.h

@@ -65,6 +65,10 @@ DEALINGS IN THE SOFTWARE.  */
 
 #define B2B_MAX_ALLELES 5
 
+#define B2B_DROP      0
+#define B2B_INC_AD    1
+#define B2B_INC_AD0   2
+
 #define PLP_HAS_SOFT_CLIP(i) ((i)&1)
 #define PLP_HAS_INDEL(i)     ((i)&2)
 #define PLP_SAMPLE_ID(i)     ((i)>>2)
@@ -74,7 +78,7 @@ DEALINGS IN THE SOFTWARE.  */
 #define PLP_SET_SAMPLE_ID(i,n)   ((i)|=(n)<<2)
 
 typedef struct __bcf_callaux_t {
-    int fmt_flag;
+    int fmt_flag, ambig_reads;
     int capQ, min_baseQ, max_baseQ, delta_baseQ;
     int openQ, extQ, tandemQ; // for indels
     uint32_t min_support, max_support; // for collecting indel candidates

doc/bcftools.txt

@@ -1975,6 +1975,14 @@ INFO/DPR    .. Deprecated in favor of INFO/AD; Number of high-quality bases for
         ont	           -B -Q5 --max-BQ 30 -I
         pacbio-ccs     -D -Q5 --max-BQ 50 -F0.1 -o25 -e1 -M99999
 
+*--ambiguous-reads* 'drop'|'incAD'|'incAD0'::
+    What to do with ambiguous indel reads that do not span an entire
+    short tandem repeat region: discard ambiguous reads from calling
+    and do not increment high-quality AD depth counters ('drop'),
+    exclude from calling but increment AD counters proportionally ('incAD'),
+    exclude from calling and increment the first value of the AD counter
+    ('incAD0') ['drop']
+
 *-e, --ext-prob* 'INT'::
     Phred-scaled gap extension sequencing error probability. Reducing 'INT'
     leads to longer indels [20]

mpileup.c

@@ -68,7 +68,7 @@ typedef struct _mplp_pileup_t mplp_pileup_t;
 // Data shared by all bam files
 typedef struct {
     int min_mq, flag, min_baseQ, max_baseQ, delta_baseQ, capQ_thres, max_depth,
-        max_indel_depth, max_read_len, fmt_flag;
+        max_indel_depth, max_read_len, fmt_flag, ambig_reads;
     int rflag_require, rflag_filter, output_type;
     int openQ, extQ, tandemQ, min_support; // for indels
     double min_frac; // for indels
@@ -842,6 +842,7 @@ static int mpileup(mplp_conf_t *conf)
     conf->bca->min_support = conf->min_support;
     conf->bca->per_sample_flt = conf->flag & MPLP_PER_SAMPLE;
     conf->bca->fmt_flag = conf->fmt_flag;
+    conf->bca->ambig_reads = conf->ambig_reads;
 
     conf->bc.bcf_hdr = conf->bcf_hdr;
     conf->bc.n  = nsmpl;
@@ -1158,11 +1159,11 @@ static void print_usage(FILE *fp, const mplp_conf_t *mplp)
         "  -o, --open-prob INT     Phred-scaled gap open seq error probability [%d]\n", mplp->openQ);
     fprintf(fp,
         "  -p, --per-sample-mF     apply -m and -F per-sample for increased sensitivity\n"
-        "  -P, --platforms STR     comma separated list of platforms for indels [all]\n");
+        "  -P, --platforms STR     comma separated list of platforms for indels [all]\n"
+        "  --ar, --ambig-reads STR   What to do with ambiguous indel reads: drop,incAD,incAD0 [drop]\n");
     fprintf(fp,
-        "      --indel-bias FLOAT  Raise to favour recall over precision [%.2f]\n"
-        "      --ambiguous-reads FLAGS  Count ambiguous indel reads as: alt,ref,drop [drop]\n"
-        "\n", mplp->indel_bias);
+        "      --indel-bias FLOAT  Raise to favour recall over precision [%.2f]\n", mplp->indel_bias);
+    fprintf(fp,"\n");
     fprintf(fp,
         "Configuration profiles activated with -X, --config:\n"
         "    1.12:        -Q13 -h100 -m1 -F0.002\n"
@@ -1208,7 +1209,7 @@ int main_mpileup(int argc, char *argv[])
     // the default to be changed in future, see also parse_format_flag()
     mplp.fmt_flag = B2B_INFO_VDB|B2B_INFO_RPB|B2B_INFO_SCB|B2B_INFO_ZSCORE;
     mplp.max_read_len = 500;
-
+    mplp.ambig_reads = B2B_DROP;
     hts_srand48(0);
 
     static const struct option lopts[] =
@@ -1270,6 +1271,8 @@ int main_mpileup(int argc, char *argv[])
         {"config", required_argument, NULL, 'X'},
         {"mwu-u", no_argument, NULL, 'U'},
         {"seed", required_argument, NULL, 13},
+        {"ambig-reads", required_argument, NULL, 14},
+        {"ar", required_argument, NULL, 14},
         {NULL, 0, NULL, 0}
     };
     while ((c = getopt_long(argc, argv, "Ag:f:r:R:q:Q:C:BDd:L:b:P:po:e:h:Im:F:EG:6O:xa:s:S:t:T:M:X:U",lopts,NULL)) >= 0) {
@@ -1410,6 +1413,12 @@ int main_mpileup(int argc, char *argv[])
             }
             break;
         case 13: hts_srand48(atoi(optarg)); break;
+        case 14:
+            if ( !strcasecmp(optarg,"drop") ) mplp.ambig_reads = B2B_DROP;
+            else if ( !strcasecmp(optarg,"incAD") ) mplp.ambig_reads = B2B_INC_AD;
+            else if ( !strcasecmp(optarg,"incAD0") ) mplp.ambig_reads = B2B_INC_AD0;
+            else error("The option to --ambig-reads not recognised: %s\n",optarg);
+            break;
         default:
             fprintf(stderr,"Invalid option: '%c'\n", c);
             return 1;

test/mpileup/indel-AD.1.out

@@ -166,9 +166,9 @@
 000000F	535	.	G	A,<*>	0	.	DP=125;I16=65,52,0,1,4309,171791,12,144,7020,421200,60,3600,2679,64385,25,625;QS=0.997223,0.00277713,0;SGB=-0.379885;RPBZ=-1.14518;BQBZ=-1.74874;SCBZ=0.0762539;FS=0;MQ0F=0	PL:AD	0,255,255,255,255,255:117,1,0
 000000F	536	.	T	G,A,<*>	0	.	DP=125;I16=65,51,0,2,4274,171298,24,288,6960,417600,120,7200,2661,64041,48,1154;QS=0.994416,0.002792,0.002792,0;VDB=0.1;SGB=-0.453602;RPBZ=-1.69957;BQBZ=-2.39373;SCBZ=-0.714873;FS=0;MQ0F=0	PL:AD	0,255,255,255,255,255,255,255,255,255:116,1,1,0
 000000F	537	.	A	<*>	0	.	DP=125;I16=65,53,0,0,4390,175290,0,0,7080,424800,0,0,2713,65375,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,255,255:118,0
-000000F	537	.	AC	A	0	.	INDEL;IDV=60;IMF=0.48;DP=125;I16=38,27,30,27,2600,104000,2280,91200,3900,234000,3420,205200,1473,35311,1298,31150;QS=0.267262,0.732738;VDB=0.0291871;SGB=-0.693147;RPBZ=-0.543708;SCBZ=0.641853;FS=0;MQ0F=0	PL:AD	255,0,26:65,57
+000000F	537	.	AC	A	0	.	INDEL;IDV=60;IMF=0.48;DP=125;I16=37,25,30,27,2480,99200,2280,91200,3720,223200,3420,205200,1399,33485,1298,31150;QS=0.260049,0.739951;VDB=0.0363603;SGB=-0.693147;RPBZ=-0.543708;SCBZ=0.641853;FS=0;MQ0F=0	PL:AD	255,0,27:62,57
 000000F	538	.	C	<*>	0	.	DP=65;I16=36,26,0,0,2195,86349,0,0,3720,223200,0,0,1432,34600,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,187,255:62,0
-000000F	538	.	CT	C	0	.	INDEL;IDV=64;IMF=0.512;DP=125;I16=32,29,36,25,2440,97600,2440,97600,3660,219600,3660,219600,1385,33095,1380,33214;QS=0.227073,0.772927;VDB=0.0341788;SGB=-0.693147;RPBZ=0.242079;SCBZ=-0.994262;FS=0;MQ0F=0	PL:AD	255,0,26:61,61
+000000F	538	.	CT	C	0	.	INDEL;IDV=64;IMF=0.512;DP=125;I16=30,26,36,25,2240,89600,2440,97600,3360,201600,3660,219600,1279,30735,1380,33214;QS=0.218762,0.781238;VDB=0.0165787;SGB=-0.693147;RPBZ=0.242079;SCBZ=-0.994262;FS=0;MQ0F=0	PL:AD	255,0,27:56,61
 000000F	539	.	T	<*>	0	.	DP=60;I16=29,26,0,0,2120,86238,0,0,3300,198000,0,0,1260,30374,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,166,255:55,0
 000000F	540	.	G	<*>	0	.	DP=124;I16=64,53,0,0,4130,161310,0,0,7020,421200,0,0,2703,65511,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,255,255:117,0
 000000F	541	.	G	<*>	0	.	DP=124;I16=64,53,0,0,4143,160525,0,0,7020,421200,0,0,2705,65703,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,255,255:117,0
@@ -286,11 +286,11 @@
 000000F	653	.	A	<*>	0	.	DP=21;I16=9,12,0,0,804,31130,0,0,1260,75600,0,0,303,5555,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,63,255:21,0
 000000F	654	.	A	<*>	0	.	DP=21;I16=9,12,0,0,659,22061,0,0,1260,75600,0,0,285,5117,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,63,255:21,0
 000000F	655	.	C	<*>	0	.	DP=21;I16=9,12,0,0,664,22342,0,0,1260,75600,0,0,266,4666,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,63,255:21,0
-000000F	655	.	CACAATACAA	CACAA	0	.	INDEL;IDV=6;IMF=0.285714;DP=21;I16=4,11,5,1,1800,216000,720,86400,900,54000,360,21600,166,2878,100,1788;QS=0.371728,0.628272;VDB=2.08672e-10;SGB=-0.616816;RPBZ=-2.81289;SCBZ=-2.52262;FS=0;MQ0F=0	PL:AD	129,0,28:15,6
+000000F	655	.	CACAATACAA	CACAA	0	.	INDEL;IDV=6;IMF=0.285714;DP=21;I16=0,2,5,1,240,28800,720,86400,120,7200,360,21600,46,1060,100,1788;QS=0.0977444,0.902256;VDB=0.00018837;SGB=-0.616816;RPBZ=-2.81289;SCBZ=-2.52262;FS=0;MQ0F=0	PL:AD	159,0,2:2,6
 000000F	656	.	A	<*>	0	.	DP=11;I16=4,7,0,0,404,15690,0,0,660,39600,0,0,141,2411,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,33,255:11,0
 000000F	657	.	C	<*>	0	.	DP=11;I16=4,7,0,0,413,15607,0,0,660,39600,0,0,131,2189,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,33,255:11,0
 000000F	658	.	A	<*>	0	.	DP=10;I16=3,7,0,0,121,1651,0,0,600,36000,0,0,122,1986,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,30,79:10,0
-000000F	658	.	AA	AAATTA	0	.	INDEL;IDV=2;IMF=0.125;DP=16;I16=7,6,1,2,1300,130000,300,30000,780,46800,180,10800,167,2727,50,902;QS=0.439462,0.560538;VDB=1.46535e-07;SGB=-0.511536;RPBZ=-1.70026;SCBZ=-1.35678;FS=0;MQ0F=0	PL:AD	70,0,24:13,3
+000000F	658	.	AA	AAATTA	0	.	INDEL;IDV=2;IMF=0.125;DP=16;I16=2,1,1,2,300,30000,300,30000,180,10800,180,10800,59,1205,50,902;QS=0.31694,0.68306;VDB=0.00155977;SGB=-0.511536;RPBZ=-1.70026;SCBZ=-1.35678;FS=0;MQ0F=0	PL:AD	99,0,38:3,3
 000000F	659	.	A	<*>	0	.	DP=10;I16=3,5,0,0,86,1088,0,0,480,28800,0,0,75,1077,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,24,63:8,0
 000000F	660	.	T	<*>	0	.	DP=2;I16=0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;QS=0,0;FS=0;MQ0F=0	PL:AD	0,0,0:0,0
 000000F	661	.	A	<*>	0	.	DP=8;I16=0,2,0,0,8,32,0,0,120,7200,0,0,26,340,0,0;QS=1,0;FS=0;MQ0F=0	PL:AD	0,6,7:2,0

test/mpileup/indel-AD.2.out

@@ -21,4 +21,4 @@
 ##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (high-quality bases)">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
 11	75	.	G	<*>	0	.	DP=68;I16=6,62,0,0,2437,87909,0,0,3770,217210,0,0,838,15940,0,0;QS=1,0;MQSBZ=0.140975;FS=0;MQ0F=0	PL:AD	0,205,255:68,0
-11	75	.	GTAAAATAAAATAAAATAAAATAAA	GTAAAATAAAATAAAATAAAATAAAATAAA	0	.	INDEL;IDV=6;IMF=0.0882353;DP=68;I16=5,57,1,5,7440,892800,720,86400,3596,212164,174,5046,691,12331,147,3609;QS=0.834915,0.165085;VDB=0.13856;SGB=-0.616816;RPBZ=-3.24592;MQBZ=-6.13241;MQSBZ=0.140975;SCBZ=-0.546919;FS=0;MQ0F=0	PL:AD	0,54,150:62,6
+11	75	.	GTAAAATAAAATAAAATAAAATAAA	GTAAAATAAAATAAAATAAAATAAAATAAA	0	.	INDEL;IDV=6;IMF=0.0882353;DP=68;I16=5,9,1,5,1680,201600,720,86400,840,50400,174,5046,244,5778,147,3609;QS=0.730233,0.269767;VDB=0.00674908;SGB=-0.616816;RPBZ=-3.24592;MQBZ=-6.13241;MQSBZ=0.140975;SCBZ=-0.546919;FS=0;MQ0F=0	PL:AD	83,0,244:14,6

test/mpileup/indel-AD.3.out

@@ -0,0 +1,24 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=11,length=150>
+##ALT=<ID=*,Description="Represents allele(s) other than observed.">
+##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
+##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of raw reads supporting an indel">
+##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of raw reads supporting an indel">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
+##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
+##INFO=<ID=RPBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Read Position Bias (closer to 0 is better)">
+##INFO=<ID=MQBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Mapping Quality Bias (closer to 0 is better)">
+##INFO=<ID=BQBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Base Quality Bias (closer to 0 is better)">
+##INFO=<ID=MQSBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Mapping Quality vs Strand Bias (closer to 0 is better)">
+##INFO=<ID=SCBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Soft-Clip Length Bias (closer to 0 is better)">
+##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
+##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
+##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
+##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
+##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (high-quality bases)">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
+11	75	.	G	<*>	0	.	DP=68;I16=6,62,0,0,2437,87909,0,0,3770,217210,0,0,838,15940,0,0;QS=1,0;MQSBZ=0.140975;FS=0;MQ0F=0	PL:AD	0,205,255:68,0
+11	75	.	GTAAAATAAAATAAAATAAAATAAA	GTAAAATAAAATAAAATAAAATAAAATAAA	0	.	INDEL;IDV=6;IMF=0.0882353;DP=68;I16=5,9,1,5,1680,201600,720,86400,840,50400,174,5046,244,5778,147,3609;QS=0.730233,0.269767;VDB=0.00674908;SGB=-0.616816;RPBZ=-3.24592;MQBZ=-6.13241;MQSBZ=0.140975;SCBZ=-0.546919;FS=0;MQ0F=0	PL:AD	83,0,244:45,23

test/mpileup/indel-AD.4.out

@@ -0,0 +1,24 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=11,length=150>
+##ALT=<ID=*,Description="Represents allele(s) other than observed.">
+##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
+##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of raw reads supporting an indel">
+##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of raw reads supporting an indel">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
+##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
+##INFO=<ID=RPBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Read Position Bias (closer to 0 is better)">
+##INFO=<ID=MQBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Mapping Quality Bias (closer to 0 is better)">
+##INFO=<ID=BQBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Base Quality Bias (closer to 0 is better)">
+##INFO=<ID=MQSBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Mapping Quality vs Strand Bias (closer to 0 is better)">
+##INFO=<ID=SCBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Soft-Clip Length Bias (closer to 0 is better)">
+##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
+##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
+##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
+##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
+##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths (high-quality bases)">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample
+11	75	.	G	<*>	0	.	DP=68;I16=6,62,0,0,2437,87909,0,0,3770,217210,0,0,838,15940,0,0;QS=1,0;MQSBZ=0.140975;FS=0;MQ0F=0	PL:AD	0,205,255:68,0
+11	75	.	GTAAAATAAAATAAAATAAAATAAA	GTAAAATAAAATAAAATAAAATAAAATAAA	0	.	INDEL;IDV=6;IMF=0.0882353;DP=68;I16=5,9,1,5,1680,201600,720,86400,840,50400,174,5046,244,5778,147,3609;QS=0.730233,0.269767;VDB=0.00674908;SGB=-0.616816;RPBZ=-3.24592;MQBZ=-6.13241;MQSBZ=0.140975;SCBZ=-0.546919;FS=0;MQ0F=0	PL:AD	83,0,244:62,6

test/test.pl

@@ -681,6 +681,8 @@
 test_mpileup($opts,in=>[qw(mpileup.3 mpileup.4)],out=>'mpileup/mpileup.11.out',args=>q[-G {PATH}/mplp.11.rgs]);
 test_mpileup($opts,in=>[qw(indel-AD.1)],out=>'mpileup/indel-AD.1.out',ref=>'indel-AD.1.fa',args=>q[-a AD]);
 test_mpileup($opts,in=>[qw(indel-AD.2)],out=>'mpileup/indel-AD.2.out',ref=>'indel-AD.2.fa',args=>q[-a AD -r 11:75]);
+test_mpileup($opts,in=>[qw(indel-AD.2)],out=>'mpileup/indel-AD.3.out',ref=>'indel-AD.2.fa',args=>q[-a AD -r 11:75 --ambig-reads incAD]);
+test_mpileup($opts,in=>[qw(indel-AD.2)],out=>'mpileup/indel-AD.4.out',ref=>'indel-AD.2.fa',args=>q[-a AD -r 11:75 --ambig-reads incAD0]);
 test_mpileup($opts,in=>[qw(mpileup-SCR)],out=>'mpileup/mpileup-SCR.out',ref=>'mpileup-SCR.fa',args=>q[-a INFO/SCR,FMT/SCR]);
 test_csq($opts,in=>'csq',out=>'csq.1.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.gff3');
 test_csq($opts,in=>'csq',out=>'csq.1.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.chr.gff3');