Mercurial > repos > bgruening > chemfp
annotate chemfp_clustering/butina_clustering.py @ 28:f874af3cf8df
ChemicalToolBoX update.
| author | Bjoern Gruening <bjoern.gruening@gmail.com> |
|---|---|
| date | Tue, 23 Jul 2013 00:54:07 +0200 |
| parents | 6c496b524b41 |
| children |
| rev | line source |
|---|---|
| 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 | |
| 9 import sys | |
| 10 import os | |
| 4 | 11 import tempfile |
| 6 | 12 import argparse |
| 13 import subprocess | |
| 11 | 14 from chemfp import search |
| 0 | 15 |
| 6 | 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 | |
| 0 | 26 |
| 6 | 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) | |
| 0 | 33 |
| 6 | 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) | |
| 0 | 39 |
| 6 | 40 def butina( args ): |
| 41 """ | |
| 42 Taylor-Butina clustering from the chemfp help. | |
| 43 """ | |
| 44 out = args.output_path | |
|
22
6c496b524b41
ChemicalToolBoX update.
Bjoern Gruening <bjoern.gruening@gmail.com>
parents:
11
diff
changeset
|
45 targets = chemfp.open( args.input_path, format='fps' ) |
|
6c496b524b41
ChemicalToolBoX update.
Bjoern Gruening <bjoern.gruening@gmail.com>
parents:
11
diff
changeset
|
46 arena = chemfp.load_fingerprints( targets ) |
| 0 | 47 |
| 6 | 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") | |
| 0 | 51 |
| 6 | 52 sorted_ids = unix_sort(results) |
| 11 | 53 |
| 6 | 54 # Determine the true/false singletons and the clusters |
| 55 true_singletons = [] | |
| 56 false_singletons = [] | |
| 57 clusters = [] | |
| 0 | 58 |
| 6 | 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() | |
| 11 | 63 #print arena.ids[ fp_idx ], [arena.ids[ m ] for m in members] |
| 6 | 64 if fp_idx in seen: |
| 65 # Can't use a centroid which is already assigned | |
| 66 continue | |
| 67 seen.add(fp_idx) | |
| 0 | 68 |
| 6 | 69 if size == 0: |
| 70 # The only fingerprint in the exclusion sphere is itself | |
| 11 | 71 true_singletons.append( fp_idx ) |
| 6 | 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) | |
| 0 | 84 |
| 6 | 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 ) | |
| 0 | 92 |
| 6 | 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) | |
| 0 | 100 |
| 6 | 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() | |
| 0 | 110 |
| 111 | |
| 6 | 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.") | |
| 0 | 122 |
| 6 | 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, | |
| 132 default=4) | |
| 133 | |
| 134 options = parser.parse_args() | |
| 135 butina( options ) |