Test Script for various duplicate markers

Dockerfile

@@ -1,4 +1,4 @@
-FROM australia-southeast1-docker.pkg.dev/cpg-common/images/cpg_hail_gcloud:0.2.134.cpg2-2
+FROM australia-southeast1-docker.pkg.dev/cpg-common/images/cpg_hail_gcloud:0.2.137.cpg1-2
 
 ENV PYTHONDONTWRITEBYTECODE=1
 ENV VERSION=0.4.5

pyproject.toml

@@ -7,7 +7,7 @@ name='align_genotype'
 description='Dragmap align & genotype workflow, using CPG-Flow'
 readme = "README.md"
 # currently cpg-flow is pinned to this version
-requires-python = ">=3.10,<3.11"
+requires-python = ">=3.10,<3.12"
 version="0.4.5"
 license={ "file" = "LICENSE" }
 classifiers=[
@@ -21,7 +21,7 @@ classifiers=[
 ]
 
 dependencies=[
-    'cpg-flow~=1.2',
+    'cpg-flow~=1.3',
     'hatchling',
     'peddy',
     'slack_sdk',

src/align_genotype/config_template.toml

@@ -2,6 +2,8 @@
 # set this to false if we want to use non-preemptible machines
 align_spot = true
 
+sequencing_type = 'genome'
+
 # exome/genome workflows choose one of these default values
 exome_realignment_shards = 1
 genome_realignment_shards = 10
@@ -53,7 +55,7 @@ genome_cram_gb = 100
 downsample_reads = 20
 
 # Shard count for HaplotypeCaller
-scatter_count_genotype = 50
+scatter_count_genotype = 30
 
 # HaplotypeCaller resource options
 # Best performance for genomes - highmem, 4 cpus, 4 pairhmm threads

src/align_genotype/jobs/genotype.py

@@ -201,7 +201,6 @@ def _haplotype_caller_one(
     HaplotypeCaller \\
     -R {reference.base} \\
     -I {cram_resource_group.cram} \\
-    --read-index {cram_resource_group['cram.crai']} \\
     {f'-L {interval} ' if interval is not None else ''} \\
     --disable-spanning-event-genotyping \\
     --dragen-mode \\

src/align_genotype/scripts/duptest_to_cram.py

@@ -0,0 +1,82 @@
+"""
+Rinse the various output files from duplicate marking into definitely-3.0-CRAM files
+"""
+
+from argparse import ArgumentParser
+from hailtop.batch import Batch
+from hailtop.batch.job import BashJob
+
+from cpg_utils import hail_batch, to_path
+
+
+# I'm embedding this directly to make target images super obvious
+IMAGES: dict[str, str] = {
+    'rust_dupmark': 'australia-southeast1-docker.pkg.dev/cpg-common/images-dev/rust_dupmarker:0.1.0-1',
+    'samblaster': 'australia-southeast1-docker.pkg.dev/cpg-common/images/samblaster:0.1.26-1',
+    'sambamba': 'australia-southeast1-docker.pkg.dev/cpg-common/images/sambamba:1.0.1-1',
+    'streammd': 'australia-southeast1-docker.pkg.dev/cpg-common/images-dev/streammd:4.3.0-1',
+    'samtools': 'australia-southeast1-docker.pkg.dev/cpg-common/images/samtools:1.21-1',
+}
+
+OUTPUTS: dict[str, str] = {
+    'rust_dupmark': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/rust_dupmark/result.cram',
+    'sambamba': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/sambamba/result.bam',
+    # 'samblaster': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/samblaster/result.cram',  # already handled
+    # 'streammd': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/streammd/result.cram',  # already handled
+}
+
+needs_sorting = ['samblaster', 'streammd']
+
+
+def make_a_job(batch: Batch, tool: str) -> BashJob:
+    """ronseal."""
+    new_job = batch.new_bash_job(name=f'Convert: {tool}')
+    new_job.image(IMAGES['samtools'])
+    new_job.cpu(2)
+    new_job.memory('highmem')
+    new_job.storage('100GiB')
+    return new_job
+
+
+def main(outdir: str) -> None:
+    batch_instance = hail_batch.get_batch(name='Duplicate Marker outputs to CRAM.')
+
+    # set with workflow.ref_fasta
+    ref_fa = hail_batch.fasta_res_group(batch_instance).base
+
+    for tool, file in OUTPUTS.items():
+        outroot = f'{outdir}/{tool}/result'
+
+        if to_path(f'{outroot}.cram').exists():
+            continue
+
+        input_file = batch_instance.read_input(file)
+
+        job = make_a_job(batch_instance, tool)
+
+        job.declare_resource_group(
+            output={
+                'cram': '{root}.cram',
+                'cram.crai': '{root}.cram.crai',
+            },
+        )
+
+        job.command(f"""
+        samtools view --write-index -@2 \\
+            -T {ref_fa} \\
+            -O cram,version=3.0 \\
+            -o {job.output.cram} \\
+            {input_file}
+        echo "samtools view finished successfully"
+        """)
+
+        batch_instance.write_output(job.output, outroot)
+
+    batch_instance.run(wait=False)
+
+
+if __name__ == '__main__':
+    parser = ArgumentParser()
+    parser.add_argument('-o', help='Output root, all outputs derived from this path', required=True)
+    args = parser.parse_args()
+    main(outdir=args.o)

src/align_genotype/scripts/run_duptest_variant_calling.py

@@ -0,0 +1,42 @@
+"""
+NOT DONE YET
+
+Secondary part of the process - now that there are some duplicate-marked files (mix of CRAM and BAM),
+run them through the variant calling process we already have, to create a gVCF per dupmark tool
+"""
+
+from cpg_utils import hail_batch, to_path
+
+from align_genotype.jobs.genotype import genotype
+
+
+CRAMS = {
+    'rust_dupmark': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/cram30/rust_dupmark/result.cram',
+    'sambamba': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/cram30/sambamba/result.cram',
+    'samblaster': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/samblaster/result.cram',
+    # 'streammd': 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/streammd/result.cram.cram',  # not completed yet
+}
+
+
+def main() -> None:
+    batch_instance = hail_batch.get_batch(name='Various Duplicate Markers test.')
+
+    # get the intervals file from main workflow outputs (copied over)
+    intervals_path = 'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/intervals'
+    intervals_as_paths = list(to_path(intervals_path).glob('*.interval_list'))
+
+    for marker, cramfile in CRAMS.items():
+        _jobs = genotype(
+            sequencing_group_name=marker,
+            tmp_prefix=to_path(f'gs://cpg-ghfm-kidgen-test-tmp/duplicate_marker_test'),
+            cram_path=cramfile,
+            output_path=to_path(f'gs://cpg-ghfm-kidgen-test/duplicate_marker_test/haplotypecaller/{marker}.g.vcf.gz'),
+            intervals=intervals_as_paths,
+            job_attrs={'tool': marker}
+        )
+
+    batch_instance.run(wait=False)
+
+
+if __name__ == '__main__':
+    main()

src/align_genotype/scripts/run_multiple_dup_markers.py

@@ -0,0 +1,210 @@
+"""
+Intention here is to run a variety of tools for duplicate marking, and evaluate their runtime performance
+Eventually this will lead to a proper evaluation, comparing the resulting alignments
+
+Making some big storage/cpu/memory assumptions, and can scale all of that back down once we observe performance
+"""
+
+from argparse import ArgumentParser
+from hailtop.batch import Batch
+from hailtop.batch.job import BashJob
+
+from cpg_utils import hail_batch
+from cpg_flow import utils as flow_utils
+
+
+IMAGES: dict[str, str] = {
+    'rust_dupmark': 'australia-southeast1-docker.pkg.dev/cpg-common/images-dev/rust_dupmarker:0.1.0-1',
+    'samblaster': 'australia-southeast1-docker.pkg.dev/cpg-common/images/samblaster:0.1.26-1',
+    'sambamba': 'australia-southeast1-docker.pkg.dev/cpg-common/images/sambamba:1.0.1-1',
+    'streammd': 'australia-southeast1-docker.pkg.dev/cpg-common/images-dev/streammd:4.3.0-1',
+}
+
+
+def make_a_job(batch: Batch, tool: str) -> BashJob:
+    """ronseal."""
+    new_job = batch.new_bash_job(name=f'Tool: {tool}')
+    new_job.image(IMAGES[tool])
+    new_job.cpu(8)
+    new_job.memory('highmem')
+    new_job.storage('360GiB')
+    return new_job
+
+
+def sort_by_qname(batch: Batch, bamfile: str, outfile: str) -> BashJob | None:
+    """Some tools require the input to be grouped by QNAME, not coordinate."""
+
+    if flow_utils.exists(outfile):
+        return None
+
+    sb_job = make_a_job(batch, 'sambamba')
+    sb_job.command(f"""
+        sambamba sort \\
+            -n \\
+            -t 4 \\
+            --tmpdir=$BATCH_TMPDIR \\
+            {bamfile} \\
+        -o {sb_job.out}
+        """)
+    batch.write_output(sb_job.out, outfile)
+    return sb_job
+
+
+def create_sambamba_job(batch: Batch, bamfile: str, outfile: str) -> None:
+    """
+    https://lomereiter.github.io/sambamba/docs/sambamba-markdup.html
+    sambamba markdup OPTIONS <input.bam> <output.bam>
+
+    only works BAM -> BAM, needs further CRAM compression step
+    """
+
+    if flow_utils.exists(outfile):
+        return
+
+    sb_job = make_a_job(batch, 'sambamba')
+
+    sb_job.command(f"""
+    sambamba markdup \\
+        -t 4 \\
+        --tmpdir=$BATCH_TMPDIR \\
+        {bamfile} {sb_job.out}
+    """)
+    batch.write_output(sb_job.out, outfile)
+
+
+def create_samblaster_job(batch: Batch, bamfile: str, outfile: str, reference: str) -> BashJob | None:
+    """ronseal."""
+    if flow_utils.exists(f'{outfile}.cram'):
+        return None
+
+    sb_job = make_a_job(batch, 'samblaster')
+
+    sb_job.declare_resource_group(
+        output={
+            'cram': '{root}.cram',
+            'cram.crai': '{root}.cram.crai',
+        },
+    )
+
+    sb_job.command(f"""
+    mkdir $BATCH_TMPDIR/sort_tmp
+    samtools view \\
+        -h \\
+        --output-fmt SAM \\
+        {bamfile} | \\
+    samblaster | \\
+    samtools sort \\
+        --reference {reference} \\
+        -O CRAM,version=3.0 \\
+        -T $BATCH_TMPDIR/sort_tmp \\
+        --write-index \\
+        -o {sb_job.output.cram} \\
+        -
+    """)
+    batch.write_output(sb_job.output, outfile)
+
+    return sb_job
+
+
+def create_streammd_job(batch: Batch, bamfile: str, outfile: str, reference: str) -> BashJob | None:
+    """
+    e.g. bwa mem ref.fa r1.fq r2.fq|streammd
+    """
+
+    if flow_utils.exists(f'{outfile}.cram'):
+        return None
+
+    streammd_job = make_a_job(batch, 'streammd')
+
+    streammd_job.declare_resource_group(
+        output={
+            'cram': '{root}.cram',
+            'cram.crai': '{root}.cram.crai',
+        },
+    )
+
+    streammd_job.command(f"""
+    mkdir $BATCH_TMPDIR/sort_tmp
+    samtools view \\
+        -h \\
+        --output-fmt SAM \\
+        {bamfile} | \\
+    streammd | \\
+    samtools sort \\
+        --reference {reference} \\
+        -O CRAM,version=3.0 \\
+        -T $BATCH_TMPDIR/sort_tmp \\
+        --write-index \\
+        -o {streammd_job.output.cram} \\
+        -
+    """)
+
+    batch.write_output(streammd_job.output, outfile)
+    return streammd_job
+
+
+def create_dupmark_job(batch: Batch, bamfile: str, reference: str, outfile: str):
+    """
+    A high-performance SAM/BAM/CRAM duplicate marker using a Bloom filter.
+    Reads from standard input and writes CRAM to standard output
+
+    Usage: dupmark [OPTIONS]
+
+    Options:
+        -e, --expected-items <EXPECTED_ITEMS>
+        Expected number of unique reads to process. This sizes the bloom filter [default: 100000000]
+        -f, --false-pos-rate <FALSE_POS_RATE>
+        Acceptable false positive rate for the bloom filter [default: 0.0001]
+        -r, --reference <REFERENCE>
+        Optional reference FASTA for CRAM compression
+    """
+
+    if flow_utils.exists(outfile):
+        return
+
+    dm_job = make_a_job(batch, 'rust_dupmark')
+    dm_job.command(f"""
+    cat {bamfile} | dupmark --reference {reference} --expected-items 700000000 > {dm_job.out}
+    """)
+    batch.write_output(dm_job.out, outfile)
+
+
+def main(bamfile: str, outdir: str) -> None:
+    batch_instance = hail_batch.get_batch(name='Various Duplicate Markers test.')
+
+    # set with workflow.ref_fasta
+    ref_fa = hail_batch.fasta_res_group(batch_instance).base
+
+    input_bam = batch_instance.read_input(bamfile)
+
+    create_sambamba_job(batch_instance, input_bam, outfile=f'{outdir}/sambamba/result.bam')
+
+    create_dupmark_job(batch_instance, input_bam, outfile=f'{outdir}/rust_dupmark/result.cram', reference=ref_fa)
+
+    # these tools explicitly require QNAME sorted files as input
+    qname_sort_out = f'{outdir}/qname_sorted/result.bam'
+    qname_sort = sort_by_qname(batch=batch_instance, bamfile=input_bam, outfile=qname_sort_out)
+
+    qname_result = batch_instance.read_input(qname_sort_out)
+
+    streammd_job = create_streammd_job(
+        batch_instance, bamfile=qname_result, outfile=f'{outdir}/streammd/result', reference=ref_fa
+    )
+    if qname_sort and streammd_job:
+        streammd_job.depends_on(qname_sort)
+
+    samblaster_job = create_samblaster_job(
+        batch_instance, bamfile=qname_result, outfile=f'{outdir}/samblaster/result', reference=ref_fa
+    )
+    if qname_sort and samblaster_job:
+        samblaster_job.depends_on(qname_sort)
+
+    batch_instance.run(wait=False)
+
+
+if __name__ == '__main__':
+    parser = ArgumentParser()
+    parser.add_argument('-i', help='Sorted BAM file', required=True)
+    parser.add_argument('-o', help='output root, all outputs derived from this path', required=True)
+    args = parser.parse_args()
+    main(bamfile=args.i, outdir=args.o)