olgabot/sourmash sig merge

.github/workflows/ci.yml

@@ -67,15 +67,18 @@ jobs:
           - "test --sketch_scaled 2,4 --sketch_num_hashes_log2 false"
           - "test_bam --barcodes_file false --rename_10x_barcodes false --save_fastas false --write_barcodes_meta_csv false"
           - "test_bam --rename_10x_barcodes false --write_barcodes_meta_csv false"
+          - "test_bam --skip_sig_merge"
           - "test_bam --write_barcodes_meta_csv false"
           - "test_bam --barcodes_file false --rename_10x_barcodes false"
           - "test_bam --rename_10x_barcodes false"
           - "test_fastas"
           - "test_protein_fastas"
           - "test_remove_ribo"
+          - "test_sig_merge"
           - "test_tenx_tgz"
           - "test_translate"
           - "test_translate_bam"
+          - "test_translate_bam --skip_sig_merge"
     steps:
       - name: Check out pipeline code
         uses: actions/checkout@v2

.vscode/settings.json

@@ -0,0 +1,3 @@
+{
+    "python.formatting.provider": "black"
+}

CHANGELOG.md

@@ -20,8 +20,10 @@ Initial release of nf-core/kmermaid, created with the [nf-core](https://nf-co.re
 barcode fastq
 * Add version printing for sencha, bam2fasta, and sourmash in Dockerfile, update versions in environment.yml
 * For processes translate, sourmash compute  add cpus=1 as they are only serial ([#107](https://github.com/nf-core/kmermaid/pull/107))
+* Add `sourmash sig merge` for aligned/unaligned signatures from bam files, and add `--skip_sig_merge` option to turn it off
 * Add `--protein_fastas` option for creating sketches of already-translated protein sequences
 * Add `--skip_compare option` to skip `sourmash_compare_sketches` process
+* Add merging of aligned/unaligned parts of single-cell data ([#117](https://github.com/nf-core/kmermaid/pull/117))
 * Add renamed package dependency orpheum (used to be known as sencha)
 
 ### `Fixed`
@@ -61,3 +63,5 @@ barcode fastq
 ### `Dependencies`
 
 ### `Deprecated`
+
+* Removed ability to specify multiple `--scaled` or `--num-hashes` values to enable merging of signatures

bin/merge_rename_sigs.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python
+import argparse
+import itertools
+import shutil
+
+import sourmash
+
+
+def merge(filenames, ksize, moltype, name=None, outsig=None):
+    """
+    merge one or more signatures.
+
+    Adapted from 'sourmash sig merge' command
+    """
+
+    first_sig = None
+    mh = None
+    total_loaded = 0
+
+    for sigfile in filenames:
+        this_n = 0
+
+        for sigobj in sourmash.load_file_as_signatures(
+            sigfile, ksize=ksize, select_moltype=moltype
+        ):
+
+            if sigobj is None:
+                error(
+                    "No signature in file {}",
+                    sigfile,
+                )
+
+            # first signature? initialize a bunch of stuff
+            if first_sig is None:
+                first_sig = sigobj
+                mh = first_sig.minhash.copy_and_clear()
+
+            try:
+                sigobj_mh = sigobj.minhash
+
+                mh.merge(sigobj_mh)
+            except:
+                error(
+                    "ERROR when merging signature '{}' ({}) from file {}",
+                    sigobj.name(),
+                    sigobj.md5sum()[:8],
+                    sigfile,
+                )
+                raise
+
+            this_n += 1
+            total_loaded += 1
+
+    merged_sigobj = sourmash.SourmashSignature(mh)
+    if name is not None:
+        merged_sigobj._name = name
+
+    if outsig is not None:
+        with open(outsig, "wt") as f:
+            sourmash.save_signatures([merged_sigobj], fp=f)
+
+    return merged_sigobj
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="""Merge signatures with same ksize and moltype"""
+    )
+    # Signature files
+    parser.add_argument("sigfiles", nargs="+")
+    parser.add_argument(
+        "--moltypes",
+        type=str,
+        help="Molecule types, comma-separated. e,g. 'protein,dayhoff'",
+        required=True,
+    )
+    parser.add_argument(
+        "--ksizes", type=str, help="K-mer sizes to combine", required=True
+    )
+
+    parser.add_argument(
+        "-o", "--outsig", type=str, help="Signature file to output to", required=True
+    )
+    parser.add_argument(
+        "-n",
+        "--name",
+        type=str,
+        help="Name of the signature to use",
+    )
+
+    args = parser.parse_args()
+
+    # Only iterate and read over the sigfiles if there is really something to merge
+    # "something to merge" = there is more than one sigfile. otherwise there's no point
+    # in reading in the files only to make the same file again
+    if len(args.sigfiles) > 1:
+        ksizes = map(int, args.ksizes.split(","))
+        moltypes = args.moltypes.split(",")
+
+        merged_sigobjs = []
+        for moltype, ksize in itertools.product(moltypes, ksizes):
+            merged = merge(
+                args.sigfiles, moltype=moltype, ksize=ksize, name=args.name
+            )
+            merged_sigobjs.append(merged)
+
+        with open(args.outsig, "wt") as f:
+            sourmash.save_signatures(merged_sigobjs, fp=f)
+    else:
+        # Otherwise, nothing to merge. Simply copy the file to the
+        # output signature location
+        shutil.copyfile(args.sigfiles[0], args.outsig)

bin/validate_sketch_value.py

@@ -0,0 +1,154 @@
+#!/usr/bin/env python
+from __future__ import print_function
+from collections import OrderedDict, defaultdict, Counter
+import logging
+import argparse
+import glob
+import os
+import sys
+
+
+# Create a logger
+logging.basicConfig(format="%(name)s - %(asctime)s %(levelname)s: %(message)s")
+logger = logging.getLogger(__file__)
+logger.setLevel(logging.INFO)
+
+
+def get_sketch_value(value, value_log2):
+    try:
+        if value:
+            if "," in value:
+                logger.exception(
+                    f"Can only provide a single number to --sketch_num_hashes or"
+                    f" --sketch_scaled. Provided '{value}"
+                )
+            sketch_value = int(value)
+        else:
+            if "," in value_log2:
+                logger.exception(
+                    f"Can only provide a single number to --sketch_num_hashes_log2 or"
+                    f" --sketch_scaled_log2. Provided '{value_log2}"
+                )
+            sketch_value = 2 ** int(value_log2)
+    except ValueError:
+        logger.exception(
+            "Can only supply a single value to --sketch_num_hashes, "
+            "--sketch_num_hashes_log2, --sketch_scaled, "
+            "--sketch_scaled_log2 "
+        )
+    return sketch_value
+
+
+def value_or_bool(value):
+    if value == "false":
+        return False
+    if value == "true":
+        logger.exception(
+            "Must set a value for the --sketch_num_hashes, "
+            "--sketch_num_hashes_log2, --sketch_scaled, "
+            "--sketch_scaled_log2 options! Cannot simply set the "
+            "flag. E.g. '--sketch_num_hashes 5' is valid but "
+            "'--sketch_num_hashes' on its own is not"
+        )
+        sys.exit(1)
+    else:
+        return value
+
+
+def main(
+    sketch_num_hashes,
+    sketch_num_hashes_log2,
+    sketch_scaled,
+    sketch_scaled_log2,
+    out,
+    sketch_style,
+):
+    sketch_num_hashes = value_or_bool(sketch_num_hashes)
+    sketch_num_hashes_log2 = value_or_bool(sketch_num_hashes_log2)
+    sketch_scaled = value_or_bool(sketch_scaled)
+    sketch_scaled_log2 = value_or_bool(sketch_scaled_log2)
+
+    using_size = sketch_num_hashes or sketch_num_hashes_log2
+    using_scaled = sketch_scaled or sketch_scaled_log2
+
+    if using_size and using_scaled:
+        logger.exception(
+            "Cannot specify both sketch scales and sizes! Can only"
+            " use one of --sketch_num_hashes, --sketch_num_hashes_log2, --sketch_scaled, "
+            "--sketch_scaled_log2. Exiting."
+        )
+        sys.exit(1)
+
+    if using_size:
+        if sketch_num_hashes and sketch_num_hashes_log2:
+            logger.exception(
+                "Cannot specify both --sketch_num_hashes and --sketch_num_hashes_log2! Exiting."
+            )
+            sys.exit(1)
+        sketch_value = get_sketch_value(sketch_num_hashes, sketch_num_hashes_log2)
+        with open(sketch_style, "w") as f:
+            f.write("size")
+    elif using_scaled:
+        if sketch_scaled and sketch_scaled_log2:
+            logger.exception(
+                "Cannot specify both --sketch_scaled and --sketch_scaled_log2! Exiting."
+            )
+            sys.exit(1)
+        sketch_value = get_sketch_value(sketch_scaled, sketch_scaled_log2)
+        with open(sketch_style, "w") as f:
+            f.write("scaled")
+
+    else:
+        logger.info(
+            "Did not specify a sketch size or scale with any of "
+            "--sketch_num_hashes, --sketch_num_hashes_log2, --sketch_scaled, --sketch_scaled_log2! "
+            "Falling back on sourmash's default of --sketch_scaled 500"
+        )
+        sketch_value = 500
+
+    with open(out, "w") as f:
+        f.write(f"{sketch_value}\n")
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="""Ensure that sketch sizes/scaleds provided are valid"""
+    )
+    parser.add_argument(
+        "--sketch_num_hashes", type=str, help="Flat size of the sketches"
+    )
+    parser.add_argument(
+        "--sketch_num_hashes_log2", type=str, help="Flat size of the sketches, log2"
+    )
+    parser.add_argument(
+        "--sketch_scaled", type=str, help="Fraction of total observed hashes to observe"
+    )
+    parser.add_argument(
+        "--sketch_scaled_log2",
+        type=str,
+        help="Fraction of total observed hashes to observe, log2",
+    )
+    parser.add_argument(
+        "-o",
+        "--output",
+        dest="output",
+        default="sketch_value.txt",
+        type=str,
+        help="file with output",
+    )
+    parser.add_argument(
+        "--sketch_style",
+        default="sketch_style.txt",
+        type=str,
+        help="file indicating 'size' or 'scaled'",
+    )
+
+    args = parser.parse_args()
+    main(
+        args.sketch_num_hashes,
+        args.sketch_num_hashes_log2,
+        args.sketch_scaled,
+        args.sketch_scaled_log2,
+        args.output,
+        args.sketch_style,
+    )