diff --git a/tools/chewtree/.shed.yml b/tools/chewtree/.shed.yml new file mode 100644 index 0000000000..d900e19334 --- /dev/null +++ b/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 diff --git a/tools/chewtree/chewtree.xml b/tools/chewtree/chewtree.xml new file mode 100644 index 0000000000..b373d59ce5 --- /dev/null +++ b/tools/chewtree/chewtree.xml @@ -0,0 +1,75 @@ + + Calculate a phylogenetic tree from chewBBACA allele profiles + + macros.xml + + + + + + + + echo @TOOL_VERSION@ + '$phantcec_tree' + ]]> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tools/chewtree/macros.xml b/tools/chewtree/macros.xml new file mode 100644 index 0000000000..596f04851e --- /dev/null +++ b/tools/chewtree/macros.xml @@ -0,0 +1,23 @@ + + 1.0 + 0 + 25.0 + + + + + + + + + numpy + biopython + + + + + 10.1099/mgen.0.000166 + + + + diff --git a/tools/chewtree/scripts/mentalist_tree.py b/tools/chewtree/scripts/mentalist_tree.py new file mode 100755 index 0000000000..15d6910522 --- /dev/null +++ b/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 \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() diff --git a/tools/chewtree/scripts/mlst-hash-distance.py b/tools/chewtree/scripts/mlst-hash-distance.py new file mode 100755 index 0000000000..85af316bba --- /dev/null +++ b/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 -o \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:]) + diff --git a/tools/chewtree/test-data/allele_profile.tsv b/tools/chewtree/test-data/allele_profile.tsv new file mode 100644 index 0000000000..60391bef42 --- /dev/null +++ b/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