fix the ONT deduplication code

Author: SHuang-Broad

docker/lr-bam-dedup/Dockerfile

@@ -0,0 +1,18 @@
+FROM python:3.9.16-slim-bullseye
+
+COPY remove_duplicate_ont_aln.py /opt/
+
+RUN pip install pysam==0.21.0
+
+ARG DEBIAN_FRONTEND=noninteractive
+RUN apt-get -qqy update --fix-missing && \
+    apt-get -qqy dist-upgrade && \
+    apt-get -qqy install --no-install-recommends \
+                         samtools && \
+    apt-get -qqy clean && \
+    rm -rf /tmp/* \
+           /var/tmp/* \
+           /var/cache/apt/* \
+           /var/lib/apt/lists/* \
+           /usr/share/man/?? \
+           /usr/share/man/??_*

docker/lr-bam-dedup/Makefile

@@ -0,0 +1,12 @@
+VERSION = 0.1.1
+TAG1 = us.gcr.io/broad-dsp-lrma/lr-bam-dedup:$(VERSION)
+TAG2 = us.gcr.io/broad-dsp-lrma/lr-bam-dedup:latest
+
+all: build push
+
+build:
+	docker build -t $(TAG1) -t $(TAG2) .
+
+push:
+	docker push $(TAG1)
+	docker push $(TAG2)

docker/lr-bam-dedup/remove_duplicate_ont_aln.py

@@ -0,0 +1,76 @@
+import argparse
+import pysam
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Remove redundant alignment records from ONT BAM file',
+                                     prog='remove_redundant_reads')
+    parser.add_argument('-p', '--prefix', type=str, default="shard", help="Output prefix")
+    parser.add_argument('-a', '--annotations', type=str, help="Annotations on (potential) duplicate reads")
+
+    parser.add_argument('bam', type=str, help="BAM")
+    args = parser.parse_args()
+
+    # create a dict of set's, a trick to avoid Hash collisions
+    guilty_dict_per_chr = dict()
+    with open(args.annotations) as f:
+        for line in f:
+            arr = line.strip().split('\t')
+            name = arr[0]
+            chrom = arr[2]
+            guilty_dict_per_chr.setdefault(chrom, set())
+            guilty_dict_per_chr[chrom].add(name)
+
+    print("chromosomes on which there are duplicate records:")
+    print(f"  {guilty_dict_per_chr}")
+
+    # Silence message about the .bai file not being found.
+    pysam.set_verbosity(0)
+
+    num_alignments, num_dropped_alignments = 0, 0
+    duplicate_signatures = []
+    bf = pysam.Samfile(args.bam, 'rb', check_sq=False)
+    with pysam.Samfile(f'{args.prefix}.bam', 'wb', header=bf.header) as out:
+        # we rely on the observation that for coordinate sorted BAM,
+        # duplicate records will appear in blocks, hence once we step off a position with duplicates, we start afresh
+        current_position = -1
+        current_signatures = set()
+        for read in bf:
+            num_alignments += 1
+
+            chrom = read.reference_name
+            n = read.query_name
+            if chrom in guilty_dict_per_chr and n in guilty_dict_per_chr[chrom]:
+
+                mq = read.mapping_quality
+                sam_flag = read.flag
+                pos = read.reference_start
+                cigar = read.cigarstring
+                signature = f"{n}-{chrom}-{pos}-{mq}-{sam_flag}-{cigar}"
+
+                if current_position != pos:  # new position, let's write and reset
+                    out.write(read)
+                    current_position = pos
+                    current_signatures = set()
+                    current_signatures.add(signature)
+                elif signature in current_signatures:  # You're a duplicate record, and not appearing for the 1st time!
+                    num_dropped_alignments += 1
+                    duplicate_signatures.append(signature)  # same signature may appear more than twice, hence list and append
+                    pass
+                else:  # you are in a new group of duplicates that map to this location
+                    out.write(read)
+                    current_signatures.add(signature)
+            else:
+                out.write(read)
+
+    print(f'num_alignments: {num_alignments}')
+    print(f'num_dropped_alignments: {num_dropped_alignments}')
+    print(f'num_kept_alignments: {num_alignments - num_dropped_alignments}')
+
+    with open(f'{args.prefix}.duplicate.signatures.txt', 'w') as out:
+        for sig in duplicate_signatures:
+            out.write(f"{sig}\n")
+
+
+if __name__ == "__main__":
+    main()

docker/lr-utils/remove_duplicate_ont_aln.py

@@ -21,6 +21,9 @@ def main():
             guilty_dict_per_chr.setdefault(chrom, set())
             guilty_dict_per_chr[chrom].add(name)
 
+    print("chromosomes on which there are duplicate records:")
+    print(f"  {guilty_dict_per_chr}")
+
     # Silence message about the .bai file not being found.
     pysam.set_verbosity(0)
 
@@ -36,12 +39,13 @@ def main():
 
             chrom = read.reference_name
             n = read.query_name
-            if n in guilty_dict_per_chr[chrom]:
+            if chrom in guilty_dict_per_chr and n in guilty_dict_per_chr[chrom]:
 
                 mq = read.mapping_quality
                 sam_flag = read.flag
                 pos = read.reference_start
-                signature = f"{n}-{chrom}-{pos}-{mq}-{sam_flag}-"
+                cigar = read.cigarstring
+                signature = f"{n}-{chrom}-{pos}-{mq}-{sam_flag}-{cigar}"
 
                 if current_position != pos:  # new position, let's write and reset
                     out.write(read)

wdl/tasks/Utility/Utils.wdl

@@ -1950,11 +1950,11 @@ task FilterBamOnTag {
 task DeduplicateBam {
 
     meta {
-        description: "Utility to drop (occationally happening) duplicate records in input BAM"
+        description: "Utility to drop (occationally happening) literal duplicate records in input BAM"
     }
 
     parameter_meta {
-        aligned_bam: "input BAM file"
+        aligned_bam: "input BAM file (must be coordinate sorted)"
         aligned_bai: "input BAM index file"
         same_name_as_input: "if true, output BAM will have the same name as input BAM, otherwise it will have the input basename with .dedup suffix"
         runtime_attr_override: "override default runtime attributes"
@@ -1975,19 +1975,38 @@ task DeduplicateBam {
     String prefix = if (same_name_as_input) then base else (base + ".dedup")
 
     command <<<
+        set -eux
+
+        samtools view -H "~{aligned_bam}" | grep "@HD" > hd.line
+        if ! grep -qF "SO:coordinate" hd.line;
+        then
+            echo "BAM must be coordinate sorted!" && echo && cat hd.line && exit 1
+        fi
         echo "==========================================================="
         echo "collecting duplicate information"
         time \
             samtools view -@ 1 "~{aligned_bam}" | \
-            awk -F "\t" 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5}' | \
+            awk -F "\t" 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6}' | \
             sort | uniq -d \
             > "~{aligned_bam}".duplicates.txt
+
+        cnt=$(wc -l "~{aligned_bam}".duplicates.txt | awk '{print $1}')
+        if [[ ${cnt} -eq 0 ]];
+        then
+            echo "No duplicates found"
+            if ! ~{same_name_as_input} ;
+            then
+                mv "~{aligned_bam}" "~{prefix}.bam"
+                mv "~{aligned_bai}" "~{prefix}.bam.bai"
+            fi
+            exit 0
+        fi
         echo "==========================================================="
         echo "de-duplicating"
         time python3 /opt/remove_duplicate_ont_aln.py \
+            "~{aligned_bam}" \
             --prefix "~{prefix}" \
-            --annotations "~{aligned_bam}".duplicates.txt \
-            "~{aligned_bam}"
+            --annotations "~{aligned_bam}".duplicates.txt
         echo "==========================================================="
         echo "DONE"
         samtools index "~{prefix}.bam"
@@ -1996,6 +2015,7 @@ task DeduplicateBam {
     output {
         File corrected_bam = "~{prefix}.bam"
         File corrected_bai = "~{prefix}.bam.bai"
+        File duplicate_record_signatures = "~{prefix}.duplicate.signatures.txt"
     }
 
     #########################
@@ -2006,7 +2026,7 @@ task DeduplicateBam {
         boot_disk_gb:       10,
         preemptible_tries:  0,
         max_retries:        1,
-        docker:             "us.gcr.io/broad-dsp-lrma/lr-utils:0.1.10"
+        docker:             "us.gcr.io/broad-dsp-lrma/lr-bam-dedup:0.1.1"
     }
     RuntimeAttr runtime_attr = select_first([runtime_attr_override, default_attr])
     runtime {