From 66f5437dc01446d0981224e027fd3513b6c78184 Mon Sep 17 00:00:00 2001 From: Brad Langhorst Date: Sun, 16 Aug 2026 09:43:49 -0400 Subject: [PATCH 1/4] supports excluding reads by a tag's value, with tests --- functional-tests.sh | 29 +++++++++++++++++++++++++++++ mosdepth.nim | 31 +++++++++++++++++++++++++++++-- tests/dup-tags.bam | Bin 0 -> 2289 bytes tests/dup-tags.bam.bai | Bin 0 -> 784 bytes 4 files changed, 58 insertions(+), 2 deletions(-) create mode 100644 tests/dup-tags.bam create mode 100644 tests/dup-tags.bam.bai diff --git a/functional-tests.sh b/functional-tests.sh index b0fb606..831900c 100755 --- a/functional-tests.sh +++ b/functional-tests.sh @@ -76,6 +76,35 @@ run length_filter $exe t tests/ovl.bam --max-frag-len 79 assert_exit_code 0 assert_equal "MT 0 16569 0" "$(zgrep ^MT t.per-base.bed.gz)" +# --exclude-tag. tests/dup-tags.bam holds 3 reads all spanning MT:0-80: one untagged +# non-duplicate, one optical duplicate (DT:Z:SQ) and one library/PCR duplicate (DT:Z:LB), so the +# depth over that span is just the number of reads that survived filtering. + +# default excludes every duplicate: only the untagged read remains +run exclude_tag_default $exe t tests/dup-tags.bam +assert_exit_code 0 +assert_equal "MT 0 80 1 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" + +# keep duplicates (-F 772) but drop the optical one: untagged + library = 2 +run exclude_tag_optical $exe t tests/dup-tags.bam -F 772 --exclude-tag DT:SQ +assert_exit_code 0 +assert_equal "MT 0 80 2 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" + +# keep duplicates and exclude nothing: all 3 +run exclude_tag_none $exe t tests/dup-tags.bam -F 772 +assert_exit_code 0 +assert_equal "MT 0 80 3 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" + +# an absent tag must never exclude: filtering on a tag no read carries changes nothing +run exclude_tag_absent $exe t tests/dup-tags.bam -F 772 --exclude-tag ZZ:nope +assert_exit_code 0 +assert_equal "MT 0 80 3 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" + +# comma-separated list drops both classes, leaving only the untagged read +run exclude_tag_multi $exe t tests/dup-tags.bam -F 772 --exclude-tag DT:SQ,DT:LB +assert_exit_code 0 +assert_equal "MT 0 80 1 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" + run bad_frag_len_filter $exe t tests/ovl.bam --min-frag-len 10 --max-frag-len 9 assert_in_stderr "--max-frag-len was lower than --min-frag-len." assert_exit_code 2 diff --git a/mosdepth.nim b/mosdepth.nim index 5a35279..be229cc 100644 --- a/mosdepth.nim +++ b/mosdepth.nim @@ -251,7 +251,8 @@ proc coverage(bam: hts.Bam, arr: var coverage_t, region: var region_t, targets: seq[Target], mapq: int = -1, min_len: int = -1, max_len: int = int.high, eflag: uint16 = 1796, iflag: uint16 = 0, read_groups: seq[string] = (@[]), fast_mode: bool = false, - fragment_mode: bool = false, last_tid: var int = -1): int = + fragment_mode: bool = false, exclude_tags: seq[string] = (@[]), + last_tid: var int = -1): int = # depth updates arr in-place and yields the tid for each chrom. # returns -1 if the chrom is not found in the bam header # returns -2 if the chrom was found in the header, but there was no data for it @@ -261,6 +262,7 @@ proc coverage(bam: hts.Bam, arr: var coverage_t, region: var region_t, mate: Record seen = newTable[string, Record]() has_read_groups = read_groups.len > 0 + has_exclude_tags = exclude_tags.len > 0 var tid = if region != nil: get_tid(targets, region.chrom, last_tid) else: -1 if tid == -1: @@ -283,6 +285,20 @@ proc coverage(bam: hts.Bam, arr: var coverage_t, region: var region_t, var t = tag[string](rec, "RG") if t.isNone or not read_groups.contains(t.get): continue + # --exclude-tag TAG:VALUE. A read is dropped only when the tag is PRESENT and matches; + # an absent tag never excludes. That asymmetry is deliberate -- the motivating case is + # duplicate-type tags (DT:Z:SQ optical vs DT:Z:LB library), where only duplicates carry the + # tag at all and the untagged majority must be kept. + if has_exclude_tags: + var drop = false + for spec in exclude_tags: + let c = spec.find(':') + if c <= 0: continue + var t = tag[string](rec, spec[0 ..< c]) + if t.isSome and t.get == spec[(c+1) .. ^1]: + drop = true + break + if drop: continue if tgt.tid != rec.b.core.tid: raise newException(OSError, "expected only a single chromosome per query") @@ -633,6 +649,10 @@ proc main(bam: hts.Bam, chrom: region_t, mapq: int, min_len: int, max_len: int, if $args["--read-groups"] != "nil": for r in ($args["--read-groups"]).split(','): read_groups.add($r) + var exclude_tags: seq[string] = @[] + if $args["--exclude-tag"] != "nil": + for r in ($args["--exclude-tag"]).split(','): + exclude_tags.add($r) var levels = get_min_levels(targets) var chrom_region_distribution = newSeq[int64](region_distribution.len) @@ -706,7 +726,7 @@ proc main(bam: hts.Bam, chrom: region_t, mapq: int, min_len: int, max_len: int, rchrom = region_t(chrom: target.name) var tid = coverage(bam, arr, rchrom, targets, mapq, min_len, max_len, eflag, iflag, read_groups = read_groups, fast_mode = fast_mode, - fragment_mode = fragment_mode, + fragment_mode = fragment_mode, exclude_tags = exclude_tags, last_tid = last_tid) if tid == -1: continue # -1 means that chrom is not even in the bam if tid != -2: # -2 means there were no reads in the bam @@ -913,6 +933,10 @@ Other options: by ','. -m --use-median output median of each region (in --by) instead of mean. -R --read-groups only calculate depth for these comma-separated read groups IDs. + --exclude-tag exclude reads carrying this string tag with this exact value, + e.g. DT:SQ to drop optical duplicates while keeping PCR + duplicates. Reads lacking the tag are never excluded. May be + given as a comma-separated list (DT:SQ,XX:YY). -h --help show help """ @@ -979,6 +1003,9 @@ Other options: if $args["--read-groups"] != "nil": opts = opts or SamField.SAM_RGAUX.int + # an arbitrary aux tag cannot be read from CRAM unless the aux block is decoded + if $args["--exclude-tag"] != "nil": + opts = opts or SamField.SAM_AUX.int discard bam.set_option(FormatOption.CRAM_OPT_REQUIRED_FIELDS, opts) discard bam.set_option(FormatOption.CRAM_OPT_DECODE_MD, 0) diff --git a/tests/dup-tags.bam b/tests/dup-tags.bam new file mode 100644 index 0000000000000000000000000000000000000000..06ee81ad7b36755894e33449f1d0b93b2c442d78 GIT binary patch literal 2289 zcmV?Rv`X=}S{dU~e1s<*0ob~gqgAOXb~48#a4$$B9-pF}}|DDGPf z#z*mjzAb{X;ES#X3I2cGrp_^jfG+F0#isHd2^isdEa&|t9LVs>y>$8n5P-v{J&vR39i`tx(D9=&Bq!f-QQGqCvF6BVc z5;b#_G0WtRU7}Wox>PcY+m)`Jp$-FOYC9$BWT+*XaBaJa=4PmA3)_T%B`R~2aV8+W zstTIH!m@9cmp;a_BR+~z97?VN$Z zQdgs7K}ndfN(q)z+=0=n6r&|&yIL$!ID`RmQ=_bcGGPz2+>2PkW~P>*>6DbowOD59 z(1s0CC?-oQhdGuOOO%u;^;mX6x$wYflq@NWaxGS6C7EvJ7L;`K9-3xJNv=&&R8~^x zL*~VhDWm$l%#yMM)2PbIwzM8nyPz0UqGd=+N(w!eeMspgvvLiY9?KL3C4{NR;w5Eq zJr*q~Tk3U{D{r}nCssj~u`H>_qD7-zCJvwXVw6SaB|Iql+LRQxU8z|WB{OQ9dZ)Nk zaNHI8)^oSy;6PC`;mj;J;W$ixXv=tTq)M%BDoai>{e)m;Jex67v*WblTy4SOI4?Pq z3*9FxIlFd96@8H<2WPfcJ1hI*K&iSEpBa}-fS$b z%_rXYa^NMY9~;e)k)-evWi%rcqG@V0o8Dm1g}170Zr3Wv$W_x4-(Ptpd%Hf_NYa}6>QsckTtsbRWC^T1%nJ7$SkZ{Vd~JN9?{*bn=D z(%zkpy|6cA`@I9<^ag%Ey58D>mU+fWBlSjs(d-*17!2BmSxf8fBSh)0A2-w33lnhM z-iRY`o%EJ57PFg-H@VZ>SXxU3#C_ZkqpRU!5cq+gbUyc~&1i6H6vBQFjnEtWojPqK(KPP+ zop-q$m_qN=G#mteyE_Vkr00c$9$fx(r8kU%K@U4J@Zv#hgn1^d%W%ma^mmSb+IYt_ z`2Ubc9(}UwB|g)`UV_2LLqA;&;Wg5)qj5cRzPUH2XodvSUNp2s+ED>x~%jkUAO z+q_Tm4qIfCrVc_AgW3{N!kHW+clWDT-gZdc1!A4QojUa)}xvkT^&C`p^TZSK| zaWuJJdFjAGVi{xa8q9EztcG2G9L3i=_;O7t;;YH%Q6FM;qXT~o)*Vq08&NXC2ws2` zG%mX84>~+2m&#;d;(}KGG3Ocs?+JCPSfeZo#y7&~Y&;r5CBb4z;Q80tfnNAV_>gz~ z>@@9!(;zU0KCB;{kG%rzjR#h#*x&X0`#s-_gX^Vfw6Ox7ieLM_KgoWV2Km&^i~=lt z*}ps;4zbHwAlRrB9<)*azwle@8*k6dJwynZfj>q%OTWJ#B~G5a^TkXOgB#WH<>T7c-+2k`bB9k04Tzag5l1mQGzj_QMN$&n^@7OGQrsT$lKbVV*n2q=5 zA}4qE-n=&x1v&r7AC6|CNj`aWCl@WU@c2h_(I$&Ox_nP2I^^}g{`x(c=#qs?FXo~o z&z^i@7DZ0Uk29Ib$-!^#pT)!;CO12aS+W<%TX(;lX?C~Cs|$qO2jby} z9(m#YbqK%k$jx~Lk^udpTZizAj=Z^7hww{|d?HY^(9rmj0$BjP@`REE`ZAxaLnh{R z=4%Rsh9{n?LD`oXalWqj5&%8=v|@+_%{`+)nAgrX6kn({^G!v|1!(Vv0+|5)Q-K6( z&3;Qsf&=uT0tqzSdQO2*>z4{dF^S_hm5A8YH@>SFqSmhz$ix=D^gRW_PIO;ZAdEQs zisB31e_w^Lt0#V_Kz-a6qBz=H<2uC9*B7RWA(lONP0_*{pL$$@(ESTvsX^Iy2zmNh zMGLc^-c%r*(9v%z5U%f$?^nbIAWZd7KUeZXUz67rUpVbQ z|5h78&epZi{Y#G=BXa;%@BO6x{0zAlkoNEJK19gB0euYVy37;+03VA81ONa4009360763o0B!(V#K6EH!2khl zHl9L^3=9kb3=9nap@4ybfuSHLF*A>Wft7(lkVv3GWv~j#l^hUPvY@*%zn~;DIWcEI z-RTk%6&%PwjB9f;lZq0HDhJxNK28h(I$2uv0s;U4ABzYC000000RIL6LPG)o8vp|U L0000000000sC8el literal 0 HcmV?d00001 diff --git a/tests/dup-tags.bam.bai b/tests/dup-tags.bam.bai new file mode 100644 index 0000000000000000000000000000000000000000..0dd065fd03fa621ecfdd5d33785ec5cb89f5126b GIT binary patch literal 784 zcmZ>A^kfWY82W$-=mk$91_nm3SJF5jKDr8~y(Td*faDn%P{o-cir^%wUbyfmV+4f& E0IsS9X#fBK literal 0 HcmV?d00001 From 24c551bfbe667b5bd7c563fb1071d499184d23c1 Mon Sep 17 00:00:00 2001 From: Brad Langhorst Date: Thu, 27 Aug 2026 07:20:53 -0400 Subject: [PATCH 2/4] validate and split tags and values1;5u early, to save the : lookup in the loop --- functional-tests.sh | 5 +++++ mosdepth.nim | 22 +++++++++++++--------- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/functional-tests.sh b/functional-tests.sh index 831900c..572b0b1 100755 --- a/functional-tests.sh +++ b/functional-tests.sh @@ -105,6 +105,11 @@ run exclude_tag_multi $exe t tests/dup-tags.bam -F 772 --exclude-tag DT:SQ,DT:LB assert_exit_code 0 assert_equal "MT 0 80 1 MT 80 16569 0 " "$(zgrep ^MT t.per-base.bed.gz | tr -s '[:space:]' ' ')" +# a spec without a colon is rejected up front rather than silently ignored +run exclude_tag_malformed $exe t tests/dup-tags.bam --exclude-tag DT +assert_exit_code 2 +assert_in_stderr "--exclude-tag expects TAG:VALUE" + run bad_frag_len_filter $exe t tests/ovl.bam --min-frag-len 10 --max-frag-len 9 assert_in_stderr "--max-frag-len was lower than --min-frag-len." assert_exit_code 2 diff --git a/mosdepth.nim b/mosdepth.nim index be229cc..24437d7 100644 --- a/mosdepth.nim +++ b/mosdepth.nim @@ -251,7 +251,8 @@ proc coverage(bam: hts.Bam, arr: var coverage_t, region: var region_t, targets: seq[Target], mapq: int = -1, min_len: int = -1, max_len: int = int.high, eflag: uint16 = 1796, iflag: uint16 = 0, read_groups: seq[string] = (@[]), fast_mode: bool = false, - fragment_mode: bool = false, exclude_tags: seq[string] = (@[]), + fragment_mode: bool = false, + exclude_tags: seq[tuple[tag: string, value: string]] = (@[]), last_tid: var int = -1): int = # depth updates arr in-place and yields the tid for each chrom. # returns -1 if the chrom is not found in the bam header @@ -291,11 +292,9 @@ proc coverage(bam: hts.Bam, arr: var coverage_t, region: var region_t, # tag at all and the untagged majority must be kept. if has_exclude_tags: var drop = false - for spec in exclude_tags: - let c = spec.find(':') - if c <= 0: continue - var t = tag[string](rec, spec[0 ..< c]) - if t.isSome and t.get == spec[(c+1) .. ^1]: + for ex in exclude_tags: + var t = tag[string](rec, ex.tag) + if t.isSome and t.get == ex.value: drop = true break if drop: continue @@ -649,10 +648,15 @@ proc main(bam: hts.Bam, chrom: region_t, mapq: int, min_len: int, max_len: int, if $args["--read-groups"] != "nil": for r in ($args["--read-groups"]).split(','): read_groups.add($r) - var exclude_tags: seq[string] = @[] + var exclude_tags: seq[tuple[tag: string, value: string]] = @[] if $args["--exclude-tag"] != "nil": - for r in ($args["--exclude-tag"]).split(','): - exclude_tags.add($r) + for spec in ($args["--exclude-tag"]).split(','): + # split once here so the per-alignment check is a plain comparison. + let c = spec.find(':') + if c <= 0: + stderr.write_line("[mosdepth] error --exclude-tag expects TAG:VALUE, got '" & spec & "'") + quit(2) + exclude_tags.add((tag: spec[0 ..< c], value: spec[(c+1) .. ^1])) var levels = get_min_levels(targets) var chrom_region_distribution = newSeq[int64](region_distribution.len) From 7e1b71c3b7769e5754b86ef2d6bbc033c40a6d08 Mon Sep 17 00:00:00 2001 From: Brad Langhorst Date: Thu, 27 Aug 2026 07:33:36 -0400 Subject: [PATCH 3/4] CI: bump Nim 1.6.18 -> 1.6.20 in build matrix nim-unicodedb 0.14.1 now requires nim >= 1.6.20, and it is pulled in transitively via docopt -> regex, so the 1.6.18 matrix legs fail at dependency resolution with "Unsatisfied dependency: nim (>= 1.6.20)". 1.6.20 is the final 1.6.x release, so this keeps the same oldest-supported-Nim intent. --- .github/workflows/build.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 72f7d3d..29ff7ad 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -19,7 +19,7 @@ jobs: matrix: os: [ubuntu-22.04, macos-14] version: - - 1.6.18 + - 1.6.20 - 2.0.2 steps: From fea92a4f71343bb212612d4216bc83fd94579f86 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 9 Feb 2026 03:32:19 +0000 Subject: [PATCH 4/4] d4 upstream is impacted by a newly added rust lint... d4binding/src/api.rs:252 declares `memcpy` as `fn(*mut c_char, *const c_char, usize)` instead of `fn(*mut c_void, *const c_void, usize) -> *mut c_void`. Recent rustc denies that by default via the `invalid_runtime_symbol_definitions` lint, so the ubuntu runners (which ship a newer preinstalled toolchain than the macOS ones) fail with exit code 101 before mosdepth is built rather than patching upstream, i think we can ignore this lint... --- .github/workflows/build.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index 29ff7ad..a218846 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -116,6 +116,13 @@ jobs: fi - name: Install d4 + env: + # d4binding declares `memcpy` with a mismatched signature + # (d4binding/src/api.rs:252). Recent rustc denies that by default via + # the `invalid_runtime_symbol_definitions` lint, which breaks the + # upstream build. The mismatch is ABI-harmless, so downgrade it here + # until 38/d4-format fixes the declaration. + RUSTFLAGS: -A invalid_runtime_symbol_definitions run: | #export HTSLIB=system git clone https://github.com/38/d4-format