From e2585dcd1180291ccebcf7ee632ab57f45c84389 Mon Sep 17 00:00:00 2001
From: Neil Horner <neil.horner@nanoporetech.com>
Date: Thu, 30 May 2024 19:40:11 +0000
Subject: [PATCH] Template update 2024 05 24

---
 bin/workflow_glue/configure_igv.py       | 205 +++++++++++++++++++++++
 bin/workflow_glue/get_max_depth_locus.py |  61 +++++++
 lib/common.nf                            |  56 +++++++
 nextflow.config                          |   2 +-
 4 files changed, 323 insertions(+), 1 deletion(-)
 create mode 100755 bin/workflow_glue/configure_igv.py
 create mode 100755 bin/workflow_glue/get_max_depth_locus.py
 create mode 100644 lib/common.nf

diff --git a/bin/workflow_glue/configure_igv.py b/bin/workflow_glue/configure_igv.py
new file mode 100755
index 0000000..45b4ebd
--- /dev/null
+++ b/bin/workflow_glue/configure_igv.py
@@ -0,0 +1,205 @@
+"""Create an IGV config file."""
+
+from itertools import zip_longest
+import json
+import sys
+
+from .util import get_named_logger, wf_parser  # noqa: ABS101
+
+
+def parse_fnames(fofn):
+    """Parse list with filenames and return them grouped as ref-, BAM-, or VCF-related.
+
+    :param fofn: File with list of file names (one per line)
+    :return: dict of reference-related filenames (with keys 'ref', 'fai', and '.gzi' and
+        `None` as default values); lists of BAM- and VCF-related filenames
+    """
+    ref_extensions = [".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fna.gz"]
+    ref_dict = {}
+    bams = []
+    bam_indices = []
+    vcfs = []
+    vcf_indices = []
+    with open(fofn, "r") as f:
+        for line in f:
+            fname = line.strip()
+            if any(fname.endswith(ext) for ext in ref_extensions):
+                ref_dict["ref"] = fname
+            elif fname.endswith(".fai"):
+                ref_dict["fai"] = fname
+            elif fname.endswith(".gzi"):
+                ref_dict["gzi"] = fname
+            elif fname.endswith(".bam"):
+                bams.append(fname)
+            elif fname.endswith(".bai"):
+                bam_indices.append(fname)
+            elif fname.endswith(".vcf") or fname.endswith(".vcf.gz"):
+                vcfs.append(fname)
+            elif fname.endswith(".csi") or fname.endswith(".tbi"):
+                vcf_indices.append(fname)
+    # do some sanity checks
+    if "ref" not in ref_dict:
+        raise ValueError(
+            "No reference file (i.e. file ending in one of "
+            f"{ref_extensions} was found)."
+        )
+    ref = ref_dict["ref"]
+    if (gzi := ref_dict.get("gzi")) is not None:
+        # since we got a '.gzi' index, make sure that the reference is actually
+        # compressed
+        if not ref_dict["ref"].endswith(".gz"):
+            raise ValueError(
+                f"Found GZI reference index '{gzi}', but the reference file "
+                f"'{ref}' appears not to be compressed."
+            )
+    if bam_indices:
+        if len(bams) != len(bam_indices):
+            raise ValueError("Got different number of BAM and BAM index files.")
+    if vcf_indices:
+        if len(vcfs) != len(vcf_indices):
+            raise ValueError("Got different number of VCF and VCF index files.")
+    if bams and vcfs:
+        if len(bams) != len(vcfs):
+            raise ValueError("Got different number of BAM and VCF files.")
+    # if we got BAM or VCF indices, pair them up with their corresponding files (and
+    # otherwise with `None`)
+    bams_with_indices = zip_longest(bams, bam_indices)
+    vcfs_with_indices = zip_longest(vcfs, vcf_indices)
+    return ref_dict, bams_with_indices, vcfs_with_indices
+
+
+def get_reference_options(ref, fai=None, gzi=None):
+    """Create dict with IGV reference options.
+
+    :param ref: reference file name
+    :param fai: name reference `.fai` index file
+    :param gzi: name of `.gzi` index file for a compressed reference
+    :return: dict with reference options
+    """
+    # initialise the options dict and add the index attributes later
+    ref_opts = {
+        "id": "ref",
+        "name": "ref",
+        "wholeGenomeView": False,
+        "fastaURL": ref,
+    }
+    if fai is not None:
+        ref_opts["indexURL"] = fai
+    if gzi is not None:
+        ref_opts["compressedIndexURL"] = gzi
+    return ref_opts
+
+
+def get_alignment_track(bam, bai=None, extra_opts=None):
+    """Create dict with options for IGV alignment track.
+
+    :param bam: name of BAM file to be displayed
+    :param bai: name of BAM index file
+    :param extra_opts: dict of extra options for the alignment track
+    :return: dict with alignment track options
+    """
+    alignment_track_dict = {
+        "name": bam,
+        "type": "alignment",
+        "format": "bam",
+        "url": bam,
+    }
+    # add the BAM index if present
+    if bai is not None:
+        alignment_track_dict["indexURL"] = bai
+    alignment_track_dict.update(extra_opts or {})
+    return alignment_track_dict
+
+
+def get_variant_track(vcf, index=None, extra_opts=None):
+    """Create dict with options for IGV variant track.
+
+    :param vcf: name of VCF file to be displayed
+    :param index: name of VCF index file (ending in `.csi` or `.tbi`)
+    :param extra_opts: dict of extra options for the variant track
+    :return: dict with variant track options
+    """
+    variant_track_dict = {
+        "name": vcf,
+        "type": "variant",
+        "format": "vcf",
+        "url": vcf,
+    }
+    # add the VCF index if we got an index extension
+    if index is not None:
+        variant_track_dict["indexURL"] = index
+    variant_track_dict.update(extra_opts or {})
+    return variant_track_dict
+
+
+def main(args):
+    """Run the entry point."""
+    logger = get_named_logger("configIGV")
+
+    # parse the FOFN
+    ref_dict, bams_with_indices, vcfs_with_indices = parse_fnames(args.fofn)
+
+    # initialise the IGV options dict with the reference options
+    json_dict = {"reference": get_reference_options(**ref_dict)}
+
+    # if we got JSON files with extra options for the alignment / variant tracks, read
+    # them
+    extra_alignment_opts = {}
+    if args.extra_bam_opts is not None:
+        with open(args.extra_bam_opts, "r") as f:
+            extra_alignment_opts = json.load(f)
+    extra_variant_opts = {}
+    if args.extra_vcf_opts is not None:
+        with open(args.extra_vcf_opts, "r") as f:
+            extra_variant_opts = json.load(f)
+
+    # now add the alignment and variant tracks
+    json_dict["tracks"] = []
+    # we use `zip_longest` to make sure that variant and alignment tracks from the same
+    # sample are added after each other
+    for (vcf, vcf_index), (bam, bam_index) in zip_longest(
+        vcfs_with_indices, bams_with_indices, fillvalue=(None, None)
+    ):
+        if vcf is not None:
+            # add an variant track for the VCF
+            json_dict["tracks"].append(
+                get_variant_track(vcf, vcf_index, extra_variant_opts)
+            )
+        if bam is not None:
+            # add an alignment track for the BAM
+            json_dict["tracks"].append(
+                get_alignment_track(bam, bam_index, extra_alignment_opts)
+            )
+
+    if args.locus is not None:
+        json_dict["locus"] = args.locus
+
+    json.dump(json_dict, sys.stdout, indent=4)
+
+    logger.info("Printed IGV config JSON to STDOUT.")
+
+
+def argparser():
+    """Argument parser for entrypoint."""
+    parser = wf_parser("configure_igv")
+    parser.add_argument(
+        "--fofn",
+        required=True,
+        help=(
+            "File with list of names of reference / BAM / VCF files and indices "
+            "(one filename per line)"
+        ),
+    )
+    parser.add_argument(
+        "--locus",
+        help="Locus string to set initial genomic coordinates to display in IGV",
+    )
+    parser.add_argument(
+        "--extra-bam-opts",
+        help="JSON file with extra options for alignment tracks",
+    )
+    parser.add_argument(
+        "--extra-vcf-opts",
+        help="JSON file with extra options for variant tracks",
+    )
+    return parser
diff --git a/bin/workflow_glue/get_max_depth_locus.py b/bin/workflow_glue/get_max_depth_locus.py
new file mode 100755
index 0000000..eaa216d
--- /dev/null
+++ b/bin/workflow_glue/get_max_depth_locus.py
@@ -0,0 +1,61 @@
+"""Find max depth window in a `mosdepth` regions BED file and write as locus string."""
+
+from pathlib import Path
+import sys
+
+import pandas as pd
+
+from .util import get_named_logger, wf_parser  # noqa: ABS101
+
+
+def main(args):
+    """Run the entry point."""
+    logger = get_named_logger("getMaxDepth")
+
+    # read the regions BED file
+    df = pd.read_csv(
+        args.depths_bed, sep="\t", header=None, names=["ref", "start", "end", "depth"]
+    )
+
+    # get the window with the largest depth
+    ref, start, end, depth = df.loc[df["depth"].idxmax()]
+
+    # get the length of the reference of that window
+    ref_length = df.query("ref == @ref")["end"].iloc[-1]
+
+    # show the whole reference in case it's shorter than the desired locus size
+    if ref_length < args.locus_size:
+        start = 1
+        end = ref_length
+    else:
+        # otherwise, show a region of the desired size around the window
+        half_size = args.locus_size // 2
+        mid = (start + end) // 2
+        start = mid - half_size
+        end = mid + half_size
+        # check if the region starts below `1` or ends beyond the end of the reference
+        if start < 1:
+            start = 1
+            end = args.locus_size
+        if end > ref_length:
+            start = ref_length - args.locus_size
+            end = ref_length
+
+    # write depth and locus string
+    sys.stdout.write(f"{depth}\t{ref}:{start}-{end}")
+
+    logger.info("Wrote locus with maximum depth to STDOUT.")
+
+
+def argparser():
+    """Argument parser for entrypoint."""
+    parser = wf_parser("check_bam_headers")
+    parser.add_argument(
+        "depths_bed",
+        type=Path,
+        help="path to mosdepth regions depth file (can be compressed)",
+    )
+    parser.add_argument(
+        "locus_size", type=int, help="size of the locus in basepairs (e.g. '2000')"
+    )
+    return parser
diff --git a/lib/common.nf b/lib/common.nf
new file mode 100644
index 0000000..9dfd5a7
--- /dev/null
+++ b/lib/common.nf
@@ -0,0 +1,56 @@
+import groovy.json.JsonBuilder
+
+process getParams {
+    label "wf_common"
+    cpus 1
+    memory "2 GB"
+    output:
+        path "params.json"
+    script:
+        def paramsJSON = new JsonBuilder(params).toPrettyString()
+    """
+    # Output nextflow params object to JSON
+    echo '$paramsJSON' > params.json
+    """
+}
+
+process configure_igv {
+    label "wf_common"
+    cpus 1
+    memory "2 GB"
+    input:
+        // the python script will work out what to do with all the files based on their
+        // extensions
+        path "file-names.txt"
+        val locus_str
+        val bam_extra_opts
+        val vcf_extra_opts
+    output: path "igv.json"
+    script:
+    // the locus argument just makes sure that the initial view in IGV shows something
+    // interesting
+    String locus_arg = locus_str ? "--locus $locus_str" : ""
+    // extra options for alignment tracks
+    def bam_opts_json_str = \
+        bam_extra_opts ? new JsonBuilder(bam_extra_opts).toPrettyString() : ""
+    String bam_extra_opts_arg = \
+        bam_extra_opts ? "--extra-bam-opts bam-extra-opts.json" : ""
+    // extra options for variant tracks
+    def vcf_opts_json_str = \
+        vcf_extra_opts ? new JsonBuilder(vcf_extra_opts).toPrettyString() : ""
+    String vcf_extra_opts_arg = \
+        vcf_extra_opts ? "--extra-vcf-opts vcf-extra-opts.json" : ""
+    """
+    # write out JSON files with extra options for the alignment and variant tracks
+    echo '$bam_opts_json_str' > bam-extra-opts.json
+    echo '$vcf_opts_json_str' > vcf-extra-opts.json
+
+    workflow-glue configure_igv \
+        --fofn file-names.txt \
+        $locus_arg \
+        $bam_extra_opts_arg \
+        $vcf_extra_opts_arg \
+    > igv.json
+    """
+}
+
diff --git a/nextflow.config b/nextflow.config
index e2805fe..33e5fd2 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -105,7 +105,7 @@ params {
       ]
     agent = null
     container_sha = "shae7c9f184996a384e99be68e790f0612f0c732867"
-	common_sha = "sha91cd87900c86f05bf36d8c77b841b8fda5ecf3aa"
+	common_sha = "sha338caea0a2532dc0ea8f46638ccc322bb8f9af48"
     }
 }