Mercurial > repos > iuc > tn93_filter

comparison tn93_cluster.py @ 0:ba95715078c9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/tn93/ commit eec640a7c26b728f8175885926fe368b0756d9e5"
author iuc
date Fri, 23 Apr 2021 03:05:08 +0000
parents
children cf50aeb956f2
comparison
equal deleted inserted replaced
-1:000000000000 0:ba95715078c9
1 import argparse
2 import json
3 import os
4 import shlex
5 import shutil
6 import subprocess
7 import sys
8
9
10 def cluster_to_fasta(json_file, fasta_file, reference_name=None):
11 with open(json_file, "r") as fh:
12 cluster_json = json.load(fh)
13 with open(fasta_file, "w") as fh2:
14 for c in cluster_json:
15 if reference_name is not None:
16 if reference_name in c['members']:
17 cc = c['centroid'].split('\n')
18 cc[0] = ">" + reference_name
19 print("\n".join(cc), file=fh2)
20 continue
21 print(c['centroid'], file=fh2)
22
23 return(os.path.getmtime(fasta_file), len(cluster_json))
24
25
26 def run_command(command):
27 proc = subprocess.Popen(shlex.split(command), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
28 stdout, stderr = proc.communicate()
29 result = proc.returncode
30 if result != 0:
31 print('Command `%s` failed with exit code %s\n' % (command, result), file=sys.stderr)
32 print('--------------------- STDOUT ---------------------')
33 print(stdout.decode().replace('\\n', '\n'))
34 print('------------------- END STDOUT -------------------')
35 print('--------------------- STDERR ---------------------', file=sys.stderr)
36 print(stderr.decode().replace('\\n', '\n'), file=sys.stderr)
37 print('------------------- END STDERR -------------------', file=sys.stderr)
38 return(int(result))
39
40
41 def main(arguments):
42 threshold = arguments.threshold
43 step = threshold * 0.25
44 shutil.copy(arguments.input, os.path.join(os.getcwd(), 'reference_msa.fa'))
45 shutil.copy(arguments.input, os.path.join(os.getcwd(), 'reference_msa.fa.bak'))
46 with open(arguments.reference) as fh:
47 for line in fh:
48 if line[0] == '>':
49 _ref_seq_name = line[1:].split(' ')[0].strip()
50 break
51 while True and threshold <= 1:
52 command = 'tn93-cluster -o clusters.json -t %g -a %s -c %s -m json -l %d -g %f reference_msa.fa' % (threshold, arguments.ambigs, arguments.cluster_type, arguments.overlap, arguments.fraction)
53 return_code = run_command(command)
54 if return_code != 0:
55 return return_code
56 input_stamp, cluster_count = cluster_to_fasta('clusters.json', 'reference_msa.fa.bak', _ref_seq_name)
57 if cluster_count <= arguments.cluster_count or threshold == 1:
58 break
59 else:
60 threshold += step
61 print('Found %d clusters at threshold %f' % (cluster_count, threshold))
62 shutil.copy('reference_msa.fa.bak', arguments.compressed)
63 shutil.copy('clusters.json', arguments.output)
64 os.remove('reference_msa.fa.bak')
65 return 0
66
67
68 if __name__ == '__main__':
69 parser = argparse.ArgumentParser(description='Combine alignments into a single file, adding a reference sequence as well')
70 parser.add_argument('--input', help='Input MSA', required=True, type=str)
71 parser.add_argument('--reference', help='Reference sequence', required=True, type=str)
72 parser.add_argument('--output', help='Input MSA', required=True, type=str)
73 parser.add_argument('--threshold', help='Threshold', required=True, type=float)
74 parser.add_argument('--ambigs', help='Handle ambigs', required=True, type=str)
75 parser.add_argument('--cluster-type', help='Cluster type', required=True, type=str)
76 parser.add_argument('--overlap', help='Overlap', required=True, type=int)
77 parser.add_argument('--fraction', help='Fraction', required=True, type=float)
78 parser.add_argument('--cluster-count', help='Max query', required=True, type=int)
79 parser.add_argument('--compressed', help='File to write compressed clusters to', required=True, type=str)
80 arguments = parser.parse_args()
81 exit(main(arguments))