Merge pull request #52 from ShawHahnLab/release-0.5.0

Version 0.5.0

Author: ressy

CHANGELOG.md

@@ -1,5 +1,31 @@
 # Changelog
 
+## 0.5.0 - 2023-01-04
+
+### Added
+
+ * `summarize` command will automatically use all available references for a
+   given species if a species is given but no references ([#50])
+ * `tree` command for creating and formatting phylogenetic trees ([#44])
+ * support for additional arguments for `getreads` command passed through to
+   bcl2fastq ([#43])
+ * `msa` command for building multiple sequence alignments with
+   [MUSCLE](https://drive5.com/muscle5/) ([#41])
+
+### Fixed
+
+ * `convert` command and underlying input/output features now handles sequence
+   descriptions ([#51])
+ * `identity` command now uses a custom sequence ID column if one is given
+   ([#49])
+
+[#51]: https://github.com/ShawHahnLab/igseq/pull/51
+[#50]: https://github.com/ShawHahnLab/igseq/pull/50
+[#49]: https://github.com/ShawHahnLab/igseq/pull/49
+[#44]: https://github.com/ShawHahnLab/igseq/pull/44
+[#43]: https://github.com/ShawHahnLab/igseq/pull/43
+[#41]: https://github.com/ShawHahnLab/igseq/pull/41
+
 ## 0.4.0 - 2022-09-17
 
 ### Added

conda/meta.yaml

@@ -1,5 +1,5 @@
 # https://docs.conda.io/projects/conda-build/en/latest/resources/define-metadata.html
-{% set version = "0.4.0" %}
+{% set version = "0.5.0" %}
 {% set build = "0" %}
 
 package:
@@ -43,13 +43,15 @@ requirements:
   # should follow the package match specifications."
   run:
     - python>=3.9.7,<4
-    - biopython>=1.79,<2
+    - biopython>=1.80,<2
     - python-newick>=1.3.2,<2
     - bwa>=0.7.17,<0.8
-    - samtools>=1.15,<2
+    - samtools>=1.16,<2
     - cutadapt>=3.7,<4
     - pear>=0.9.6,<1
     - igblast>=1.19,<2
+    - muscle>=5.1,<6
+    - fasttree>=2.1,<3
 
 test:
   source_files:

igseq/__main__.py

@@ -19,6 +19,8 @@
 from . import vdj_match
 from . import convert
 from . import identity
+from . import msa
+from . import tree
 from . import show
 from .util import IgSeqError
 from .version import __version__
@@ -77,10 +79,9 @@ def main(arglist=None):
         try:
             if args_extra:
                 # If there were unparsed arguments, see if we're in one of the
-                # commands (currently just igblast) that can take extra
-                # pass-through arguments.  If so pass them along, but if not,
-                # error out.
-                if args.func in [_main_igblast]:
+                # commands that can take extra pass-through arguments.  If so
+                # pass them along, but if not, error out.
+                if args.func in [_main_igblast, _main_getreads]:
                     args.func(args, args_extra)
                 else:
                     parser.parse_args(args_extra)
@@ -111,13 +112,14 @@ def main(arglist=None):
         except BrokenPipeError:
             os.dup2(devnull, sys.stderr.fileno())
 
-def _main_getreads(args):
+def _main_getreads(args, extra_args=None):
     if args.no_counts:
         args.countsfile = None
     getreads.getreads(
         path_input=args.input,
         dir_out=args.outdir,
         path_counts=args.countsfile,
+        extra_args=extra_args,
         threads_load=args.threads_load,
         threads_proc=args.threads,
         dry_run=args.dry_run)
@@ -240,6 +242,30 @@ def _main_identity(args):
         colmap=colmap,
         dry_run=args.dry_run)
 
+def _main_msa(args):
+    colmap = args_to_colmap(args)
+    msa.msa(
+        path_in=args.input,
+        path_out=args.output,
+        fmt_in=args.input_format,
+        fmt_out=args.output_format,
+        colmap=colmap,
+        dry_run=args.dry_run)
+
+def _main_tree(args):
+    colmap = args_to_colmap(args)
+    tree.tree(
+        paths_in=args.input,
+        path_out=args.output,
+        fmt_in=args.input_format,
+        fmt_out=args.output_format,
+        aligned=args.aligned,
+        pattern=args.set_pattern,
+        lists=args.set_list,
+        colors=args.set_color,
+        colmap=colmap,
+        dry_run=args.dry_run)
+
 def _setup_log(verbose, quiet, prefix):
     # Handle warnings via logging
     logging.captureWarnings(True)
@@ -306,6 +332,14 @@ def __setup_arg_parser():
         help="Calculate pairwise identities",
         description=rewrap(identity.__doc__),
         formatter_class=argparse.RawDescriptionHelpFormatter)
+    p_msa = subps.add_parser("msa",
+        help="Create multiple sequence alignments",
+        description=rewrap(msa.__doc__),
+        formatter_class=argparse.RawDescriptionHelpFormatter)
+    p_tree = subps.add_parser("tree",
+        help="Create and format phylogenetic trees",
+        description=rewrap(tree.__doc__),
+        formatter_class=argparse.RawDescriptionHelpFormatter)
     p_show = subps.add_parser("show",
         help="show file contents",
         description=rewrap(show.__doc__),
@@ -519,6 +553,54 @@ def __setup_arg_parser():
         help="Name of column containing sequences (for tabular input/output)")
     p_identity.set_defaults(func=_main_identity)
 
+    __add_common_args(p_msa)
+    p_msa.add_argument("input",
+        help="input file path, or a literal '-' for standard input")
+    p_msa.add_argument("output",
+        help="output file path, or a literal '-' for standard output")
+    p_msa.add_argument("--input-format",
+        help="format of input "
+        "(default: detect from input filename if possible)")
+    p_msa.add_argument("--output-format",
+        help="format of output "
+        "(default: detect from output filename if possible)")
+    p_msa.add_argument("--col-seq-id",
+        help="Name of column containing sequence IDs (for tabular input/output)")
+    p_msa.add_argument("--col-seq",
+        help="Name of column containing sequences (for tabular input/output)")
+    p_msa.add_argument("--col-seq-desc",
+        help="Name of column containing sequence descriptions (for tabular input/output)")
+    p_msa.set_defaults(func=_main_msa)
+
+    __add_common_args(p_tree)
+    p_tree.add_argument("input", nargs="+",
+        help="input file path, or a literal '-' for standard input")
+    p_tree.add_argument("output",
+        help="output file path, or a literal '-' for standard output")
+    p_tree.add_argument("--input-format",
+        help="format of input "
+        "(default: detect from input filename if possible)")
+    p_tree.add_argument("--output-format",
+        help="format of output "
+        "(default: detect from output filename if possible)")
+    p_tree.add_argument("--aligned", action=argparse.BooleanOptionalAction,
+        help="Explicitly specify if input is aligned or not, for sequence input "
+        "(default: guess from lengths)")
+    p_tree.add_argument("--col-seq-id",
+        help="Name of column containing sequence IDs (for tabular input)")
+    p_tree.add_argument("--col-seq",
+        help="Name of column containing sequences (for tabular input)")
+    p_tree.add_argument("--set-pattern", "-P",
+        help="regular expression to define set membership, with zero or one capture groups. "
+        "If a capture group is given, only that text is used to define the set names.")
+    p_tree.add_argument("--set-list", "-L", action="append",
+        help="filename containing a list of sequence IDs for a set.  "
+        "This can be given multiple times for multiple set/filename pairs.")
+    p_tree.add_argument("--set-color", "-C", action="append",
+        help="setname=colorcode, like set1=#ff0000, to override automatic set colors. "
+        "This can be given multiple times for multiple set/color pairs.")
+    p_tree.set_defaults(func=_main_tree)
+
     return parser
 
 def __add_common_args(obj):

igseq/colors.py

@@ -0,0 +1,45 @@
+
+# from SONAR's display_tree, which in turn got these from the iwanthue tool
+COLORS = ["#BE4229", "#E74721", "#8C431C", "#CB6A27", "#E98C25", "#946E13", "#D0A620", "#8B8A22", "#A9B71D", "#555C10", "#839A21", "#82B531", "#4D8426", "#55C42A", "#34992A", "#2B6B2E", "#4FC456", "#33B66D", "#4296CB", "#5A8CE2", "#3E5988", "#656CE2", "#524EA0", "#8F83CC", "#A57CE4", "#8E46AD", "#C056EB", "#CA6BE4", "#7B4D87", "#D186D7"]
+
+def merge_colors(colors, scale=0):
+    """Take an average of a list of colors and shift toward black.
+
+    More colors results in a darker result, up to the integer value given for
+    scale.  If scale is less than the number of colors this scaling is skipped.
+    """
+    result = [0, 0, 0]
+    if not colors:
+        return result
+    if len(colors) == 1:
+        return colors[0]
+    for color in colors:
+        for idx in range(3):
+            result[idx] += color[idx]
+    # not quite right, should rotate, really, not move directly toward the
+    # middle... but it'll do for now
+    if scale < len(colors):
+        scaling = 1
+    else:
+        scaling = ((scale - len(colors))/scale)**0.3
+    for idx in range(3):
+        result[idx] = result[idx] / len(colors)
+        result[idx] = int(result[idx] * scaling)
+    return result
+
+def color_str_to_trio(color_txt):
+    """Convert hex color string to trio of 0:255 ints."""
+    color_txt = color_txt.removeprefix("#")
+    # e.g. "ff0000"
+    if len(color_txt) == 6:
+        color = [int(color_txt[idx:(idx+2)], 16) for idx in range(0, 6, 2)]
+    # e.g. "f00" = "ff0000"
+    elif len(color_txt) == 3:
+        color = [int(color_txt[idx:(idx+1)]*2, 16) for idx in range(0, 3)]
+    else:
+        raise ValueError
+    return color
+
+def color_trio_to_str(color):
+    """Convert trio of 0:255 ints to hex color string."""
+    return "#" + "".join([f"{c:02x}" for c in color])

igseq/data/environment.yml

@@ -6,10 +6,12 @@ channels:
   - defaults
 dependencies:
   - python>=3.9.7,<4
-  - biopython>=1.79,<2
+  - biopython>=1.80,<2
   - python-newick>=1.3.2,<2
   - bwa>=0.7.17,<0.8
-  - samtools>=1.15,<2
+  - samtools>=1.16,<2
   - cutadapt>=3.7,<4
   - pear>=0.9.6,<1
   - igblast>=1.19,<2
+  - muscle>=5.1,<6
+  - fasttree>=2.1,<3

igseq/data/examples/inputs/msa/seqs.fasta

@@ -0,0 +1,4 @@
+>seq1
+ATCGCAAGACCCCTTGCTTGAGTATATCCC
+>seq2
+ATCGCAAGACCCCTTCTTGAGTATATCCC

igseq/data/examples/inputs/tree/seqs.fasta

@@ -0,0 +1,14 @@
+>wk16-001
+AGGGTCGTCACCATTTCCGCTAAGATCGTAACGACCATCGATGCAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGGATCTGCC
+>wk16-002
+AGGGTCGTCACCATTGCCGCTAAGATCGTAACGACCATCGATGCAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGGATCTGCC
+>wk16-003
+AGGGTCGTCACCATTTCCGCTAAGATCGTAACGACCATCGATGCAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGTCCTGGATCTGCC
+>wk20-001
+AGGGTCGTCACCATTTCCGCTAAGATCGTAACGACCATCGATGTAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGGATCTGCC
+>wk20-002
+AGGGTCGTCACCATTGCCGCTAAGATCGTAACGACCATCGATGCAGACGTATGGTATATGTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGGATCTGCC
+>wk20-003
+AGGGTCGTCACCATTTCCGCTACGATCGTAACGACCATCGATGCAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGTCCTGGATCTGCC
+>somethingelse
+AGGGCCGTCAGCATTACCGCTAAGGTCGTAATGACCATCTATGCACACGTATGGTATAGTTATTCGCCCACCAGCAGCTCGTGCAGAGCCTGGTTCTGCC

igseq/data/examples/inputs/tree/seqs2.fasta

@@ -0,0 +1,8 @@
+>wk24-001
+AGGGTCGTCACTATTTCCGCTAAGATCGTAACGACCATCGATGTAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGGATCTGCC
+>wk24-002
+AGGGTCGTCACCATTTCCGCTAAGATCGTAACGACCATCGATGTAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGGCCTGCATCTGCC
+>wk24-003
+AGGCTCGTCACCATTGCCGCTAAGATCGTAACGACCATCGATGCAGACGTATGGTATATGTATTCGCCGACCAGAAGCTCATGCAGGGCCTGGATCTGCC
+>wk24-004
+AGGCTCGTCACCATTTCCGCTACGATCGTAACGACCATCGATGCAGACGTATGGTATATTTATTCGCCGACCAGAAGCTCTTGCAGGTCCTGGATCTGCC

igseq/data/examples/inputs/tree/wk16.txt

@@ -0,0 +1,3 @@
+wk16-001
+wk16-002
+wk16-003

igseq/data/examples/inputs/tree/wk20.txt

@@ -0,0 +1,3 @@
+wk20-001
+wk20-002
+wk20-003