|
0
|
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
|
|
6
|
9 from chemfp import search
|
|
0
|
10 import sys
|
|
|
11 import os
|
|
4
|
12 import tempfile
|
|
6
|
13 import argparse
|
|
|
14 import subprocess
|
|
0
|
15
|
|
|
16
|
|
6
|
17 def unix_sort(results):
|
|
|
18 temp_unsorted = tempfile.NamedTemporaryFile(delete=False)
|
|
|
19 for (i,indices) in enumerate( results.iter_indices() ):
|
|
|
20 temp_unsorted.write('%s %s\n' % (len(indices), i))
|
|
|
21 print i, indices
|
|
|
22 temp_unsorted.close()
|
|
|
23 temp_sorted = tempfile.NamedTemporaryFile(delete=False)
|
|
|
24 temp_sorted.close()
|
|
|
25 p = subprocess.Popen(['sort', '-n', '-r', '-k', '1,1'], stdin=open(temp_unsorted.name), stdout=open(temp_sorted.name, 'w+'))
|
|
|
26 stdout, stderr = p.communicate()
|
|
|
27 return_code = p.returncode
|
|
0
|
28
|
|
6
|
29 if return_code:
|
|
|
30 sys.stdout.write(stdout)
|
|
|
31 sys.stderr.write(stderr)
|
|
|
32 sys.stderr.write("Return error code %i from command:\n" % return_code)
|
|
|
33 temp_sorted.close()
|
|
|
34 os.remove(temp_unsorted.name)
|
|
0
|
35
|
|
6
|
36 for line in open(temp_sorted.name):
|
|
|
37 size, fp_idx = line.strip().split()
|
|
|
38 yield (int(size), int(fp_idx))
|
|
|
39
|
|
|
40 os.remove(temp_sorted.name)
|
|
0
|
41
|
|
|
42
|
|
6
|
43 def butina( args ):
|
|
|
44 """
|
|
|
45 Taylor-Butina clustering from the chemfp help.
|
|
|
46 """
|
|
0
|
47
|
|
6
|
48 # make sure that the file ending is fps
|
|
|
49 temp_file = tempfile.NamedTemporaryFile()
|
|
|
50 temp_link = "%s.%s" % (temp_file.name, 'fps')
|
|
|
51 temp_file.close()
|
|
|
52 os.symlink(os.path.realpath(args.input_path), temp_link)
|
|
|
53 #os.system('ln -s %s %s' % (args.input_path, temp_link) )
|
|
|
54
|
|
|
55 out = args.output_path
|
|
|
56 arena = chemfp.load_fingerprints( temp_link )
|
|
0
|
57
|
|
6
|
58 chemfp.set_num_threads( args.processors )
|
|
|
59 results = search.threshold_tanimoto_search_symmetric(arena, threshold = args.tanimoto_threshold)
|
|
|
60 results.reorder_all("move-closest-first")
|
|
|
61 print [r.get_indices() for r in results]
|
|
0
|
62
|
|
6
|
63 # TODO: more memory efficient search?
|
|
|
64 # Reorder so the centroid with the most hits comes first.
|
|
|
65 # (That's why I do a reverse search.)
|
|
|
66 # Ignore the arbitrariness of breaking ties by fingerprint index
|
|
|
67 """
|
|
|
68 results = sorted( ( (len(indices), i, indices)
|
|
|
69 for (i,indices) in enumerate(results.iter_indices()) ),
|
|
|
70 reverse=True)
|
|
|
71 """
|
|
|
72 sorted_ids = unix_sort(results)
|
|
|
73
|
|
|
74 # Determine the true/false singletons and the clusters
|
|
|
75 true_singletons = []
|
|
|
76 false_singletons = []
|
|
|
77 clusters = []
|
|
0
|
78
|
|
6
|
79 seen = set()
|
|
|
80 #for (size, fp_idx, members) in results:
|
|
|
81 for (size, fp_idx) in sorted_ids:
|
|
|
82 members = results[fp_idx].get_indices()
|
|
|
83 print 'indices (s: %s, fp_idx:%s) -> %s' % (size, fp_idx, members)
|
|
|
84 print 'scores (s: %s, fp_idx:%s) -> %s' % (size, fp_idx, results[fp_idx].get_scores())
|
|
|
85 if fp_idx in seen:
|
|
|
86 # Can't use a centroid which is already assigned
|
|
|
87 continue
|
|
|
88 seen.add(fp_idx)
|
|
0
|
89
|
|
6
|
90 if size == 0:
|
|
|
91 # The only fingerprint in the exclusion sphere is itself
|
|
|
92 true_singletons.append(fp_idx)
|
|
|
93 continue
|
|
|
94
|
|
|
95 # Figure out which ones haven't yet been assigned
|
|
|
96 unassigned = set(members) - seen
|
|
|
97
|
|
|
98 if not unassigned:
|
|
|
99 false_singletons.append(fp_idx)
|
|
|
100 continue
|
|
|
101
|
|
|
102 # this is a new cluster
|
|
|
103 clusters.append( (fp_idx, unassigned) )
|
|
|
104 seen.update(unassigned)
|
|
0
|
105
|
|
6
|
106 len_cluster = len(clusters)
|
|
|
107 #out.write( "#%s true singletons: %s\n" % ( len(true_singletons), " ".join(sorted(arena.ids[idx] for idx in true_singletons)) ) )
|
|
|
108 #out.write( "#%s false singletons: %s\n" % ( len(false_singletons), " ".join(sorted(arena.ids[idx] for idx in false_singletons)) ) )
|
|
|
109
|
|
|
110 out.write( "#%s true singletons\n" % len(true_singletons) )
|
|
|
111 out.write( "#%s false singletons\n" % len(false_singletons) )
|
|
|
112 out.write( "#clusters: %s\n" % len_cluster )
|
|
0
|
113
|
|
6
|
114 # Sort so the cluster with the most compounds comes first,
|
|
|
115 # then by alphabetically smallest id
|
|
|
116 def cluster_sort_key(cluster):
|
|
|
117 centroid_idx, members = cluster
|
|
|
118 return -len(members), arena.ids[centroid_idx]
|
|
|
119
|
|
|
120 clusters.sort(key=cluster_sort_key)
|
|
0
|
121
|
|
6
|
122 for centroid_idx, members in clusters:
|
|
|
123 centroid_name = arena.ids[centroid_idx]
|
|
|
124 out.write("%s\t%s\t%s\n" % (centroid_name, len(members), " ".join(arena.ids[idx] for idx in members)))
|
|
|
125 #ToDo: len(members) need to be some biggest top 90% or something ...
|
|
|
126
|
|
|
127 for idx in true_singletons:
|
|
|
128 out.write("%s\t%s\n" % (arena.ids[idx], 0))
|
|
|
129
|
|
|
130 out.close()
|
|
|
131 os.remove( temp_link )
|
|
0
|
132
|
|
|
133
|
|
6
|
134 if __name__ == "__main__":
|
|
|
135 parser = argparse.ArgumentParser(description="""Taylor-Butina clustering for fps files.
|
|
|
136 For more details please see the original publication or the chemfp documentation:
|
|
|
137 http://www.chemomine.co.uk/dbclus-paper.pdf
|
|
|
138 https://chemfp.readthedocs.org
|
|
|
139 """)
|
|
|
140
|
|
|
141 parser.add_argument("-i", "--input", dest="input_path",
|
|
|
142 required=True,
|
|
|
143 help="Path to the input file.")
|
|
0
|
144
|
|
6
|
145 parser.add_argument("-o", "--output", dest="output_path", type=argparse.FileType('w'),
|
|
|
146 default=sys.stdout,
|
|
|
147 help="Path to the output file.")
|
|
|
148
|
|
|
149 parser.add_argument("-t", "--threshold", dest="tanimoto_threshold", type=float,
|
|
|
150 default=0.8,
|
|
|
151 help="Tanimoto threshold [0.8]")
|
|
|
152
|
|
|
153 parser.add_argument('-p', '--processors', type=int,
|
|
|
154 default=4)
|
|
|
155
|
|
|
156 options = parser.parse_args()
|
|
|
157 butina( options )
|
|
0
|
158
|
|
|
159
|
