Repository: samtools/samtools (c)

Issue:
Minimum depth threshold for "samtools coverage"?
#### Is your feature request related to a problem? Please specify.

I'd like a simple way to calculate the coverage stats for a bam file. The closest out-of-the-box thing I can find is `samtools coverage`. I run the command like so:

```
samtools coverage input.bam
```

The simple summary stats are great, and what I need, _except_ for the fact that the 'covbases' statistic is defined as 'Number of covered bases with depth >= 1'. Instead, I would prefer to have the statistic of 'Number of covered bases with depth >= X', where X is a minimum depth defined by a command line flag.

#### Describe the solution you would like.
I would like to be able to use a command like:

```
samtools coverage --min-depth 20 input.bam
```
With the output statistics being calculated (in this case) only for bases covered with depth >=20. Would this be possible? So far I use a different command like:

```
bedtools genomecov -ibam input.bam -d | awk '$3 >= 20 { count++ } END { print (count/NR)*100 }'
```

But I'd prefer to have the much simpler `samtools` command and simplify the dependencies of my pipeline.

Would this be possible? Thanks!




Patch:
diff --git a/coverage.c b/coverage.c
index 92f497dc8..7fcb34f7e 100644
--- a/coverage.c
+++ b/coverage.c
@@ -1,7 +1,7 @@
 /* coverage.c -- samtools coverage subcommand
 
     Copyright (C) 2018,2019 Florian Breitwieser
-    Portions copyright (C) 2019-2021, 2023-2024 Genome Research Ltd.
+    Portions copyright (C) 2019-2021, 2023-2025 Genome Research Ltd.
 
     Author: Florian P Breitwieser <florian.bw@gmail.com>
 
@@ -119,6 +119,8 @@ static int usage(void) {
             "  -d, --depth INT         maximum allowed coverage depth [1000000].\n"
             "                          If 0, depth is set to the maximum integer value,\n"
             "                          effectively removing any depth limit.\n"
+            "      --min-depth INT     minimum coverage depth below which a position \n"
+            "                          to be ignored [1]\n"
             "Output options:\n"
             "  -m, --histogram         show histogram instead of tabular output\n"
             "  -D, --plot-depth        plot depth instead of tabular output\n"
@@ -317,7 +319,7 @@ int main_coverage(int argc, char *argv[]) {
     int fail_flags = (BAM_FUNMAP | BAM_FSECONDARY | BAM_FQCFAIL | BAM_FDUP); // Default fail flags
     int required_flags = 0;
     int print_value_warning = 0;
-
+    int mindepth = 1;   //min depth in file, below which pos is ignored
     int *n_plp = NULL;
     sam_hdr_t *h = NULL; // BAM header of the 1st input
 
@@ -351,6 +353,7 @@ int main_coverage(int argc, char *argv[]) {
         {"region", required_argument, NULL, 'r'},
         {"help", no_argument, NULL, 'h'},
         {"depth", required_argument, NULL, 'd'},
+        {"min-depth", required_argument, NULL, 3},
         { NULL, 0, NULL, 0 }
     };
 
@@ -367,6 +370,11 @@ int main_coverage(int argc, char *argv[]) {
                 if ((fail_flags = bam_str2flag(optarg)) < 0) {
                     fprintf(stderr,"Could not parse --ff %s\n", optarg); return EXIT_FAILURE;
                 }; break;
+            case 3: //min-depth
+                if ((i = atoi(optarg)) > 0) {
+                    mindepth = i;
+                }
+                break;
             case 'o': opt_output_file = optarg; opt_full_width = false; break;
             case 'l': opt_min_len = atoi(optarg); break;
             case 'q': opt_min_mapQ = atoi(optarg); break;
@@ -607,6 +615,7 @@ int main_coverage(int argc, char *argv[]) {
         }
 
         bool count_base = false;
+        unsigned long long summed_baseQ = 0, quality_bases = 0, depth = 0;
         for (i = 0; i < n_bam_files; ++i) { // base level filters have to go here
             int depth_at_pos = n_plp[i];
             for (j = 0; j < n_plp[i]; ++j) {
@@ -618,8 +627,8 @@ int main_coverage(int argc, char *argv[]) {
                     if (bam_get_qual(p->b)[p->qpos] < opt_min_baseQ) {
                         --depth_at_pos; // low base quality
                     } else {
-                        stats[tid].summed_baseQ += bam_get_qual(p->b)[p->qpos];
-                        stats[tid].quality_bases++;
+                        summed_baseQ += bam_get_qual(p->b)[p->qpos];
+                        ++quality_bases;
                     }
                 } else {
                     print_value_warning = 1; // no quality at position
@@ -628,14 +637,18 @@ int main_coverage(int argc, char *argv[]) {
 
             if (depth_at_pos > 0) {
                 count_base = true;
-                stats[tid].summed_coverage += depth_at_pos;
+                depth += depth_at_pos;
             }
 
             if(current_bin < n_bins && opt_plot_coverage) {
                 hist[current_bin] += depth_at_pos;
             }
         }
-        if (count_base) {
+        if (count_base && depth >= mindepth) {
+            stats[tid].summed_coverage += depth;
+            stats[tid].summed_baseQ += summed_baseQ;
+            stats[tid].quality_bases += quality_bases;
+
             stats[tid].n_covered_bases++;
             if (opt_print_histogram && current_bin < n_bins && !opt_plot_coverage)
                 ++(hist[current_bin]); // Histogram based on breadth of coverage
diff --git a/doc/samtools-coverage.1 b/doc/samtools-coverage.1
index 41cfc80b6..55b16a1cc 100644
--- a/doc/samtools-coverage.1
+++ b/doc/samtools-coverage.1
@@ -3,7 +3,7 @@
 .SH NAME
 samtools coverage \- produces a histogram or table of coverage per chromosome
 .\"
-.\" Copyright (C) 2019, 2021, 2023 Genome Research Ltd.
+.\" Copyright (C) 2019, 2021, 2023, 2025 Genome Research Ltd.
 .\"
 .\" Author: James Bonfield <jkb@sanger.ac.uk>
 .\"
@@ -95,6 +95,9 @@ Filter flags: skip reads with mask bits set
 .BI -d,\ --depth \ INT
 Maximum allowed coverage depth [1000000]. If 0, depth is set to the maximum
 integer value effectively removing any depth limit.
+.TP
+.BI --min-depth\  INT
+Minimum coverage depth, below which a position is ignored [1]
 
 .PP
 Output options:


Test patch:
diff --git a/test/coverage/1.expected b/test/coverage/1.expected
new file mode 100644
index 000000000..35cd07c2a
--- /dev/null
+++ b/test/coverage/1.expected
@@ -0,0 +1,3 @@
+#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
+T1	1	40	6	19	47.5	0.725	7.62	36.7
+T2	1	40	6	8	20	0.4	7.12	49
diff --git a/test/coverage/2.expected b/test/coverage/2.expected
new file mode 100644
index 000000000..b4b909dc3
--- /dev/null
+++ b/test/coverage/2.expected
@@ -0,0 +1,3 @@
+#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
+T1	1	40	6	10	25	0.5	7.5	36.7
+T2	1	40	6	6	15	0.35	7.43	49
diff --git a/test/coverage/3.expected b/test/coverage/3.expected
new file mode 100644
index 000000000..4165972a6
--- /dev/null
+++ b/test/coverage/3.expected
@@ -0,0 +1,3 @@
+#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
+T2	1	40	6	1	2.5	0.05	8	49
+T1	1	40	0	0	0	0	0	0
diff --git a/test/coverage/4.expected b/test/coverage/4.expected
new file mode 100644
index 000000000..cd66a4283
--- /dev/null
+++ b/test/coverage/4.expected
@@ -0,0 +1,3 @@
+#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
+T1	1	40	11	19	47.5	1.3	7.46	36.8
+T2	1	40	12	8	20	0.8	7.12	49
diff --git a/test/coverage/5.expected b/test/coverage/5.expected
new file mode 100644
index 000000000..adf91a007
--- /dev/null
+++ b/test/coverage/5.expected
@@ -0,0 +1,3 @@
+#rname	startpos	endpos	numreads	covbases	coverage	meandepth	meanbaseq	meanmapq
+T1	1	40	11	4	10	0.4	5.25	36.8
+T2	1	40	12	6	15	0.7	7.43	49
diff --git a/test/dat/sample.sam b/test/dat/sample.sam
new file mode 100644
index 000000000..0eaa44af7
--- /dev/null
+++ b/test/dat/sample.sam
@@ -0,0 +1,30 @@
+@HD	VN:1.17	SO:coordinate
+@SQ	SN:T1	LN:40
+@SQ	SN:T2	LN:40
+@CO	@SQ SN* LN* AH AN AS DS  M5 SP TP UR
+@CO	@RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM
+@CO	@PG ID* PN CL PP DS VN
+@CO	this is a dummy alignment file to demonstrate different abilities of hts apis
+@CO	QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]…
+@CO	1234567890123456789012345678901234567890
+@CO	AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT	T1
+@CO	TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT	T2
+@CO	ITR1-ITR2M, ITR2-ITR2M are proper pairs in T1 and T2, UNMP1 is partly mapped and pair is unmapped, UNMP2 & 3 are unmapped
+@CO	A1-A2, A4-A3 are proper pairs with A4-A3 in different read order. A5 is secondary alignment
+@PG	ID:samtools	PN:samtools	VN:1.22-6-g3036eb9	CL:./samtools sort -o /tmp/s.sam ../htslib/samples/sample.sam
+ITR1	99	T1	5	40	4M	=	33	10	ACTG	()()
+UNMP1	73	T1	21	40	3M	*	0	5	GGG	&&1
+A1	99	T1	25	35	6M	=	31	8	ACTGTT	******
+A5	355	T1	25	55	4M	=	33	5	ACTG	PPPP
+B1	99	T1	25	35	6M	=	31	8	GCTATT	******
+B5	355	T1	25	35	4M	=	33	5	AGTG	PPPP
+A2	147	T1	31	33	6M	=	25	-8	ACTGTT	()()()
+ITR1M	147	T1	33	37	4M	=	5	-10	ACTG	$$$$
+A4	99	T2	12	50	3M	=	23	5	GAA	()(
+B4	99	T2	12	50	3M	=	23	5	GAT	()(
+ITR2	147	T2	23	49	2M	=	35	-10	TT	**
+A3	147	T2	23	47	2M1X	=	12	-5	TTG	(((
+B3	147	T2	23	47	2M1X	=	12	-5	TAG	(((
+ITR2M	99	T2	35	51	2M	=	23	10	AA	&&
+UNMP2	141	*	0	0	*	*	0	7	AA	&&
+UNMP3	77	*	0	0	*	*	0	5	GGG	&&2
diff --git a/test/test.pl b/test/test.pl
index d183c7e9d..eb06ec99c 100755
--- a/test/test.pl
+++ b/test/test.pl
@@ -79,6 +79,7 @@
 test_reset($opts);
 test_checksum($opts);
 test_checksum($opts, threads=>2);
+test_coverage($opts);
 
 print "\nNumber of tests:\n";
 printf "    total            .. %d\n", $$opts{nok}+$$opts{nfailed}+$$opts{nxfail}+$$opts{nxpass};
@@ -3984,3 +3985,22 @@ sub test_checksum
     }
     test_cmd($opts, out=>"checksum/chk2.2.expected", cmd=>"$$opts{bin}/samtools $chk -m $$opts{path}/checksum/chk2-*.tmp.chk | sed 's/\\(# Checksum[^:]*:\\).*/\\1/'");
 }
+
+
+sub test_coverage
+{
+    my ($opts, %args) = @_;
+
+    #basic / existing
+    test_cmd($opts, out=>"coverage/1.expected", cmd=>"$$opts{bin}/samtools coverage $$opts{path}/dat/sample.sam");
+    #coverage --min-depth 1
+    test_cmd($opts, out=>"coverage/1.expected", cmd=>"$$opts{bin}/samtools coverage --min-depth 1 $$opts{path}/dat/sample.sam");
+    #coverage --min-depth 2
+    test_cmd($opts, out=>"coverage/2.expected", cmd=>"$$opts{bin}/samtools coverage --min-depth 2 $$opts{path}/dat/sample.sam");
+    #coverage --min-depth 2 -Q 8 -q 45
+    test_cmd($opts, out=>"coverage/3.expected", cmd=>"$$opts{bin}/samtools coverage --min-depth 2 -Q 8 -q 45 $$opts{path}/dat/sample.sam");
+    #shows coverage is based on all inputs
+    cmd("cat '$$opts{path}/dat/sample.sam' | sed '/A1/d' > $$opts{tmp}/sample1.sam");
+    test_cmd($opts, out=>"coverage/4.expected", cmd=>"$$opts{bin}/samtools coverage --min-depth 1 $$opts{path}/dat/sample.sam $$opts{tmp}/sample1.sam");
+    test_cmd($opts, out=>"coverage/5.expected", cmd=>"$$opts{bin}/samtools coverage --min-depth 4 $$opts{path}/dat/sample.sam $$opts{tmp}/sample1.sam");
+}