1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
|
require 'digest/md5'
ENV["JAVA_HOME"] ||= "/usr/lib/jvm/java-7-openjdk"
BABEL_3D_CACHE_DIR = File.join(File.dirname(__FILE__),"..",'/babel_3d_cache')
# TODO store 3D structures in mongodb
# TODO store descriptors in mongodb
module OpenTox
module Algorithm
class Descriptor
include OpenTox
JAVA_DIR = File.join(File.dirname(__FILE__),"..","java")
CDK_JAR = Dir[File.join(JAVA_DIR,"cdk-*jar")].last
JOELIB_JAR = File.join(JAVA_DIR,"joelib2.jar")
LOG4J_JAR = File.join(JAVA_DIR,"log4j.jar")
JMOL_JAR = File.join(JAVA_DIR,"Jmol.jar")
obexclude = ["cansmi","cansmiNS","formula","InChI","InChIKey","s","smarts","title"]
OBDESCRIPTORS = Hash[OpenBabel::OBDescriptor.list_as_string("descriptors").split("\n").collect do |d|
name,description = d.split(/\s+/,2)
["Openbabel."+name,description] unless obexclude.include? name
end.compact.sort{|a,b| a[0] <=> b[0]}]
cdk_desc = YAML.load(`java -classpath #{CDK_JAR}:#{JAVA_DIR} CdkDescriptorInfo`)
CDKDESCRIPTORS = Hash[cdk_desc.collect { |d| ["Cdk."+d[:java_class].split('.').last.sub(/Descriptor/,''), d[:description]] }.sort{|a,b| a[0] <=> b[0]}]
CDKDESCRIPTOR_VALUES = cdk_desc.collect { |d| prefix="Cdk."+d[:java_class].split('.').last.sub(/Descriptor/,''); d[:names].collect{ |name| prefix+"."+name } }.flatten
# exclude Hashcode (not a physchem property) and GlobalTopologicalChargeIndex (Joelib bug)
joelibexclude = ["MoleculeHashcode","GlobalTopologicalChargeIndex"]
# strip Joelib messages from stdout
JOELIBDESCRIPTORS = Hash[YAML.load(`java -classpath #{JOELIB_JAR}:#{LOG4J_JAR}:#{JAVA_DIR} JoelibDescriptorInfo | sed '0,/---/d'`).collect do |d|
name = d[:java_class].sub(/^joelib2.feature.types./,'')
# impossible to obtain meaningful descriptions from JOELIb, see java/JoelibDescriptors.java
["Joelib."+name, "no description available"] unless joelibexclude.include? name
end.compact.sort{|a,b| a[0] <=> b[0]}]
DESCRIPTORS = OBDESCRIPTORS.merge(CDKDESCRIPTORS.merge(JOELIBDESCRIPTORS))
DESCRIPTOR_VALUES = OBDESCRIPTORS.keys + CDKDESCRIPTOR_VALUES + JOELIBDESCRIPTORS.keys
require_relative "unique_descriptors.rb"
def self.description descriptor
lib = descriptor.split('.').first
case lib
when "Openbabel"
OBDESCRIPTORS[descriptor]
when "Cdk"
name = descriptor.split('.')[0..-2].join('.')
CDKDESCRIPTORS[name]
when "Joelib"
JOELIBDESCRIPTORS[descriptor]
when "lookup"
"Read feature values from a dataset"
end
end
def self.smarts_match compounds, smarts, count=false
bad_request_error "Compounds for smarts_match are empty" unless compounds
bad_request_error "Smarts for smarts_match are empty" unless smarts
compounds = parse compounds
obconversion = OpenBabel::OBConversion.new
obmol = OpenBabel::OBMol.new
obconversion.set_in_format('inchi')
smarts_pattern = OpenBabel::OBSmartsPattern.new
#fingerprint = {}
smarts = [smarts] unless smarts.is_a? Array
fingerprint = Array.new(compounds.size){Array.new(smarts.size,false)}
compounds.each_with_index do |compound,c|
obconversion.read_string(obmol,compound.inchi)
smarts.each_with_index do |smart,s|
smarts_pattern.init(smart)
if smarts_pattern.match(obmol)
count ? value = smarts_pattern.get_map_list.to_a.size : value = 1
else
value = 0
end
fingerprint[c][s] = value
end
end
fingerprint
end
def self.smarts_count compounds, smarts
smarts_match compounds,smarts,true
end
def self.physchem compounds, descriptors=UNIQUEDESCRIPTORS
compounds = parse compounds
des = {}
descriptors.each do |d|
lib, descriptor = d.split(".",2)
lib = lib.downcase.to_sym
des[lib] ||= []
des[lib] << descriptor
end
result = {}
des.each do |lib,d|
send(lib, compounds, d).each do |compound,values|
result[compound] ||= {}
result[compound].merge! values
end
end
result
end
def self.openbabel compounds, descriptors
compounds = parse compounds
$logger.debug "compute #{descriptors.size} openbabel descriptors for #{compounds.size} compounds"
obdescriptors = descriptors.collect{|d| OpenBabel::OBDescriptor.find_type d}
obmol = OpenBabel::OBMol.new
obconversion = OpenBabel::OBConversion.new
obconversion.set_in_format 'inchi'
fingerprint = {}
compounds.each do |compound|
obconversion.read_string obmol, compound.inchi
fingerprint[compound] = {}
obdescriptors.each_with_index do |descriptor,i|
fingerprint[compound]["Openbabel."+descriptors[i]] = fix_value(descriptor.predict(obmol))
end
end
fingerprint
end
def self.cdk compounds, descriptors
compounds = parse compounds
$logger.debug "compute #{descriptors.size} cdk descriptors for #{compounds.size} compounds"
sdf = sdf_3d compounds
# use java system call (rjb blocks within tasks)
# use Tempfiles to avoid "Argument list too long" error
run_cmd "java -classpath #{CDK_JAR}:#{JAVA_DIR} CdkDescriptors #{sdf} #{descriptors.join(" ")}"
fingerprint = {}
YAML.load_file(sdf+"cdk.yaml").each_with_index do |calculation,i|
$logger.error "Descriptor calculation failed for compound #{compounds[i].inchi}." if calculation.empty?
descriptors.each do |descriptor|
fingerprint[compounds[i]] = calculation
end
end
FileUtils.rm sdf+"cdk.yaml"
fingerprint
end
def self.joelib compounds, descriptors
compounds = parse compounds
$logger.debug "compute #{descriptors.size} joelib descriptors for #{compounds.size} compounds"
# use java system call (rjb blocks within tasks)
# use Tempfiles to avoid "Argument list too long" error
sdf = sdf_3d compounds
run_cmd "java -classpath #{JOELIB_JAR}:#{JMOL_JAR}:#{LOG4J_JAR}:#{JAVA_DIR} JoelibDescriptors #{sdf} #{descriptors.join(' ')}"
fingerprint = {}
YAML.load_file(sdf+"joelib.yaml").each_with_index do |calculation,i|
$logger.error "Descriptor calculation failed for compound #{compounds[i].inchi}." if calculation.empty?
descriptors.each do |descriptor|
fingerprint[compounds[i]] = calculation
end
end
FileUtils.rm sdf+"joelib.yaml"
fingerprint
end
def self.lookup compounds, features, dataset
compounds = parse compounds
fingerprint = []
compounds.each do |compound|
fingerprint << []
features.each do |feature|
end
end
end
def self.run_cmd cmd
cmd = "#{cmd} 2>&1"
$logger.debug "running external cmd: '#{cmd}'"
p = IO.popen(cmd) do |io|
while line = io.gets
$logger.debug "> #{line.chomp}"
end
io.close
raise "external cmd failed '#{cmd}' (see log file for error msg)" unless $?.to_i == 0
end
end
def self.sdf_3d compounds
compounds = parse compounds
obconversion = OpenBabel::OBConversion.new
obmol = OpenBabel::OBMol.new
obconversion.set_in_format 'inchi'
obconversion.set_out_format 'sdf'
digest = Digest::MD5.hexdigest compounds.collect{|c| c.inchi}.inspect
sdf_file = "/tmp/#{digest}.sdf"
if File.exists? sdf_file # do not recreate existing 3d sdfs
$logger.debug "re-using cached 3d structures from #{sdf_file}"
else
tmp_file = Tempfile.new('sdf')
# create 3d sdf file (faster in Openbabel than in CDK)
# MG: moreover, CDK 3d generation is faulty
# MG: WARNING: Openbabel 3d generation is not deterministic
# MG: WARNING: Openbabel 3D generation does not work for mixtures
c = 0
compounds.each do |compound|
c += 1
cmp_file = File.join(BABEL_3D_CACHE_DIR,Digest::MD5.hexdigest(compound.inchi)+".sdf")
cmp_sdf = nil
if File.exists? cmp_file
$logger.debug "read cached 3d structure for compound #{c}/#{compounds.size}"
cmp_sdf = File.read(cmp_file)
else
$logger.debug "compute 3d structure for compound #{c}/#{compounds.size}"
obconversion.read_string obmol, compound.inchi
sdf_2d = obconversion.write_string(obmol)
error = nil
if compound.inchi.include?(";") # component includes multiple compounds (; in inchi, . in smiles)
error = "OpenBabel 3D generation failes for multi-compound #{compound.inchi}, trying to calculate descriptors from 2D structure."
else
OpenBabel::OBOp.find_type("Gen3D").do(obmol)
sdf_3d = obconversion.write_string(obmol)
error = "3D generation failed for compound #{compound.inchi}, trying to calculate descriptors from 2D structure." if sdf_3d.match(/.nan/)
end
if error
$logger.warn error
# TODO
# @feature_dataset[RDF::OT.Warnings] ? @feature_dataset[RDF::OT.Warnings] << error : @feature_dataset[RDF::OT.Warnings] = error
cmp_sdf = sdf_2d
else
cmp_sdf = sdf_3d
File.open(cmp_file,"w") do |f|
f.write(cmp_sdf)
end
end
end
tmp_file.write cmp_sdf
end
tmp_file.close
File.rename(tmp_file, sdf_file)
end
sdf_file
end
def self.parse compounds
case compounds.class.to_s
when "OpenTox::Compound"
compounds = [compounds]
when "Array"
compounds
when "OpenTox::Dataset"
compounds = compounds.compounds
else
bad_request_error "Cannot calculate descriptors for #{compounds.class} objects."
end
end
def self.fix_value val
val = val.first if val.is_a? Array and val.size == 1
if val.numeric?
val = Float(val)
val = nil if val.nan? or val.infinite?
else
val = nil if val == "NaN"
end
val
end
private_class_method :sdf_3d, :fix_value, :parse, :run_cmd
end
end
end
|