Merge pull request #40 from ShawHahnLab/release-0.4.0

Release 0.4.0

Author: ressy

CHANGELOG.md

@@ -1,5 +1,29 @@
 # Changelog
 
+## 0.4.0 - 2022-09-17
+
+### Added
+
+ * Automatic usage of all available references for a given species in `igblast`
+   command ([#39])
+ * `identity` command for calculating pairwise identity between arbitrary
+   queries and references ([#31], [#37])
+ * Support for showing basic tree topology for Newick-format files in `show`
+   command ([#33])
+
+### Fixed
+
+ * Uppercase file extensions are now supported by the `show` command ([#35])
+ * broken pipes (such as from `igseq something | something else`) are now
+   handled gracefully ([#30])
+
+[#39]: https://github.com/ShawHahnLab/igseq/pull/39
+[#37]: https://github.com/ShawHahnLab/igseq/pull/37
+[#35]: https://github.com/ShawHahnLab/igseq/pull/35
+[#33]: https://github.com/ShawHahnLab/igseq/pull/33
+[#31]: https://github.com/ShawHahnLab/igseq/pull/31
+[#30]: https://github.com/ShawHahnLab/igseq/pull/30
+
 ## 0.3.0 - 2022-07-14
 
 ### Added

README.md

@@ -1,5 +1,7 @@
 # IgSeq Utilities
 
+[![CircleCI Build Status](https://circleci.com/gh/ShawHahnLab/igseq/tree/dev.svg?style=svg)](https://circleci.com/gh/ShawHahnLab/igseq/tree/dev)
+
 **Work in Progress**
 
 A command-line tool for various common Ig-Seq tasks.  These are heavily biased

conda/meta.yaml

@@ -1,5 +1,5 @@
 # https://docs.conda.io/projects/conda-build/en/latest/resources/define-metadata.html
-{% set version = "0.3.0" %}
+{% set version = "0.4.0" %}
 {% set build = "0" %}
 
 package:
@@ -44,11 +44,12 @@ requirements:
   run:
     - python>=3.9.7,<4
     - biopython>=1.79,<2
+    - python-newick>=1.3.2,<2
     - bwa>=0.7.17,<0.8
-    - samtools>=1.13,<2
-    - cutadapt>=3.5,<4
+    - samtools>=1.15,<2
+    - cutadapt>=3.7,<4
     - pear>=0.9.6,<1
-    - igblast>=1.17.1,<2
+    - igblast>=1.19,<2
 
 test:
   source_files:

igseq/__main__.py

@@ -2,6 +2,7 @@
 Utilities for common IgSeq tasks.
 """
 
+import os
 import sys
 import argparse
 import logging
@@ -17,6 +18,7 @@
 from . import vdj_gather
 from . import vdj_match
 from . import convert
+from . import identity
 from . import show
 from .util import IgSeqError
 from .version import __version__
@@ -72,22 +74,42 @@ def main(arglist=None):
         prefix = "[DRYRUN] "
     _setup_log(args.verbose, args.quiet, prefix)
     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]:
-                args.func(args, args_extra)
+        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]:
+                    args.func(args, args_extra)
+                else:
+                    parser.parse_args(args_extra)
             else:
-                parser.parse_args(args_extra)
-        else:
-            args.func(args)
-    except IgSeqError as err:
-        sys.stderr.write(
-            f"\nigseq failed because: {err.message}\n"
-            "Considering adding -v or -vv to the command if the problem isn't clear.\n")
-        sys.exit(1)
+                args.func(args)
+        except IgSeqError as err:
+            sys.stderr.write(
+                f"\nigseq failed because: {err.message}\n"
+                "Considering adding -v or -vv to the command if the problem isn't clear.\n")
+            sys.exit(1)
+        sys.stdout.flush()
+        sys.stderr.flush()
+    except BrokenPipeError:
+        # If stdout and/or stderr were writing to a pipe and that pipe is now
+        # closed, we'll swap in /dev/null for whichever it is to handle this
+        # quietly and to prevent it from arising again when Python tries to
+        # flush file handles on exit.
+        # Adapted from
+        # https://stackoverflow.com/questions/26692284
+        # https://docs.python.org/3/library/signal.html#note-on-sigpipe
+        devnull = os.open(os.devnull, os.O_WRONLY)
+        try:
+            sys.stdout.flush()
+        except BrokenPipeError:
+            os.dup2(devnull, sys.stdout.fileno())
+        try:
+            sys.stderr.flush()
+        except BrokenPipeError:
+            os.dup2(devnull, sys.stderr.fileno())
 
 def _main_getreads(args):
     if args.no_counts:
@@ -207,6 +229,17 @@ def _main_convert(args):
         dummyqual=args.dummy_qual,
         dry_run=args.dry_run)
 
+def _main_identity(args):
+    colmap = args_to_colmap(args)
+    identity.identity(
+        path_in=args.input,
+        path_out=args.output,
+        path_ref=args.reference,
+        fmt_in=args.input_format,
+        fmt_in_ref=args.ref_format,
+        colmap=colmap,
+        dry_run=args.dry_run)
+
 def _setup_log(verbose, quiet, prefix):
     # Handle warnings via logging
     logging.captureWarnings(True)
@@ -269,6 +302,10 @@ def __setup_arg_parser():
         help="Convert FASTA/FASTQ/CSV/TSV",
         description=rewrap(convert.__doc__),
         formatter_class=argparse.RawDescriptionHelpFormatter)
+    p_identity = subps.add_parser("identity",
+        help="Calculate pairwise identities",
+        description=rewrap(identity.__doc__),
+        formatter_class=argparse.RawDescriptionHelpFormatter)
     p_show = subps.add_parser("show",
         help="show file contents",
         description=rewrap(show.__doc__),
@@ -463,6 +500,25 @@ def __setup_arg_parser():
         'as text (e.g. use "I" for 40)')
     p_convert.set_defaults(func=_main_convert)
 
+    __add_common_args(p_identity)
+    p_identity.add_argument("input",
+        help="input file path, or a literal '-' for standard input")
+    p_identity.add_argument("output",
+        help="output file path, or a literal '-' for standard output")
+    p_identity.add_argument("-r", "--reference",
+            help="optional reference file path (default: use first query as ref)")
+    p_identity.add_argument("--input-format",
+        help="format of input "
+        "(default: detect from input filename if possible)")
+    p_identity.add_argument("--ref-format",
+        help="format of reference "
+        "(default: detect from reference filename if possible)")
+    p_identity.add_argument("--col-seq-id",
+        help="Name of column containing sequence IDs (for tabular input/output)")
+    p_identity.add_argument("--col-seq",
+        help="Name of column containing sequences (for tabular input/output)")
+    p_identity.set_defaults(func=_main_identity)
+
     return parser
 
 def __add_common_args(obj):

igseq/data/environment.yml

@@ -7,8 +7,9 @@ channels:
 dependencies:
   - python>=3.9.7,<4
   - biopython>=1.79,<2
+  - python-newick>=1.3.2,<2
   - bwa>=0.7.17,<0.8
-  - samtools>=1.13,<2
-  - cutadapt>=3.5,<4
+  - samtools>=1.15,<2
+  - cutadapt>=3.7,<4
   - pear>=0.9.6,<1
-  - igblast>=1.17.1,<2
+  - igblast>=1.19,<2

igseq/data/examples/identity.sh

@@ -0,0 +1,46 @@
+#!/usr/bin/env bash
+
+[ -v EXAMPLES ] || EXAMPLES=$(python -c 'import igseq.util; print(igseq.util.DATA)')/examples
+
+# A few example query sequences derived from the igblast example
+QUERY=$EXAMPLES/inputs/identity/query.fasta
+# AIRR TSV table for those sequences
+QUERY_AIRR=$EXAMPLES/inputs/identity/query_airr.tsv
+# A contrived reference (slightly altered version of the first query)
+REF=$EXAMPLES/inputs/identity/ref.fasta
+# Just the junction region from that sequence, as NT and AA
+REF_JUNCTION=$EXAMPLES/inputs/identity/ref_junction.fasta
+REF_JUNCTION_AA=$EXAMPLES/inputs/identity/ref_junction_aa.fasta
+# AIRR TSV for the reference
+REF_AIRR=$EXAMPLES/inputs/identity/ref_airr.tsv
+
+# Check identity of each query to the one reference, with each row in the CSV
+# output giving query ID, reference ID, and fraction of matching bases
+igseq identity -r $REF $QUERY idents.csv
+# use "igseq show" command for a human-readable view of the CSV file
+igseq show idents.csv
+
+# With no reference specified the first query sequence will be used as the reference
+igseq identity $QUERY idents_first_query.csv
+
+# To compare all-to-all, explicitly give the query path as the reference too
+igseq identity -r $QUERY $QUERY idents_all.csv
+
+# We can use a "-" to read from/to standard input and/or output, but we need to
+# specify the input format for that case
+igseq identity -r $REF --input-format fa - - < $QUERY > idents_stdout.csv
+
+# We can work with the formats supported by the convert command, like AIRR TSV
+# for example.  As with that command, the defaults are to use the sequence_id
+# and sequence columns.
+igseq identity -r $REF $QUERY_AIRR idents_from_airr.csv
+
+# We can specify other columns to use instead
+igseq identity -r $REF_JUNCTION --col-seq junction $QUERY_AIRR idents_from_airr_junctions.csv
+
+# AA works too
+igseq identity -r $REF_JUNCTION_AA --col-seq junction_aa $QUERY_AIRR idents_from_airr_junctions_aa.csv
+
+# We can also give the reference in a format other than FASTA, though the same
+# column settings are used for query and reference
+igseq identity -r $REF_AIRR --col-seq junction_aa $QUERY_AIRR idents_from_airr_junctions_aa_2.csv

igseq/data/examples/inputs/identity/query.fasta

@@ -0,0 +1,24 @@
+>query1
+CAGCTGCAGCTGCAGGAGTCGGGCCCAGGACTGGTGAAGCCTTCGGAGACCCTGTCCCTC
+ACCTGCGCTGTCTCTGGTGGCTCCATCAGCAGTAACTACTGGAGCTGGATCCGCCAGCCC
+CCAGGGAAGGGACTGGAGTGGATTGGACGTATCTCTGGTAGTGGTGGGAGCCCCGACTAC
+AACCCCTCCCTCAAGAGTCGAGTCACCATTTCAACAGACACGTCCAAGAGCCAGTTCTCC
+CTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATACT
+TACAGTAACATCCCACCCAACTTTGACTACTGGGGCCAGGGAGTCCTGGTCACCGTCTCC
+TCAG
+>query2
+CAGCTGCAGCTGCAGGAGTCGGGCCCAGGACTGGTGAAGCCTTCGGAGACCCTGTCCCTC
+ACCTGCGCTGTCTCTGTTGGCTCCATCAGCAGTAACTACTGGAGCTGGATCCGCCAGCCC
+CCAGGGAAGGGACTGGAGTGGATTGCACGTATCTCTGGTAGTGGTGGGAGCCCCGACTAC
+AACCCCTCCCGCAAGAGTCGAGTCACCATTTCAACAGACACGTCCAAGAGCCAGTTCTCC
+CTGAAGCTGAGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTATGCGAGAGATACT
+TACAGTACCATCCCACCCAACTTTGACTACTGGGGCCAGGGAGTCCTGGTCACCGTGTCC
+TCAG
+>query3
+CAGCTGCAGCTGCAGGAGTAGGGCCCAGGACTGGTGAAGCCTTCGGAGACGCTGTCCCTC
+ACCTGCGCTGTCTCTGGTGGCTCCATCAGCAGTAAGTACTGGAGCTGGATCCGCCAGCCC
+CCAGGGAATGGACTGGAGTGGATTGGACGTATCTCTGGTAGTGGTGGGAGCCCGGACTAC
+AACCCCTCCCTCAAGAGTCGAGTCACCATTTCAACAGACACGTCCATGAGCCAGTTCTCC
+CTGAAGCTGCGCTCTGTGACCGCCGCGGACACGGCCGTGTATTACTGTGCGAGAGATACT
+TACAGTAACATCCCACCCAACTTTGACTACTGGGGCCAGGGAGTCCTGGTCAGCGTCTCC
+TCAG