Question: Dear all,

Today I'm trying to clarify my understanding of `--max-depth` argument from `bcftools mpileup` command, which is currently described as:

> At a position, read maximally INT reads per input file. Note that the original samtools mpileup command had a minimum value of 8000/n where n was the number of input files given to mpileup. This means that in samtools mpileup the default was highly likely to be increased and the -d parameter would have an effect only once above the cross-sample minimum of 8000. This behavior was problematic when working with a combination of single- and multi-sample bams, therefore in bcftools mpileup the user is given the full control (and responsibility), and an informative message is printed instead [250]

My question is: I get that this argument sets the maximum number of reads used to study each position per input file, but is there any minimum threshold (e.g. 100) below which user-specified value has no effects?

For example, I'm currently using bcftools' `mpileup` + `call` [pipeline](https://samtools.github.io/bcftools/howtos/variant-calling.html) for variant calling of a 30x WGS sample:

```
$ fasta=/Users/sbslee/OneDrive/TestFiles/genome.fa
$ bam=/Users/sbslee/OneDrive/PyPGx/GeT-RM/grch37-bam/NA18519_PyPGx.sorted.markdup.recal.bam
$ region=chr2:234663648-234663651
$ bcftools mpileup -Ou -f $fasta -r $region $bam | bcftools call -mv -Ov
```

Which gives:

```
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[W::hts_idx_load3] The index file is older than the data file: /Users/sbslee/OneDrive/PyPGx/GeT-RM/grch37-bam/NA18519_PyPGx.sorted.markdup.recal.bai
[mpileup] 1 samples in 1 input files
[mpileup] maximum number of reads per input file set to -d 250
##fileformat=VCFv4.2
##FILTER=
##bcftoolsVersion=1.15+htslib-1.15
##bcftoolsCommand=mpileup -Ou -f /Users/sbslee/OneDrive/TestFiles/genome.fa -r chr2:234663648-234663651 /Users/sbslee/OneDrive/PyPGx/GeT-RM/grch37-bam/NA18519_PyPGx.sorted.markdup.recal.bam
##reference=file:///Users/sbslee/OneDrive/TestFiles/genome.fa
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##ALT=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FORMAT=
##FORMAT=
##INFO=
##INFO=
##INFO=
##INFO=
##bcftools_callVersion=1.15+htslib-1.15
##bcftools_callCommand=call -mv -Ov; Date=Sat Apr 2 16:16:56 2022
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA18519_PyPGx
chr2	234663649	.	C	G	225.417	.	DP=35;VDB=0.241828;SGB=-0.693136;MQSBZ=0;FS=0;MQ0F=0;AC=2;AN=2;DP4=0,0,17,18;MQ=60	GT:PL	1/1:255,105,0
```

Because I didn't change `--max-depth`, `mpileup` used the default value of 250 and the sample had a total of 35 reads for `chr2-234663649-C-G` (i.e. `DP=35`).

Now, by setting `--max-depth 10` I would expect the variant to be still called but with `DP=10` instead of `DP=35`. Let's run this:

```
$ bcftools mpileup -Ou -f $fasta -r $region --max-depth 10 $bam | bcftools call -mv -Ov
```

Which gives:

```
Note: none of --samples-file, --ploidy or --ploidy-file given, assuming all sites are diploid
[W::hts_idx_load3] The index file is older than the data file: /Users/sbslee/OneDrive/PyPGx/GeT-RM/grch37-bam/NA18519_PyPGx.sorted.markdup.recal.bai
[mpileup] 1 samples in 1 input files
[mpileup] maximum number of reads per input file set to -d 10
Note: The maximum per-sample depth with -d 10 is 10.0x
##fileformat=VCFv4.2
##FILTER=
##bcftoolsVersion=1.15+htslib-1.15
##bcftoolsCommand=mpileup -Ou -f /Users/sbslee/OneDrive/TestFiles/genome.fa -r chr2:234663648-234663651 -d 10 /Users/sbslee/OneDrive/PyPGx/GeT-RM/grch37-bam/NA18519_PyPGx.sorted.markdup.recal.bam
##reference=file:///Users/sbslee/OneDrive/TestFiles/genome.fa
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##contig=
##ALT=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FORMAT=
##FORMAT=
##INFO=
##INFO=
##INFO=
##INFO=
##bcftools_callVersion=1.15+htslib-1.15
##bcftools_callCommand=call -mv -Ov; Date=Sat Apr 2 16:30:10 2022
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA18519_PyPGx
chr2	234663649	.	C	G	225.417	.	DP=33;VDB=0.184924;SGB=-0.693127;MQSBZ=0;FS=0;MQ0F=0;AC=2;AN=2;DP4=0,0,16,17;MQ=60	GT:PL	1/1:255,99,0
```

As you can see, the total depth decreased from 35 to 33 but it's nowhere close to our expected 10. Can someone please enlighten me what I'm missing here?

FYI, I'm using the latest version of bcftools available in bioconda (1.15):

```
$ conda list | grep "bcftools"
bcftools 1.15 hd1821b6_2 bioconda
```

P.S. This question has been cross-posted onto [bcftools GitHub repository issue](https://github.com/samtools/bcftools/issues/1694).

Answer: Sharing a satisfying answer I received from a contributor of bcftools:

> It's not reading ahead further in the file to work out whether a read which may currently cause over-depth will be needed later, so sadly this artefact does still happen, but I believe it may attempt to mitigate it to a certain degree by limiting the number of alignments starting within close proximity (which is the usual pattern for a high spike caused by misalignment of repeats). Hence the filtering isn't as strict as it may sound, an is perhaps closer to an average maximum depth threshold rather than something to take as a literal limit.

You can read the full discussion [here](https://github.com/samtools/bcftools/issues/1694).