Mercurial > repos > bgruening > chemfp
comparison butina_clustering.py @ 2:70b071de9bee draft
planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/chemicaltoolbox/chemfp commit 01da22e4184a5a6f6a3dd4631a7b9c31d1b6d502
author | bgruening |
---|---|
date | Sat, 20 May 2017 08:31:44 -0400 |
parents | |
children | 3b14765c22ee |
comparison
equal
deleted
inserted
replaced
1:43a9e7d9b24f | 2:70b071de9bee |
---|---|
1 #!/usr/bin/env python | |
2 """ | |
3 Modified version of code examples from the chemfp project. | |
4 http://code.google.com/p/chem-fingerprints/ | |
5 Thanks to Andrew Dalke of Andrew Dalke Scientific! | |
6 """ | |
7 | |
8 import chemfp | |
9 import sys | |
10 import os | |
11 import tempfile | |
12 import argparse | |
13 import subprocess | |
14 from chemfp import search | |
15 | |
16 def unix_sort(results): | |
17 temp_unsorted = tempfile.NamedTemporaryFile(delete=False) | |
18 for (i,indices) in enumerate( results.iter_indices() ): | |
19 temp_unsorted.write('%s %s\n' % (len(indices), i)) | |
20 temp_unsorted.close() | |
21 temp_sorted = tempfile.NamedTemporaryFile(delete=False) | |
22 temp_sorted.close() | |
23 p = subprocess.Popen(['sort', '-n', '-r', '-k', '1,1'], stdin=open(temp_unsorted.name), stdout=open(temp_sorted.name, 'w+')) | |
24 stdout, stderr = p.communicate() | |
25 return_code = p.returncode | |
26 | |
27 if return_code: | |
28 sys.stdout.write(stdout) | |
29 sys.stderr.write(stderr) | |
30 sys.stderr.write("Return error code %i from command:\n" % return_code) | |
31 temp_sorted.close() | |
32 os.remove(temp_unsorted.name) | |
33 | |
34 for line in open(temp_sorted.name): | |
35 size, fp_idx = line.strip().split() | |
36 yield (int(size), int(fp_idx)) | |
37 | |
38 os.remove(temp_sorted.name) | |
39 | |
40 def butina( args ): | |
41 """ | |
42 Taylor-Butina clustering from the chemfp help. | |
43 """ | |
44 out = args.output_path | |
45 targets = chemfp.open( args.input_path, format='fps' ) | |
46 arena = chemfp.load_fingerprints( targets ) | |
47 | |
48 chemfp.set_num_threads( args.processors ) | |
49 results = search.threshold_tanimoto_search_symmetric(arena, threshold = args.tanimoto_threshold) | |
50 results.reorder_all("move-closest-first") | |
51 | |
52 sorted_ids = unix_sort(results) | |
53 | |
54 # Determine the true/false singletons and the clusters | |
55 true_singletons = [] | |
56 false_singletons = [] | |
57 clusters = [] | |
58 | |
59 seen = set() | |
60 #for (size, fp_idx, members) in results: | |
61 for (size, fp_idx) in sorted_ids: | |
62 members = results[fp_idx].get_indices() | |
63 #print arena.ids[ fp_idx ], [arena.ids[ m ] for m in members] | |
64 if fp_idx in seen: | |
65 # Can't use a centroid which is already assigned | |
66 continue | |
67 seen.add(fp_idx) | |
68 | |
69 if size == 0: | |
70 # The only fingerprint in the exclusion sphere is itself | |
71 true_singletons.append( fp_idx ) | |
72 continue | |
73 | |
74 # Figure out which ones haven't yet been assigned | |
75 unassigned = set(members) - seen | |
76 | |
77 if not unassigned: | |
78 false_singletons.append(fp_idx) | |
79 continue | |
80 | |
81 # this is a new cluster | |
82 clusters.append( (fp_idx, unassigned) ) | |
83 seen.update(unassigned) | |
84 | |
85 len_cluster = len(clusters) | |
86 #out.write( "#%s true singletons: %s\n" % ( len(true_singletons), " ".join(sorted(arena.ids[idx] for idx in true_singletons)) ) ) | |
87 #out.write( "#%s false singletons: %s\n" % ( len(false_singletons), " ".join(sorted(arena.ids[idx] for idx in false_singletons)) ) ) | |
88 | |
89 out.write( "#%s true singletons\n" % len(true_singletons) ) | |
90 out.write( "#%s false singletons\n" % len(false_singletons) ) | |
91 out.write( "#clusters: %s\n" % len_cluster ) | |
92 | |
93 # Sort so the cluster with the most compounds comes first, | |
94 # then by alphabetically smallest id | |
95 def cluster_sort_key(cluster): | |
96 centroid_idx, members = cluster | |
97 return -len(members), arena.ids[centroid_idx] | |
98 | |
99 clusters.sort(key=cluster_sort_key) | |
100 | |
101 for centroid_idx, members in clusters: | |
102 centroid_name = arena.ids[centroid_idx] | |
103 out.write("%s\t%s\t%s\n" % (centroid_name, len(members), " ".join(arena.ids[idx] for idx in members))) | |
104 #ToDo: len(members) need to be some biggest top 90% or something ... | |
105 | |
106 for idx in true_singletons: | |
107 out.write("%s\t%s\n" % (arena.ids[idx], 0)) | |
108 | |
109 out.close() | |
110 | |
111 | |
112 if __name__ == "__main__": | |
113 parser = argparse.ArgumentParser(description="""Taylor-Butina clustering for fps files. | |
114 For more details please see the original publication or the chemfp documentation: | |
115 http://www.chemomine.co.uk/dbclus-paper.pdf | |
116 https://chemfp.readthedocs.org | |
117 """) | |
118 | |
119 parser.add_argument("-i", "--input", dest="input_path", | |
120 required=True, | |
121 help="Path to the input file.") | |
122 | |
123 parser.add_argument("-o", "--output", dest="output_path", type=argparse.FileType('w'), | |
124 default=sys.stdout, | |
125 help="Path to the output file.") | |
126 | |
127 parser.add_argument("-t", "--threshold", dest="tanimoto_threshold", type=float, | |
128 default=0.8, | |
129 help="Tanimoto threshold [0.8]") | |
130 | |
131 parser.add_argument('-p', '--processors', type=int, default=4) | |
132 | |
133 options = parser.parse_args() | |
134 butina( options ) |