Repository: whatshap/whatshap

Issue:
haplotag adds a second @PG line in the header with ID:whatshap when run on a file that was previously haplotagged.
We ran whatshap haplotag on bam A using a regular phased vcf file generating bam B
We then ran whatshap haplotag on bam B using a different phased vcf file generating bam C

IGV complained about having multiple PG lines with the same ID - when adding its PG line whatshap does not check for a previous existing PG line with the same ID.

- WhatsHap and Python version (`python -V` and `whatshap --version`)
Python 3.9.17
whatshap 2.1
- How you installed the tool (conda, pip)
pip install

Here are the PG lines from a bam file where whatshap haplotag was run on it twice.  Note that the first time the bam was split into 3 parts,  and we used samtools merge to combine them together.  Samtools merge dealt with the duplicate PG IDs by appending a random suffix.

> @PG     PN:whatshap     ID:whatshap     VN:1.4  CL:whatshap haplotag --tag-supplementary --output-threads 4 --reference /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-whatshap_haplotag/shard-0/inputs/p1gmcstorprimary01/smrtlink/pb-hosted-dataset/GRCh38/human_GRCh38_no_alt_analysis_set.fasta --output cmh006654-01A.m84052_230925_140850_s1.hifi_reads.bc2093.GRCh38.aligned.haplotagged.bam /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-phase_vcf/phase_vcf/12a93a5e-f0a6-4324-ab8d-46bff5419499/call-bcftools_concat/execution/cmh006654-01A.GRCh38.deepvariant.phased.vcf.gz /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-pbmm2_align/shard-0/execution/cmh006654-01A.m84052_230925_140850_s1.hifi_reads.bc2093.GRCh38.aligned.bam      m5:70b2b4e4583912cb7df227aaba920b8f
>  
> @PG     PN:whatshap     ID:whatshap-3BA06A1C    VN:1.4  CL:whatshap haplotag --tag-supplementary --output-threads 4 --reference /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-whatshap_haplotag/shard-1/inputs/p1gmcstorprimary01/smrtlink/pb-hosted-dataset/GRCh38/human_GRCh38_no_alt_analysis_set.fasta --output cmh006654-01A.m84052_230925_143914_s2.hifi_reads.bc2093.GRCh38.aligned.haplotagged.bam /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-phase_vcf/phase_vcf/12a93a5e-f0a6-4324-ab8d-46bff5419499/call-bcftools_concat/execution/cmh006654-01A.GRCh38.deepvariant.phased.vcf.gz /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-pbmm2_align/shard-1/execution/cmh006654-01A.m84052_230925_143914_s2.hifi_reads.bc2093.GRCh38.aligned.bam      m5:70b2b4e4583912cb7df227aaba920b8f
>  
> @PG     PN:whatshap     ID:whatshap-416AB9B1    VN:1.4  CL:whatshap haplotag --tag-supplementary --output-threads 4 --reference /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-whatshap_haplotag/shard-2/inputs/p1gmcstorprimary01/smrtlink/pb-hosted-dataset/GRCh38/human_GRCh38_no_alt_analysis_set.fasta --output cmh006654-01A.m84052_231005_211216_s3.hifi_reads.bc2093.GRCh38.aligned.haplotagged.bam /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-phase_vcf/phase_vcf/12a93a5e-f0a6-4324-ab8d-46bff5419499/call-bcftools_concat/execution/cmh006654-01A.GRCh38.deepvariant.phased.vcf.gz /cromwell-executions/humanwgs/29b3a689-17f3-4d15-b6b1-55f4c7fbb457/call-sample_analysis/shard-2/sample_analysis/1313bac2-e41f-45d6-8e46-a0799e675127/call-pbmm2_align/shard-2/execution/cmh006654-01A.m84052_231005_211216_s3.hifi_reads.bc2093.GRCh38.aligned.bam      m5:70b2b4e4583912cb7df227aaba920b8f
>  
> @PG     PN:whatshap     ID:whatshap     VN:2.1  CL:whatshap haplotag -o RESULT/aligned/cmh006654-01A.GRCh38.haplotagged.pedigree.bam --reference /reference/GRCh38/GRCh38_no_alt_analysis_set.decoy.fna /cluster/home/wacheung/analysis/scratch/cmh006654/cmh006654.cmh006654-01A.GRCh38.deepvariant.glnexus.whatshap-ped.vcf.gz /cluster/home/wacheung/analysis/scratch/cmh006654/cmh006654-01A.GRCh38.haplotagged.bam m5:37994ea46121627ad1938fd44bf76726

Discussion:

Patch:
diff --git a/whatshap/cli/haplotag.py b/whatshap/cli/haplotag.py
index bce2acc2..23922d4e 100644
--- a/whatshap/cli/haplotag.py
+++ b/whatshap/cli/haplotag.py
@@ -384,6 +384,13 @@ def open_output_alignment_file(aln_output, reference, vcf_md5, bam_header, threa
         "m5": vcf_md5,
     }
     if "PG" in bam_header:
+        # Make sure we don't have duplicate PG entries
+        nr = 1
+        whatshap_ids = {pg["ID"] for pg in bam_header["PG"] if pg["ID"].startswith("whatshap")}
+        while PG_entry["ID"] in whatshap_ids:
+            PG_entry["ID"] = f"whatshap.{nr}"
+            nr += 1
+
         bam_header["PG"].append(PG_entry)
     else:
         bam_header["PG"] = [PG_entry]


Test patch:
diff --git a/tests/test_run_haplotag.py b/tests/test_run_haplotag.py
index 97055a41..3c8fd281 100644
--- a/tests/test_run_haplotag.py
+++ b/tests/test_run_haplotag.py
@@ -578,3 +578,30 @@ def test_haplotag_duplicates_are_tagged(tmp_path):
             assert r1.get_tag("PS") == r2.get_tag("PS")
             count += 1
     assert count > 0
+
+
+def test_haplotag_run_twice(tmp_path):
+    outbam = tmp_path / "output.bam"
+    run_haplotag(
+        variant_file="tests/data/haplotag_sample.vcf.gz",
+        alignment_file="tests/data/haplotag_sample.bam",
+        given_samples=["mother"],
+        output=outbam,
+    )
+    # Index bam file
+    pysam.index(str(outbam))
+
+    outbam2 = tmp_path / "output2.bam"
+    run_haplotag(
+        variant_file="tests/data/haplotag_sample.vcf.gz",
+        alignment_file=outbam,
+        given_samples=["mother"],
+        output=outbam2,
+    )
+
+    # Check that there are two PG unique entries for whatshap
+    with pysam.AlignmentFile(outbam2) as f:
+        pg_entries = f.header.get("PG")
+        whatshap_ids = [entry["ID"] for entry in pg_entries if entry["ID"].startswith("whatshap")]
+        assert len(whatshap_ids) == 2
+        assert len(set(whatshap_ids)) == 2