add chewtree tool: phylogenetic tree from chewBBACA allele profiles

tools/chewtree/.shed.yml

@@ -0,0 +1,11 @@
+name: chewtree
+owner: bgruening
+description: Builds phylogenetic trees from allelic profiles using ChewBBACA
+homepage_url: https://github.com/B-UMMI/ChewBBACA
+long_description: |
+  ChewTree is part of the ChewBBACA suite. It builds phylogenetic trees from allelic profiles computed by ChewBBACA's AlleleCall algorithm.
+remote_repository_url: https://github.com/bgruening/galaxytools/tree/master/tools/chewtree
+type: unrestricted
+categories:
+- Phylogenetics
+- Sequence Analysis

tools/chewtree/chewtree.xml

@@ -0,0 +1,75 @@
+<tool id="chewtree" name="chewTree" version="@TOOL_VERSION@+galaxy@VERSION_SUFFIX@" profile="@PROFILE@" license="MIT">
+    <description>Calculate a phylogenetic tree from chewBBACA allele profiles</description>
+    <macros>
+        <import>macros.xml</import>
+    </macros>
+    <expand macro="creator"/>
+    <expand macro="requirements"/>
+    <required_files>
+        <include path="scripts/mlst-hash-distance.py"/>
+        <include path="scripts/mentalist_tree.py"/>
+    </required_files>
+    <version_command>echo @TOOL_VERSION@</version_command>
+    <command detect_errors="aggressive"><![CDATA[
+        #if str($input_type) == "allele":
+            python '$__tool_directory__/scripts/mlst-hash-distance.py'
+                -i '$input'
+                -o '$phantcec_dm' &&
+        #else:
+            cp '$input' '$phantcec_dm' &&
+        #end if
+        python '$__tool_directory__/scripts/mentalist_tree.py' '$phantcec_dm' > '$phantcec_tree'
+    ]]></command>
+    <inputs>
+        <param name="input" format="tabular" type="data" label="chewBBACA allele profile TSV file" help="Tab-separated file with allele calls; first column sample names, subsequent columns loci, last two columns (CC, clonal_complex) excluded"/>
+        <param name="input_type" type="select" label="Input matrix type">
+            <option value="allele" selected="true">Allele profile - calculate Hamming distance matrix</option>
+            <option value="distance">Distance matrix - use pre-computed distance matrix directly</option>
+        </param>
+    </inputs>
+    <outputs>
+        <data format="tabular" name="phantcec_dm" label="${tool.name} on ${on_string}: distance matrix"/>
+        <data format="txt" name="phantcec_tree" label="${tool.name} on ${on_string}: phylogenetic tree (Newick)"/>
+    </outputs>
+    <tests>
+        <test expect_num_outputs="2">
+            <param name="input" value="allele_profile.tsv" ftype="tabular"/>
+            <param name="input_type" value="allele"/>
+            <output name="phantcec_dm" ftype="tabular">
+                <assert_contents>
+                    <has_text text="sample1"/>
+                    <has_text text="sample2"/>
+                    <has_text text="sample3"/>
+                </assert_contents>
+            </output>
+            <output name="phantcec_tree" ftype="txt">
+                <assert_contents>
+                    <has_text text=";"/>
+                </assert_contents>
+            </output>
+        </test>
+    </tests>
+    <help format="markdown"><![CDATA[
+**chewTree** - Build a phylogenetic tree from chewBBACA allele profiles
+
+-----
+
+This tool calculates a Hamming distance matrix from a chewBBACA allele profile,
+then constructs a neighbor-joining tree using Bio.Phylo.
+
+**Input**
+
+- **TSV file**: chewBBACA allele profile with samples as rows and loci as columns.
+  The last two columns (CC, clonal_complex) are ignored automatically.
+
+- **Input type**:
+  - *Allele profile*: Calculate pairwise Hamming distances directly from allele calls.
+  - *Distance matrix*: Use a pre-computed distance matrix to build the tree.
+
+**Output**
+
+- **Distance matrix** (tabular): Pairwise Hamming distance matrix.
+- **Phylogenetic tree** (Newick): Neighbor-joining tree built from the distance matrix.
+    ]]></help>
+    <expand macro="citations"/>
+</tool>

tools/chewtree/macros.xml

@@ -0,0 +1,23 @@
+<macros>
+    <token name="@TOOL_VERSION@">1.0</token>
+    <token name="@VERSION_SUFFIX@">0</token>
+    <token name="@PROFILE@">25.0</token>
+    <xml name="creator">
+        <creator>
+            <person givenName="Paul" familyName="Zierep" url="https://github.com/paulzierep" identifier="https://orcid.org/0000-0003-2982-388X"/>
+            <person givenName="Björn" familyName="Grüning" url="https://github.com/bgruening" identifier="https://orcid.org/0000-0002-3079-6586"/>
+        </creator>
+    </xml>
+    <xml name="requirements">
+        <requirements>
+            <requirement type="package" version="2.5.2">numpy</requirement>
+            <requirement type="package" version="1.88">biopython</requirement>
+        </requirements>
+    </xml>
+    <xml name="citations">
+        <citations>
+            <citation type="doi">10.1099/mgen.0.000166</citation>
+            <yield/>
+        </citations>
+    </xml>
+</macros>

tools/chewtree/scripts/mentalist_tree.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+
+import sys
+import csv
+import numpy as np
+
+import Bio.Phylo
+from Bio.Phylo.TreeConstruction import DistanceMatrix, DistanceTreeConstructor
+
+def usage():
+  print("usage: mentalist_tree <input.tsv>\n")
+
+def process_input_matrix(input_matrix):
+    """ Converts an array-of-arrays containting sample IDs and distances
+        into a BioPython DistanceMatrix object
+    """
+    input_matrix.pop(0)
+    sample_names = [row[0] for row in input_matrix]
+    for row in input_matrix:
+        row.pop(0)
+    distance_matrix = []
+    for input_matrix_row in input_matrix:
+        distance_matrix.append([int(i) for i in input_matrix_row])
+    """ np.tril() converts a matrix like this: [[0 1 2]
+                                                [1 0 1]
+                                                [2 1 0]]
+        ...into this: [[0 0 0]
+                       [1 0 0]
+                       [2 1 0]]
+        ...but what we need to pass to DistanceMatrix() is this: [[0]
+                                                                  [1 0]
+                                                                  [2 1 0]]
+        ...so that's what the (somewhat cryptic) code below does.
+    """
+    distance_matrix = np.tril(np.array(distance_matrix))
+    num_rows = distance_matrix.shape[0]
+    """ masking the distance matrix with tril_indices gives a linearized
+        distance matrix [0 1 0 2 1 0] that we need to re-construct into [[0], [1, 0], [2, 1, 0]]
+    """
+    lower_triangular_idx_mask = np.tril_indices(num_rows)
+    linear_distance_matrix = distance_matrix[lower_triangular_idx_mask]
+    distance_matrix = []
+    min = 0
+    max = 1
+    for i in range(num_rows):
+        distance_matrix.append(linear_distance_matrix[min:max].tolist())
+        min = max
+        max = max + (i + 2) 
+    distance_matrix = DistanceMatrix(names=sample_names, matrix=distance_matrix)
+    return distance_matrix
+
+def main():
+    if len(sys.argv) < 2:
+        usage()
+        sys.exit(1)
+
+    input_file = sys.argv[1]
+    reader = csv.reader(open(input_file, "r"), delimiter="\t")
+    input_matrix = list(reader)
+    # Don't build a tree with fewer than 3 samples, just produce an empty file
+    if len(input_matrix) < 4:
+      print('();')
+      sys.exit(0)
+    distance_matrix = process_input_matrix(input_matrix)
+    constructor = DistanceTreeConstructor()
+    tree = constructor.nj(distance_matrix)
+    Bio.Phylo.write(tree, sys.stdout, 'newick')
+
+if __name__ == '__main__':
+    main()

tools/chewtree/scripts/mlst-hash-distance.py

@@ -0,0 +1,114 @@
+# flake8: noqa
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+import sys, getopt
+import os
+
+
+def mlst_calls(call_file):
+  # returns an MLST call matrix, with samples on the rows and loci on the columns.
+  # header row is excluded and last two columns (CC, clonal_complex) are excluded
+  with open(call_file) as file_in:
+    lines = []
+    for line in file_in:
+      lines.append(line.split('\t')[:-2])
+  data_mlst = lines[1:]
+  return data_mlst
+
+
+def compare_alleles(allele_1, allele_2):
+  # compare alleles, return 1 if different, 0 if equal or not to be compared
+  if allele_1 == allele_2:
+    comparison = 0
+  elif allele_1 == 0 or allele_2 == 0:
+    # allele not found
+    comparison = 0
+  elif allele_1 == '0' or allele_2 == '0':
+    # allele not found
+    comparison = 0
+  elif allele_1 == 'N' or allele_2 == 'N':
+    # new allele inferred
+    comparison = 0
+  elif '+' in allele_1 or '+' in allele_2:
+    # partial allele found
+    comparison = 0
+  elif 'INF' in allele_1 or 'INF' in allele_2:
+    # new allele inferred
+    comparison = 0
+  elif 'LNF' in allele_1 or 'LNF' in allele_2:
+    # allele not found
+    comparison = 0
+  elif 'PLOT' in allele_1 or 'PLOT' in allele_2:
+    # partial allele found
+    comparison = 0
+  elif 'NIPH' in allele_1 or 'NIPH' in allele_2:
+    # partial allele found
+    comparison = 0
+  elif 'ALM' in allele_1 or 'ALM' in allele_2:
+    # partial allele found
+    comparison = 0
+  elif 'ASM' in allele_1 or 'ASM' in allele_2:
+    # partial allele found
+    comparison = 0
+  else:
+    comparison = 1
+  return comparison
+
+
+def mlst_distance(mlst):
+  # a profile file was given, substitute the allele numbers with the allele sequence hashes
+  rows = len(mlst)
+  cols = len(mlst[0])
+  D = [ [0]*(rows) for _ in range(rows) ]
+  h = []
+  for row in range(0, rows):
+    h.append(mlst[row][0])
+    for row2 in range(row+1, rows):
+      dist = 0
+      for col in range(1, cols):
+        dist = dist + compare_alleles(mlst[row][col], mlst[row2][col])
+      D[row][row2] = dist
+      D[row2][row] = dist
+  return D, h
+
+
+def main(argv):
+  input = ''
+  strusage = 'mlst-hash-distance.py -i <input.tsv> -o <output.tsv>\n'
+  numloci = 0
+  try:
+    opts, args = getopt.getopt(argv,"hi:o:",["input=","output="])
+  except getopt.GetoptError:
+    print (strusage)
+    sys.exit(2)
+  for opt, arg in opts:
+    if opt == '-h':
+      print (strusage)
+      sys.exit()
+    elif opt in ("-i", "--input"):
+      input = arg
+    elif opt in ("-o", "--output"):
+      output = arg
+  if os.path.isfile(input):
+    print ('input file is "', input, '"')
+  else: 
+    print ('input file is "', input, '" but does not exist')
+    sys.exit(0)
+  print ('output file is "', output, '"')
+  mlst = mlst_calls(input)
+  dist_mat, head_mat = mlst_distance(mlst)
+  header = '\t'.join(head_mat)
+  with open(output, "w") as output:
+    output.write('\t' + header + '\n')
+    i = 0
+    for row in dist_mat:
+      output.write(head_mat[i])
+      i = i + 1
+      for elem in row:
+        output.write('\t' + str(elem))
+      output.write('\n')
+
+
+if __name__ == "__main__":
+   main(sys.argv[1:])
+

tools/chewtree/test-data/allele_profile.tsv

@@ -0,0 +1,4 @@
+	 locus1	 locus2	 CC	 clonal_complex
+sample1	 1	 1	 1	 1
+sample2	 1	 2	 1	 1
+sample3	 2	 2	 1	 1