diff --git a/Dockerfile b/Dockerfile index d6244bd..7c61e66 100644 --- a/Dockerfile +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 3773a89..cbdfc27 100644 --- a/pyproject.toml +++ b/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', diff --git a/src/align_genotype/config_template.toml b/src/align_genotype/config_template.toml index e7da8a6..2e296ab 100644 --- a/src/align_genotype/config_template.toml +++ b/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 diff --git a/src/align_genotype/jobs/genotype.py b/src/align_genotype/jobs/genotype.py index 38b7cfb..c6fc795 100644 --- a/src/align_genotype/jobs/genotype.py +++ b/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 \\ diff --git a/src/align_genotype/scripts/duptest_to_cram.py b/src/align_genotype/scripts/duptest_to_cram.py new file mode 100644 index 0000000..ee4cb84 --- /dev/null +++ b/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) diff --git a/src/align_genotype/scripts/run_duptest_variant_calling.py b/src/align_genotype/scripts/run_duptest_variant_calling.py new file mode 100644 index 0000000..9327a97 --- /dev/null +++ b/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() diff --git a/src/align_genotype/scripts/run_multiple_dup_markers.py b/src/align_genotype/scripts/run_multiple_dup_markers.py new file mode 100644 index 0000000..41f1189 --- /dev/null +++ b/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 + + 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 number of unique reads to process. This sizes the bloom filter [default: 100000000] + -f, --false-pos-rate + Acceptable false positive rate for the bloom filter [default: 0.0001] + -r, --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)